https://orcid.org/0000-0001-6365-6849
ABSTRACT
Real-time flood forecasting requires computationally efficient models, and implicit flood models are often preferred due to their ability to use larger time steps, enhancing computational efficiency. However, these models may suffer from spurious checkerboard oscillations, causing unrealistic predictions of water depth and velocity. While some models employ the Rhie and Chow interpolation technique to mitigate these oscillations, this approach depends on the under-relaxation factor, affecting convergence and computation time. This study introduces IITMFLO-2D, an implicit 2D flood model developed using the SIMPLE-Consistent (SIMPLEC) algorithm on a non-staggered grid system, with two novel modifications to reduce checkerboard oscillations. The first modification uses an explicit equation, similar to Roe-type averaging, for the momentum weighting factor based on local flow depths. The second introduces a Courant–Friedrichs–Lewy (CFL)-based weighting factor to determine interfacial water levels, adapting to local flow dynamics. These adaptive weighting factors vary spatially and temporally. Performance of IITMFLO-2D is verified through analytical test cases, hypothetical flood simulations, and experimental events, showing significant reductions in computation time and improved accuracy compared to other models such as HECRAS-2D and TELEMAC-2D. These innovations make IITMFLO-2D suitable for real-time flood forecasting and applicable to other numerical models with complex meshes.
HIGHLIGHTS
Development of IITMFLO-2D, an efficient flood forecasting model using SIMPLEC algorithm on a non-staggered grid system.
Introduction of adaptive weighting factors to compute fluxes and interfacial water levels, reducing dependency on under-relaxation factors.
Verification through analytical tests and real-world scenarios, demonstrating improved computational efficiency.
INTRODUCTION
Over the years, researchers have developed various hydraulic models to simulate surface water flows, ranging from one-dimensional (1D) to three-dimensional (3D) models. Among these, two-dimensional (2D) flow models are extensively used to simulate free surface flows in environments such as open channels, rivers, floodplains, and coastal estuaries (Daoud et al. 2008; Villanueva et al. 2008; Wu & Lin 2015; Shipton et al. 2018). These 2D models are based on the shallow water equations (SWEs), derived by integrating the 3D Navier–Stokes equations over the flow depth. Several numerical methods have been developed to solve the SWEs, including the finite difference method (FDM) (Fletcher 1991), the finite element method (FEM) (Chung 1978; Zienkiewicz et al. 2013), and the finite-volume method (FVM) (Patankar 1980; Ferziger & Perić 2002; Wu 2007). Of these, FVM has become the most popular and widely used due to its better mass and momentum conservation properties and its ability to effectively capture flow discontinuities (Leighton 2005; Bauer & Cotter 2018; Shipton et al. 2018).
The FVM is typically divided into two categories: explicit and implicit. In the explicit FVM, the SWEs are solved directly for each cell, while the implicit FVM solves the SWEs simultaneously across the domain. Explicit methods often utilize shock-capturing schemes such as Roe (1981), HLL (Harten, Lax, and van Leer) (Toro et al. 1994), and HLLC (HLL Central wave) (Batten et al. 1997) Riemann solvers, which have been widely recognized for their effectiveness in handling flow discontinuities (Leon et al. 2008; Marras et al. 2018). These schemes are well-suited for modeling rapid flow changes, such as dam breaks, levee breaches, and flash floods. However, they are computationally intensive, requiring more computations per time step and being limited-by the Courant–Friedrichs–Lewy (CFL) condition, which constrains the time step size and increases overall runtime (Tavelli & Dumbser 2017). On the other hand, the implicit FVM 2D models offer a more efficient approach for simulating slow-rising floods, such as those resulting from heavy rainfall over large basins. By allowing larger time steps, these models reduce computation time while still providing adequate lead time for emergency planning (Daoud et al. 2008; Darwish et al. 2009). As a result, they are ideal for situations where efficient flood prediction is needed.
In the implicit solution of the SWEs, the momentum equations link velocity to the pressure gradient, and the continuity equation serves solely as an additional constraint on the velocity field without a direct connection to pressure. This weak linkage significantly influences the convergence and stability of numerical solutions, making the coupling between pressure and velocity fields crucial. When variables are stored at the geometric center of control volumes and linear interpolation is applied for internodal variation, it can result in nonphysical checkerboard oscillations (Walters & Carey 1983; Ferziger & Perić 2002; Javan et al. 2007; Versteeg & Malalasekera 2007; Larmaei et al. 2010). To mitigate these oscillations, one effective method is to employ a staggered grid. Then the coupling of pressure and velocity (i.e., SWEs) can be solved using iterative algorithms such as SIMPLE (Semi-Implicit Method for Pressure-Linked Equation) (Patankar & Spalding 1970), SIMPLER (SIMPLE-Revised) (Patankar 1981), SIMPLEC (SIMPLE-Consistent) (Van Doormaal & Raithby 1984), PISO (Pressure Implicit with Splitting of Operators) (Issa 1986) and IDEAL (Inner Doubly-iterative Efficient Algorithm for Linked-equations) (Sun et al. 2008). The higher computational cost of these pressure–velocity coupling was earlier solved by assuming a hydrostatic pressure distribution, which, however, fails to address flows over complex bed topographies. Non-hydrostatic pressure distribution was hence adopted, introducing the z-coordinate approach by Casulli & Zanolli (2002). This was further modified using the sigma coordinate approach by Hu et al. (2013).
The staggered grid arrangement stores primitive variables in different control volumes, naturally preventing checkerboard oscillations. This arrangement generally offers better stability and robustness. However, it poses challenges when implementing solution algorithms on unstructured and non-orthogonal grids, often needed to represent irregular topographical features accurately. The staggered grid also requires more memory storage, making it computationally expensive (Zijlema et al. 1995; Daoud et al. 2008; Walters et al. 2009; Zou et al. 2016). Alternatively, the momentum interpolation technique proposed by Rhie & Chow (1983) can be applied on a non-staggered grid, allowing for the previously mentioned algorithms to be employed for the coupling of velocity and pressure. The Rhie–Chow interpolation (momentum or pressure-weighted interpolation) method involves calculating an intermediate velocity that excludes the pressure effect and then interpolating this velocity from cell centers to cell faces. The false oscillations are hence eliminated by this refinement of the interpolation of the pressure gradient from neighboring grid nodes. Incorporating the Rhie and Chow interpolation, Javan et al. (2007) developed a time-splitting technique to simulate 2D non-hydrostatic free surface flow in a vertical plane.
The SIMPLEC algorithm, when combined with the Rhie and Chow interpolation technique in an implicit manner, offers significant computational advantages over semi-implicit approaches or time-splitting methods, which can help mitigate oscillations. Darwish et al. (2009) demonstrated the effectiveness of this combination in solving SWEs, highlighting the influence of the under-relaxation factor on both convergence rates and solution accuracy. Many studies have adopted the Rhie and Chow interpolation technique with the SIMPLEC algorithm to solve free surface flows (Shen et al. 2001; Yu et al. 2002; Wu et al. 2011; Wu & Lin 2015; Moguen et al. 2019). Wu et al. (2011) further refined this approach by identifying an optimal under-relaxation factor of 0.8 through trial and error. The under-relaxation factor is typically a numerical approach to enhance the convergence and stabilize the solutions of iteratively solved non-linear equations. It defines the fraction to which the iterative solution is updated with the computed values. Despite successfully utilizing the Rhie and Chow method for flux computation, challenges remain, particularly concerning oscillations arising from simple averaging of flow variables for flux computation. The use of momentum weighting factors, typically fixed at 0.5 for regular grids, fails to account for local flow dynamics, which can undermine the stability and efficiency of the numerical scheme. Furthermore, it is crucial to avoid conventional methods such as simple linear interpolation or central-difference techniques when determining water surface levels at cell interfaces, as they neglect local flow dynamics and compromise model stability and accuracy. Therefore, this study proposes two novel ideas – (i) the use of adaptive weights based on Roe’s (1981) variable averaging method for velocity and (ii) a Courant number-based weighting factor for estimating water surface elevation at cell interfaces. These weighting factors are integrated into the Rhie and Chow interpolation technique within a non-staggered grid framework, aiming to mitigate nonphysical solutions and ultimately improve model stability, accuracy, and computational efficiency.
This research thus presents a novel 2D numerical model that employs regular, non-staggered computational grids, facilitating the direct use of digital elevation models (DEMs) for effective flood simulations. Subsequent sections provide a comprehensive formulation of the model, outlining the governing equations, the numerical discretization schemes utilized, the solutions to the discretized equations, and the methods for coupling velocity with water levels. The enhancements of the developed model are compared with those in Wu et al. (2011) and Wu & Lin (2015) by evaluating a diverse set of test cases, including laboratory experiments conducted in a river-network-floodplain system (Mali & Kuiry 2020; Mali et al. 2020).
MODEL FORMULATION
The detailed approach of the model formulation emphasizing the improvements over the existing solutions is briefed in the subsections.
Governing equations for 2D shallow water flow
The depth-averaged 2D SWEs consist of continuity and momentum equations in the x- and y- directions. The governing equations in the Cartesian coordinate system are expressed as in Cunge (1980):
is the water surface elevation above the reference level (m), in which h is the local water depth (m), and 
is the bed elevation from a horizontal datum (m); 
is the density of water (kg/m3); 
and 
are shear stresses (N/m2) at the interface of bed and flow in the x- and y-directions, respectively and are determined as 
, 
where, 
. The n = Manning’s coefficient (m-1/3s) is an efficient approach to represent bed shear stress in which the value of n depends on the local bed roughness.Numerical discretization
(m) and 
(m) in the x and y directions, respectively. If the grids are exactly square like in a DEM grid, then 
. 
Non-staggered grid arrangement over a regular computational mesh.
define the current and next-time steps during the computation. The Green’s theorem is then applied to the second term of Equation (4) to yield the discretized form of the continuity equation as follows:
, 
, 
, and 
. In the above equation, 
is the time interval (s) and 
is the projected plan area (m2) of the control volume, Ω.
, 
, 
, 
and 
are coefficients, 
, 
, and 
includes all the remaining terms, defined as the source term. Further manipulation of Equation (6) results in
represents the coefficients of neighboring nodes.
represents the coefficients of neighbouring nodes, 
and 
. For regular grids, 
can be used.The convection term in Equations (2) and (3) can be discretized using several schemes. In this study, the forward Euler scheme is used to discretize the unsteady term, and the upwind scheme is applied to the convective term due to its strong unwinding capability. It needs to be noted that the coefficient matrices of the discretized equations (Equations (7) and (8)) lead to the same due to the non-staggered grid arrangement. Here, the non-staggered grid approach reduces the computer memory and computational cost requirement compared to the staggered grid arrangement.
Coupling of velocity and water level
and 
are the u- and v-velocities of the previous iteration step, 
and 
represent the first term on the right-hand side of Equations (7) and (8), respectively, and 
and 
are the relaxation factors. The value of the under-relaxation coefficient is arbitrarily chosen, but Wu et al. (2011) and Wu & Lin (2015) have shown that a value of 0.8 is suitable for flood flow studies.
and 
as follows:
and 
are the approximate values of u- and v-velocities, respectively and 
and 
denote the approximate values of 
and 
, respectively.
is the velocity correction field for the x-direction.
is the velocity correction field for the y-direction and the velocity correction is given by
and 
and 
. 
and 
are the momentum weighing factors.From the above derivation, it is clear that the Rhie and Chow technique is composed of momentum (first term on the right side of Equations (17) and (18)) and linear (second term on the right side of Equations (17) and (18)) interpolation terms. It is evident that when the under-relaxation factors are smaller enough, the main contribution will be from the second term, and spurious oscillations may re-appear. On the other hand, when the under-relaxation factors are chosen to be high enough, the main contribution will be from the first term, and the final solution might not converge. In this way, the under-relaxation factor plays an important role in the Rhie and Chow technique to reach a final converged solution in the domain (Majumdar 1988; Martínez et al. 2017).
and 
are the momentum weighing factors. These factors are used to interpolate the values of the variables at an interface from the corresponding neighboring cell values. The coefficients, 
and 
, stand for 
and 
when applying in Equations (7) and (8) on the cells centered at W and S (Figure 1), respectively.
, 
, 
and 
are the approximate fluxes (m3/s) evaluated at the interfaces w and s using the equations 
and 
, in which 
and 
are the approximate velocity components evaluated at the cell faces.
Adopted methodology.
Improvements over the Wu-2D model
The present study proposes two important modifications in the 2D model developed by Wu et al. (2011) (hereafter, the model is referred to as the Wu-2D model) in order to enhance numerical stability and computational efficiency. The improvements are demonstrated on regular grids, and they are discussed below:
- So far, the momentum weighting factor is defined using the cell geometry (i.e., as a geometrical parameter), and this turns out to be a fixed value of 0.5 for regular grids. As a result, the constant momentum weighting factor for the entire computational domain and throughout the simulation time does not consider the effect of local flow dynamics. The scaling effect due to the momentum correction factor should consider the local flow dynamics so that the momentum fluxes are computed accurately. This study, therefore, proposes that the momentum weighting factor be calculated using local flow dynamics. Roe (1981) proposed an averaging scheme for estimating relatively smooth interface values of flow variables for simulating flow with discontinuity. The same scheme was successfully adopted by Kuiry et al. (2010). The Roe averaging is extended herein to derive the momentum weighting factor to be used in the Rhie and Chow interpolation method. The factors are now expressed in terms of local flow depths and vary spatially as well as temporally:
- So far, the water surface at an interface is calculated using the various interpolation methods using the geometrical parameters of the neighboring cells. Ying et al. (2004) proposed a Courant number-based weighting factor for shock-capturing schemes to avoid numerical instability arising from the interface values. A similar expression is derived in this study in order to compute the water surface elevation at an interface; for example, the w-interface of the P-cell can be written as:where(24)
, and 
. The weighted 
values at the interfaces of a cell are, therefore, used in the momentum equations (Equations (11) and (12)) to compute the velocity components. It can be observed from Equation (24) that the proposed weighting factor depends on local velocity, flow depth, and computational time step. In this way, the weighting factor varies spatially and temporally as the solution advances with time. Therefore, the proposed water surface weighting factor is dynamic in nature. Sridharan et al. (2020) demonstrated that using such a dynamic weighting factor for estimating the numerical flux at the interface of a local-inertial model could significantly improve the model stability.
Solution of the discretized equations
Each control volume generates three equations for three unknowns (h, u, and v for mass and momentum equations), and hence, the number of algebraic equations and unknowns is equal to the number of control volumes. Eventually, a penta-diagonal matrix is formulated that has non-zero coefficients in five diagonals. The matrix is then solved iteratively using the strongly implicit procedure (SIP) proposed by Stone (1968). The SIP treats all unknowns at five points in a strongly implicit manner, and thus, the coupling between the unknowns at each grid point is greatly enhanced (Azevedo et al. 1988; Daoud et al. 2008). The discretized equations (Equations (11), (12), (17), (18), (23), (14), (16), (19) and (20)) are solved sequentially using the SIMPLEC algorithm as described by Wu et al. (2011).
Wetting and drying algorithm
The waterfront can move over a dry bed, and a dry bed can reappear during a receding flood. As a result, the cells at the moving waterfront can become wet or dry during a flood event. This process is called the wetting and drying phenomenon. The threshold depth approach, as described in Wu et al. (2011) and Wu & Lin (2015), is adopted in this study to capture the moving waterfront on a dry bed. The threshold flow depth of 
m is used to judge whether a cell is dry or wet. If the flow depth at the centroid of a cell is lower than the prescribed threshold value, then that cell is considered a dry cell; otherwise, it is a wet cell. A dry cell can only receive water and cannot donate water to its neighbors. This simple algorithm effectively captures the wetting and drying phenomenon.
MODEL PERFORMANCE
The performance of the proposed model is evaluated through simulations of three cases: 1D steady-state flow over a hump without friction, a hypothetical flood in the Severn River, UK, and flood propagation in an experimental river-network floodplain system. The 1D steady-state flow over a hump without friction case assesses the model’s accuracy against the analytical solution and demonstrates its rate of convergence. The hypothetical flood in the Severn River showcases the effectiveness of the interpolation scheme in generating non-oscillatory solutions compared to existing models. The flood propagation in an experimental river-network-floodplain system evaluates the model’s accuracy and computational efficiency by comparing the simulated water depth and inundation extent. Additionally, to validate the proposed advancements, the authors implement the Wu-2D model following Wu et al. (2011) and Wu & Lin (2015), maintaining similar code structures for comparative analysis. In the absence of observed data, results from existing models such as TELEMAC-2D (http://www.opentelemac.org) and HEC-RAS-2D (https://www.hec.usace.army.mil) are used for comparison purposes, alongside analytical solutions and observed data where available.
The TELEMAC-2D model is selected due to its capability to reproduce complex flow dynamics (Horritt & Bates 2002; Costi et al. 2018; Kaveh et al. 2019; Li et al. 2019). The TELEMAC-2D model uses the FEM and FVM solvers to solve the 2D depth-averaged SWEs on an unstructured triangular mesh. TELEMAC-2D is developed by the National Hydraulics and Environment Laboratory of the Research and Development Directorate of the French Electricity Board (EDF-R&D). The full description of the model can be found in Hervouet (2007). The HECRAS-2D model (version 5.0.5) also solves the 2D depth-averaged SWEs (Equations (1), (2) and (3)) using a hybrid discretization scheme combining FDM and FVM (Brunner 2016). FDM is used to discretize the governing equations when the grid is locally orthogonal, whereas FVM is used to discretize the governing equations when the grid is not locally orthogonal. The HECRAS-2D is developed by the US Army Corps of Engineers, Hydrologic Engineering Center with US Federal Government resources (https://www.hec.usace.army.mil/). The HECRAS-2D model can even handle the sudden rush of water into an area (Brunner 2016). Many authors have applied the model (Farooq et al. 2019; Rangari et al. 2019) for simulating various flooding scenarios. The accuracy of HECRAS-2D is reported to be within the range of reasonable to excellent. The models are evaluated in terms of computation time using a desktop computer with the following configuration: IntelⓇCore
i5-4570 CPU @ 3.20GHz, 8GB RAM, 4 cores, 4 logical processors, and a 64-bit operating system.
Case I: steady state 1D flow over a hump without friction
) of the hump is provided below.
) is calculated by
and 
are the local water depths (m) at the current (n) and previous time levels (
), respectively. The 
value is used for convergence test for a steady state solution and given as 
.For subcritical flow, the discharge per unit width of 
= 4.22 m2/s is specified at upstream of the channel, and a water depth of h = 2 m is specified at its downstream. The channel grid is generated with an equal spacing of 0.25 m in both directions.
Comparison of steady sub-critical flow over a hump without friction: (a) water surface elevation, (b) relative convergence error, and (c) total energy gradient.
Comparison of steady sub-critical flow over a hump without friction: (a) water surface elevation, (b) relative convergence error, and (c) total energy gradient.
Case II: hypothetical flood in the Severn River, UK
Flood simulation in the Severn River: (a) topography along with virtual gauges, and (b) hypothetical flood hydrograph.
Flood simulation in the Severn River: (a) topography along with virtual gauges, and (b) hypothetical flood hydrograph.
A uniform Manning’s n value of 0.035 m
s is defined throughout the flow domain (Aricò et al. 2016). For TELEMAC-2D, we have chosen the edge-by-edge implementation of upwind explicit finite volume and mass conservative schemes for computing velocity and water depth, respectively, with no turbulence model (Hervouet 2007). Since this case is completely hypothetical on a real bathymetry, TELEMAC-2D and HECRAS-2D model results are used as reference solutions on two different types of grids. The Wu-2D model is used to compare improvements directly due to the modifications proposed in the IITMFLO-2D model.
(Wu et al. 2011), and the results are presented in Figure 5. In addition, the results for several under-relaxation factors (
) are shown in Figure 6 to investigate its effect on the accuracy. Overall, all the improvements in the model are seen to follow a similar trend as predicted by TELEMAC-2D and HECRAS-2D. The prediction of flood peaks at different gauges is in good agreement with each other. Still, at some locations, the trends exhibit slight deviations along the rising and falling limbs of the hydrographs. The small differences might have arisen from the numerical approaches used in the models. However, clear improvements in results can be observed in the case of IITMFLO-2D due to the proposed modifications, as discussed below. It can be seen that the Wu-2D model produces some oscillations in the flow depth patterns, especially at the gauge P6 (Figure 5f). Meanwhile, the same Wu-2D model with water surface weighting factors (w) as proposed in this study (i.e., Rhie and Chow method + w) resulted in a smooth flow depth pattern at this location. Therefore, this result confirms that estimating the interfacial water surface level using the proposed weighting factors (Equation 24) offers greater model stability than the original Wu-2D model. The possible reason for such improvement may be attributed to the modified weighting factors, which consider the local flow dynamics, and the values are automatically adjusted accordingly to compute momentum fluxes.
Second improvement over the Wu-2D model is the introduction of Roe-based momentum weighting factors (f). The effect of this improvement is now realized in combination with the first improvement (i.e., Rhie and Chow method + w+ f). This improvement can simulate flow depths similar to the first improvement at all the stations. Surprisingly, significant improvement is observed in the overall computation time, as summarized in Table 1. This improvement allows us to take a larger time step size, and hence, total computation time is reduced in comparison to the Wu-2D model.
Implementation of the Rhie and Chow interpolation technique in the Wu-2D model limits the selection of the
value for a stable solution. In this case, a stable range is found to be between 0.4 and 0.9 depending on the value of Reynolds number (Majumdar 1988). However, it is suggested that the optimum 
value should be selected by trial and error (Wu et al. 2011) for better accuracy. However, such a process needs observed data and can be cumbersome for solving real-world problems. The proposed improvements (i.e., Rhie and Chow method + w + f) in the Wu-2D model are found to provide oscillation-free and accurate results (Figure 6) for any value in a wide range of 
(
). The proposed improvements also significantly reduce the overall computation time, as observed in Table 1. Therefore, the Rhie and Chow interpolation method with other two improvements (i.e., Rhie and Chow method + w + f) can be used to reduce the dependency on the 
value in a wide range with a significant computational advantage over the Wu-2D model. It is also clear from Table 1 that IITMFLO-2D is computationally more efficient than HECRAS-2D and TELEMAC-2D models. Thus, the proposed modifications can avoid the trial and error method for choosing the 
value for stability and accuracy. Henceforth, the IITMFLO-2D model refers to the combination of (Rhie and Chow method + w + f) in the Wu-2D model in all the applications and analyses.Accuracy assessment is performed using the root mean square error (RMSE) by comparing the simulated results derived from IITMFLO-2D and Wu-2D models with that of TELEMAC-2D and HECRAS-2D model results (Figure 7). It was found that the RMSE was lower for both models when compared with the TELEMAC-2D results. Hence, the time series of water depth closely matches the TELEMAC-2D results compared to the HECRAS-2D results. Figures 5 and 6 show that adopting the Rhie and Chow method + w+ f in IITMFLO-2D provides a similar level of accuracy to that of the Rhie and Chow method in the Wu-2D model. This is due to the fact that both models maintain a similar level of accuracy at all the gauges except at one gauge. However, IITMFLO-2D is observed to maintain accuracy for a wide range of
values, which can avoid trial and error, unlike the Wu-2D model. In addition, the proposed modifications in IITMFLO-2D provide a boost to run faster than the Wu-2D model.
Comparison of time-series water depth at the virtual gauges over the floodplain and in the river: (a) P1, (b) P2, (c) P3, (d) P4, (e) P5 and (f) P6.
Comparison of time-series water depth at the virtual gauges over the floodplain and in the river: (a) P1, (b) P2, (c) P3, (d) P4, (e) P5 and (f) P6.
Comparison of time-series water depth at the virtual gauges for different 
values. (a) P1, (b) P2, (c) P3, (d) P4, (e) P5 and (f) P6.
Comparison of time-series water depth at the virtual gauges for different values. (a) P1, (b) P2, (c) P3, (d) P4, (e) P5 and (f) P6.
Comparison of RMSE of simulated water depths derived from IITMFLO-2D and Wu-2D models with respect to TELEMAC-2D and HECRAS-2D models: (a) P1 and (b) P2.
Comparison of RMSE of simulated water depths derived from IITMFLO-2D and Wu-2D models with respect to TELEMAC-2D and HECRAS-2D models: (a) P1 and (b) P2.
Relative computation time taken for the Severn River simulation
| Numerical schemes | Rhie and Chow | Rhie and Chow + w | Rhie and Chow +  | TELEMAC-2D | HECRAS-2D |
|---|---|---|---|---|---|
| Computation time | 1.35 | 1.37 | 1.00 | 2.10 | 1.21 |
| Numerical schemes | Rhie and Chow | Rhie and Chow + w | Rhie and Chow + | TELEMAC-2D | HECRAS-2D |
|---|---|---|---|---|---|
| Computation time | 1.35 | 1.37 | 1.00 | 2.10 | 1.21 |
Case III: flood propagation in an experimental river-network-floodplain system
Measurement of water depth for steady peak floods.
The set-up bathymetry is defined in the form of a raster grid of resolution 3 cm DEM as described in Mali & Kuiry (2018). The IITMFLO-2D and Wu-2D models directly take advantage of the raster grid to generate a computational mesh of 110,160 cells. On the other hand, for TELEMAC-2D, the set-up bathymetry is described by 111,655 triangles with local refinement along the channels and a smooth transition is maintained from finer grids to the coarser grids on the floodplains. Therefore, it is evident that the topographical features of the set-up are better described in TELEMAC-2D compared to the DEM used by IITMFLO-2D and Wu-2D models. The edge-by-edge implementation of upwind explicit finite volume and mass conservative schemes are used in TELEMAC-2D to compute velocity and water depth. It should be noted that the HECRAS-2D model could not be used with the least allowable time step of less than 0.1 s, which does not provide converged solutions for this case due to the requirement of fine grids to describe the bathymetry of the set-up. A single value of Manning’s roughness coefficient 
is considered for both channels and floodplains as the set-up is made up of the same material with a similar kind of surface finishing. Therefore, a common value of 
m
s for the bed surface is used for simulations (Mali & Kuiry 2020; Mali et al. 2020). IITMFLO-2D and Wu-2D models are run for the same value of 
for all the cases because the Wu-2D model provides the best accuracy for this value. However, it must be reiterated that the IITMFLO-2D model maintains a similar level of accuracy for any value of 
within a wide range (
).
Steady flood flow in the river-network-floodplain system
The flood inundation due to peak flow is simulated for two experimental observations. Inflow discharge is prescribed at the inlet, and observed water surface elevations at three outlets are imposed as boundary conditions (Table 2). The initial condition in the channel network is specified as a constant water surface elevation of 6.13 m from a datum, i.e., on average, a constant flow depth of 0.13 m. The inflow rate is increased progressively to the actual inflow discharge to avoid any model instability. The simulations are performed until the steady-state condition throughout the set-up is attained.
Boundary conditions for steady peak flood in the set-up
| Boundary condition | ||||
|---|---|---|---|---|
| Downstream (m) | ||||
| Discharge | Upstream inflow (m3/s) | Left ( ) | Middle ( ) | Right ( ) |
 | 0.088 | 0.216 | 0.205 | 0.184 |
 | 0.098 | 0.223 | 0.219 | 0.192 |
| Boundary condition | ||||
|---|---|---|---|---|
| Downstream (m) | ||||
| Discharge | Upstream inflow (m3/s) | Left () | Middle () | Right () |
| | 0.088 | 0.216 | 0.205 | 0.184 |
| | 0.098 | 0.223 | 0.219 | 0.192 |
The location of boundary conditions and various point gauges to measure flow depth are shown in Figure 8. The flow depths are measured at 13 channel cross-sections and at 35 locations on the floodplains and are marked with filled circles (
) (Figure 8).
Comparison of simulated and observed water depths for two peak flood inflows. (a) 
Channel water depth, (b) 
Channel water depth, (c) 
Floodplain water depth and (d) 
Floodplain water depth.
Comparison of simulated and observed water depths for two peak flood inflows. (a) Channel water depth, (b) Channel water depth, (c) Floodplain water depth and (d) Floodplain water depth.
However, the RMSE values of IITMFLO-2D model results are in the range of RMSEs of other models. TELEMAC-2D has been found to provide better RMSE than others. This may be attributed to the better representation of channels and floodplain topography by the locally refined unstructured mesh used in TELEMAC-2D. Nevertheless, the IITMFLO-2D and Wu-2D models have given satisfactory results. It can be inferred that the high-inertia flow in the channel network is not modeled as well as on the floodplain due to the presence of convergence, divergence, and curvature. These places are commonly known to have 3D-complex flow patterns. In addition, the 2D SWEs are based on the assumption of hydrostatic pressure, which can not be held true at these locations. For example, at a curvature or bend, the geometry-driven secondary currents predominantly affect the longitudinal velocity, leading to the formation of complex velocity patterns and, hence, the water surface elevation. The basic flow structure at a river bend is a helix. The higher momentum near-surface flow attempts to move toward the outer bank because of centrifugal acceleration, and conversely, the slow-moving near-bed flow moves toward the inner bank of the river (Ferguson et al. 2003). Thus, the 3D flow structure at these locations becomes challenging for the 2D models to capture appropriately. Despite that, the predictions of all the 2D models are noticeably appreciable.
inflow discharge, while a close prediction is observed for 
inflow discharge. The accuracy of inundation extent is further estimated by using the measure of fit (f) (Horritt & Bates 2001; Mali & Kuiry 2020), as given below.
and 
represent the sets of cells or areas or pixels classified as wet predicted by the model and in the observed inundation map, respectively. The 
(.) gives a number of members of the set. The f shows the performance of the model as the range can lie between 
for a perfect prediction and 
for no region correctly predicted with respect to the observed data. Comparison of inundation area and measure of fit (f) during steady peak floods in the set-up
| Discharge | Inundation area (m2) | Measure of fit (f) | |||||
|---|---|---|---|---|---|---|---|
| IITMFLO-2D | Wu-2D | TELEMAC-2D | Observed | IITMFLO-2D | Wu-2D | TELEMAC-2D | |
| Q1 | 50.1 | 49.1 | 49.1 | 46.52 | 0.89 | 0.89 | 0.91 |
| Q2 | 53.45 | 52.3 | 52.9 | 53.27 | 0.92 | 0.91 | 0.92 |
| Discharge | Inundation area (m2) | Measure of fit (f) | |||||
|---|---|---|---|---|---|---|---|
| IITMFLO-2D | Wu-2D | TELEMAC-2D | Observed | IITMFLO-2D | Wu-2D | TELEMAC-2D | |
| Q1 | 50.1 | 49.1 | 49.1 | 46.52 | 0.89 | 0.89 | 0.91 |
| Q2 | 53.45 | 52.3 | 52.9 | 53.27 | 0.92 | 0.91 | 0.92 |
Comparison of water depth variation within the set-up for 
peak flood: (a) IITMFLO-2D, (b) Wu-2D, (c) TELEMAC-2D and (d) observed.
Comparison of water depth variation within the set-up for peak flood: (a) IITMFLO-2D, (b) Wu-2D, (c) TELEMAC-2D and (d) observed.
Comparison of water depth variation within the set-up for 
peak flood: (a) IITMFLO-2D, (b) Wu-2D, (c) TELEMAC-2D and (d) observed.
Comparison of water depth variation within the set-up for peak flood: (a) IITMFLO-2D, (b) Wu-2D, (c) TELEMAC-2D and (d) observed.
From Table 3, it can be seen that both IITMFLO-2D and Wu-2D exhibit a similar level of accuracy (f) as compared to the TELEMAC-2D model. However, the Wu-2D model needs to be iterated for 
value for better accuracy.
Unsteady flood flow in the river-network-floodplain system
Boundary conditions: (a) unsteady state stepped-discharge-hydrograph at the upstream and flow depth evolution at the downstream of the (b) left channel, (c) middle channel, and (d) right channel.
Boundary conditions: (a) unsteady state stepped-discharge-hydrograph at the upstream and flow depth evolution at the downstream of the (b) left channel, (c) middle channel, and (d) right channel.
) on the ortho-image (Figure 13), are used to validate the simulated results. The observed flood depths are available at every 2 s interval. 
Measurement of time-series water depth for unsteady flow hydrographs.
to investigate its sensitivity. The following observations are made during the unsteady flow: All the figures provide a detailed account of the propagation of flood waves, the rise and fall of water levels, and the time to reach peak water levels at different gauges. In general, it can be seen from the figures that most of the features of the water levels are wellmodeled, notably producing similar patterns as seen in the observed data.
In the channels, IITMFLO-2D and Wu-2D models appear to be slightly overpredicting at most of the gauges except P12 and P16 (Figures 14 and 16). Except for P9, TELEMAC-2D predicts water levels reasonably better compared to all the other models. It is interesting to note that overprediction by all the models is found at the gauges where water overflows onto floodplains and in the vicinity of junctions. Such gauges are identified as P7 and P9. The more pronounced fluctuations in the channel can be seen at P9 (Figures 14c and 16c), which is located just after the divergence D1. Also, at P9, the observed water level initially rises gradually, and then the rate of rising decreases for a certain period of time, followed by an abrupt rise to the peak. A very similar trend can also be seen on the recession limb. This may be due to the fact that the water depth on the downstream side of the divergence at this location varies with the flow rate after striking the wall at the divergence. Thus, the wave with a certain height moves back and forth, and hence, the sudden increase and decrease in the water level are observed at this location. However, such features could not be reproduced by the 2D models.
The observed and simulated maximum flow depths in the channels and over the floodplains depend strongly on the location of the gauges. This is due to the complex flow dynamics and wave patterns caused by the uneven topography. For example, 21 and 22-cm flow depths are observed at P7, where the water level rises after striking the divergence. At P7 (Figures 14b and 16b), TELEMAC-2D and Wu-2D models over-predict the water level, whereas the peak water level is well predicted for low inflows by TELEMAC-2D. For higher inflow, i.e., for the second case, the peak is underpredicted. The variations may be attributed to the complex flow dynamics at a channel-floodplain interface and junction. The flow at these locations is 3D in nature on account of the lateral and vertical exchange of momentum between the flow layers and the presence of dominant secondary circulation. These features can be linked to the vertical velocity component, which is completely absent in the 2D models. On the other hand, it can be observed that the same 2D models are able to simulate similar peaks and patterns of water levels at other locations, where the flow dynamics are relatively simple (Figures 14d, 14e, 16d and 16e).
Overall, the patterns of water levels at different gauges over the floodplain are reasonably well predicted by all the models. However, some gauges (e.g., P1, P11, P13, P14, and P15) exhibit slight discrepancies in modeling the water level variations. The movement of waterfronts is in close conformation with the observed values at some gauges (e.g., P5, P8, P10) as depicted in Figures 15 and 17. However, the modeled waterfronts enter earlier than the actual arrival time onto the floodplains and islands at some locations, for example, at P6, P11, P13, and P14. On the other hand, the arrival of the waterfronts predicted by TELEMAC-2D is delayed at the gauges P1, P4, and P6 located over the left floodplain, as shown in Figure 17. This suggests that the overflow onto the left floodplain is delayed. Conversely, the other models have predicted the early arrival of the waterfront to the left floodplain. The actual water levels at the time of the peak of the hydrographs are captured by the models at most of the locations (e.g., P5, P8, P10, P13). Among all the models, it is seen that the results of IITMFLO-2D are in good agreement with the observed data.
Different
values are tried for IITMFLO-2D and Wu-2D models for simulating the flood dynamics in the set-up (Figures 16 and 17). IITMFLO-2D and Wu-2D models are able to simulate the flood events for various values of 
and follow a similar trend as that of the observed data. It is worth mentioning that IITMFLO-2D simulates the events for a wide range of 
(0.2–0.9) values compared to the Wu-2D model (0.4–0.9). It is thus found that IITMFLO-2D works for a wide range of 
values, but a minute variation in the simulation results cannot be avoided. Moreover, IITMFLO-2D is able to cope with oscillations arising from the linear interpolation term while using a small enough 
value. The model is also able to converge while using a high enough 
value.
Comparison of modeled and observed flow depth evolution in the channels at the gauges for Q1 hydrograph: (a) P3, (b) P7, (c) P9, (d) P12, and (e) P16.
Comparison of modeled and observed flow depth evolution in the channels at the gauges for Q1 hydrograph: (a) P3, (b) P7, (c) P9, (d) P12, and (e) P16.
Comparison of modeled and observed flow depth evolution over the floodplain at the gauges for Q1 hydrograph: (a) P1, (b) P2, (c) P4, (d) P5, (e) P6, (f) P8, (g) P10, (h) P11, (i) P13, (j) P14, and (k) P15.
Comparison of modeled and observed flow depth evolution over the floodplain at the gauges for Q1 hydrograph: (a) P1, (b) P2, (c) P4, (d) P5, (e) P6, (f) P8, (g) P10, (h) P11, (i) P13, (j) P14, and (k) P15.
Comparison of modeled and observed flow depth evolution in the channels at the gauges for Q2 hydrograph. (a) P3, (b) P7, (c) P9, (d) P12 and (e) P16.
Comparison of modeled and observed flow depth evolution in the channels at the gauges for Q2 hydrograph. (a) P3, (b) P7, (c) P9, (d) P12 and (e) P16.
Comparison of modeled and observed flow depth evolution over the floodplain at the gauges for Q2 hydrograph. (a) P1, (b) P2, (c) P4, (d) P5, (e) P6, (f) P8, (g) P10, (h) P11, (i) P13, (j) P14 and (k) P15.
Comparison of modeled and observed flow depth evolution over the floodplain at the gauges for Q2 hydrograph. (a) P1, (b) P2, (c) P4, (d) P5, (e) P6, (f) P8, (g) P10, (h) P11, (i) P13, (j) P14 and (k) P15.
values can also be seen through RMSE values (Figure 19). The small 
values show a bit higher RMSEs for both the models and vice-versa. Thus, the 
values appear to potentially influence the simulated flow dynamics to a certain extent. 
RMSE calculated between observed and simulated water depth for Q1 hydrograph: (a) Channel and (b) floodplain.
RMSE calculated between observed and simulated water depth for Q1 hydrograph: (a) Channel and (b) floodplain.
RMSE calculated between observed and simulated water depth for Q2 hydrograph: (a) Channel and (b) floodplain.
RMSE calculated between observed and simulated water depth for Q2 hydrograph: (a) Channel and (b) floodplain.
Flood inundation extent extracted using visual interpretation technique for Q
hydrograph at the time instants: (a) 
s, (b) 
s and (c) 
s.
Flood inundation extent extracted using visual interpretation technique for Q hydrograph at the time instants: (a) s, (b) s and (c) s.
Flood inundation extent extracted using visual interpretation technique for Q
hydrograph at the time instants: (a) 
s, (b) 
s and (c) 
s.
Flood inundation extent extracted using visual interpretation technique for Q hydrograph at the time instants: (a) s, (b) s and (c) s.
The area of dynamic inundation extents and the measure of fit during the passage of flood hydrographs are estimated and presented in Table 4. One can easily observe that the inundation area predicted by IITMFLO-2D and Wu-2D models is rather closer to the observed area except at the initial stage when some of the channels have just started overflowing. The largest difference in the inundation area by the models is observed at 
s due to the overspreading of water on the floodplains. At all other times, the accuracy of the simulated flood extents is found to be reasonably accurate with 
for all the cases.
Comparison of area inundated and measure of fit derived from the 2D models
| Hydro-graph | Time instant (mins) | Inundation area (m2) | Measure of fit (f) | |||||
|---|---|---|---|---|---|---|---|---|
| IITMFLO-2D | Wu-2D | TELEMAC-2D | Observed | IITMFLO-2D | Wu-2D | TELEMAC-2D | ||
| Q1 | 20 | 22.5 | 24.2 | 20.07 | 16.35 | 0.71 | 0.64 | 0.77 |
| 40 | 49.22 | 50.1 | 48.12 | 47.8 | 0.9 | 0.9 | 0.91 | |
| 60 | 39.75 | 40.5 | 35.07 | 38.61 | 0.78 | 0.76 | 0.87 | |
| Q2 | 20 | 23.1 | 25.8 | 21.9 | 15.9 | 0.7 | 0.65 | 0.78 |
| 40 | 50.62 | 51.4 | 46.71 | 50.25 | 0.9 | 0.9 | 0.85 | |
| 60 | 47.2 | 48.3 | 45.5 | 46.22 | 0.89 | 0.89 | 0.89 | |
| Hydro-graph | Time instant (mins) | Inundation area (m2) | Measure of fit (f) | |||||
|---|---|---|---|---|---|---|---|---|
| IITMFLO-2D | Wu-2D | TELEMAC-2D | Observed | IITMFLO-2D | Wu-2D | TELEMAC-2D | ||
| Q1 | 20 | 22.5 | 24.2 | 20.07 | 16.35 | 0.71 | 0.64 | 0.77 |
| 40 | 49.22 | 50.1 | 48.12 | 47.8 | 0.9 | 0.9 | 0.91 | |
| 60 | 39.75 | 40.5 | 35.07 | 38.61 | 0.78 | 0.76 | 0.87 | |
| Q2 | 20 | 23.1 | 25.8 | 21.9 | 15.9 | 0.7 | 0.65 | 0.78 |
| 40 | 50.62 | 51.4 | 46.71 | 50.25 | 0.9 | 0.9 | 0.85 | |
| 60 | 47.2 | 48.3 | 45.5 | 46.22 | 0.89 | 0.89 | 0.89 | |
Relative computation time for the models during unsteady flow runs
| Models | IITMFLO-2D | Wu-2D | TELEMAC-2D | |||||
|---|---|---|---|---|---|---|---|---|
 |  |  |  |  |  |  | ||
| Computation time | 1.07 | 1.04 | 1.00 | 1.00 | 1.51 | 1.48 | 1.47 | 2.51 |
| Models | IITMFLO-2D | Wu-2D | TELEMAC-2D | |||||
|---|---|---|---|---|---|---|---|---|
| | | | | | | ||
| Computation time | 1.07 | 1.04 | 1.00 | 1.00 | 1.51 | 1.48 | 1.47 | 2.51 |
Thus, from RMSE of time-series of water levels at different gauges and measure of fit, f, values for dynamic inundation extents, it can be stated that all the models are performing reasonably well and providing almost a similar level of accuracy. However, TELEMAC-2D is found to be slightly more accurate than the others. It may be attributed to the mesh refinement performed for TELEMAC-2D simulations. In contrast, the DEM of 3 cm resolution mesh in IITMFLO-2D and Wu-2D is found inadequate to represent some topographical features accurately. Nevertheless, the accuracy of the IITMFLO-2D and Wu-2D models is quite satisfactory. From the above comparisons, it can be concluded that the proposed IITMFLO-2D model is able to reconstruct the unsteady state floods observed in the river-network-floodplain set-up.
It should be stated that the proposed IITMFLO-2D does not notably improve accuracy against the Wu-2D model, even after some modifications are implemented. However, the developed model outperforms both Wu-2D and TELEMAC-2D models in terms of computation time (Table 5). For unsteady-state runs, IITMFLO-2D is observed to be faster than Wu-2D and TELEMAC-2D models by 
1.47 and 
2.51 times, respectively, which turn out to be 
47% and 
151%, respectively, improvement in overall computation time. Such an improvement in overall computation time will definitely play an important role while using the models for real-time flood forecasting.
In a nutshell, the accuracy of both IITMFLO-2D and Wu-2D models is found to be in good agreement with the observed data. However, some discrepancies are seen in the simulated and observed data in the set-up. The discrepancies can be attributed to the wide range of uncertainties associated with the numerical and experimental approaches. The numerical solution is the approximate approach to the solution of the non-linear flow dynamic equations. The numerical schemes rely on various issues, such as the grid size, time step, and calibration procedure. Thus, the discrepancies between the observed and simulated data can always be expected. Firstly, the river-network-floodplain topography of the set-up is highly complex for the given geometrical scale and the flow dynamics due to the presence of junctions, islands, bends, or curves in the channels, highly uneven floodplains, flow interaction at the interface of channel and floodplain, etc. Also, the 2D SWEs cannot accurately represent the 3D flow field, as discussed before. Secondly, although data collection was performed using highly sophisticated instruments and sensors, still some levels of manufacturing uncertainties could not be avoided in the experimental data (Mali & Kuiry 2020; Mali et al. 2020). Therefore, the limitations in the numerical approach and certain degree of errors in experimental data are also the sources of discrepancies between the observed and simulated data. Despite the limitations, the proposed IITMFLO-2D and Wu-2D models are able to reproduce the experimental observations on floods with sufficient accuracy. The proposed improvements in IITMFLO-2D have prompted the model to run faster than the Wu-2D model.
MODEL LIMITATIONS AND FUTURE RECCOMMENDATIONS
The presented 2D model can be enhanced by implementing a block-structured or unstructured quadrilateral mesh, which allows for a more precise representation of hydraulically significant topographical features. Additionally, the model’s performance could be further optimized by parallelizing the FORTRAN code, utilizing multi-threading and/or multiprocessing through distributed memory parallelism techniques such as MPI (Message Passing Interface). Moreover, recoding the model for a GPU architecture could significantly decrease simulation times, making it more practical for studying the long-term evolution of flood events. It should be noted that the developed model is currently tailored for slow-rising floods where the flow remains primarily subcritical. In contrast, floods caused by events like dam breaks or levee breaches can involve transcritical flow regimes, which would necessitate the use of a shock-capturing numerical scheme.
CONCLUSIONS
IITMFLO-2D, a 2D flood simulation model, is developed based on the Wu-2D model, using a Cartesian grid that allows direct use of DEM as a computational mesh. The model solves the depth-averaged 2D SWEs using the FVM and SIMPLEC algorithm, with all primitive variables stored at cell centroids to prevent the decoupling of velocity and water level fields. Two key improvements over the Wu-2D model include the introduction of a CFL-based weighting factor for interfacial water level computation and an explicit expression for the momentum weighting factor using local flow dynamics. The model’s accuracy is validated through analytical test cases and simulations of a hypothetical flood in the Severn River, UK, and experimental flood events in a river-network-floodplain set-up, comparing its performance with TELEMAC-2D and HECRAS-2D models. The following specific conclusions are drawn from the study:
It is found that the simulated results of all the models are in good agreement with the observed data. The accuracy of the models is evaluated through RMSE between the simulated and observed flood depths and a measure of fit between simulated and observed flood extents.
IITMFLO-2D and Wu-2D exhibit comparable accuracy levels, while TELEMAC-2D provides the highest accuracy due to its local mesh refinement capabilities.
IITMFLO-2D offers improved stability as compared to the Wu-2D model. This is due to incorporating the CFL-based weighting factor instead of the simple linear interpolation used in the Wu-2D model.
IITMFLO-2D surpasses the performance of the other models in simulating experimental flood events, with significant reductions in computation time. It is approximately 1.47 times faster than Wu-2D and 2.51 times faster than TELEMAC-2D. The increased computation speed is due to the model’s dynamic momentum weighting factor, allowing for larger time steps.
The simulations indicate that IITMFLO-2D can effectively handle a wide range of under-relaxation factors (
) with minimal impact on accuracy, unlike the Wu-2D model. This flexibility adds to the robustness of IITMFLO-2D.The enhancements in IITMFLO-2D significantly boost computation speed, making it suitable for real-time flood forecasting, early warning dissemination, and flood mitigation planning and management. The model’s improvements can be applied to other numerical models, including those utilizing quadrilateral or unstructured computational meshes, to capture topographical features better.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.
CONFLICT OF INTEREST
The authors declare there is no conflict.
%20cropped.png?versionId=7142)