Comparing Combined 1D/2D and 2D Hydraulic Simulations Using High-Resolution Topographic Data: Examples from Sri Lanka—Lower Kelani River Basin

The application of numerical models to understand the behavioural pattern of a flood is widely found in the literature. However, the selection of an appropriate hydraulic model is highly essential to conduct reliable predictions. Predicting flood discharges and inundation extents are the two most important outcomes of flood simulations to stakeholders. Precise topographical data and channel geometries along a suitable hydraulic model are required to accurately predict floods. One-dimensional (1D) hydraulic models are now replaced by two-dimensional (2D) or combined 1D/2D models for higher performances. The Hydraulic Engineering Centre’s River Analysis System (HEC-RAS) has been widely used in all three forms for predicting flood characteristics. However, comparison studies among the 1D, 2D to 1D/2D models are limited in the literature to identify the better/best approach. Therefore, this research was carried out to identify the better approach using an example case study of the Kelani River basin in Sri Lanka. Two flood events (in 2016 and 2018) were separately simulated and tested for their accuracy using observed inundations and satellite-based inundations. It was found that the combined 1D/2D HEC-RAS hydraulic model outperforms other models for the prediction of flows and inundation for both flood events. Therefore, the combined model can be concluded as the better hydraulic model to predict flood characteristics of the Kelani River basin in Sri Lanka. With more flood studies, the conclusions can be more generalized.


Introduction
Among many other natural disasters, floods are very common throughout the world. They are dangerous and can impact a significant area. Annually, USD 662 billion of damage is estimated from these floods over the world [1]. The adverse impacts of floods range from direct impacts such as the loss of human life, damage to property, destruction of crops, loss of livestock to indirect impacts such as the spread of waterborne diseases, deterioration in water quality, etc. Therefore, the ability to predict the nature and the extent of a flood is very important to local decision makers as this would enable them to plan for such adverse impacts and minimize the damage.
It is essential to develop appropriate mitigation measures to minimize flood risk and flood damage. Structural and non-structural flood mitigation approaches are sometimes coupled together for better solutions; however, each approach has its own merits and demerits. The construction of dykes, levees, dams, flood control reservoirs, diversions, flood-ways, etc., are a few examples of structural measures. However, these structural measures might be costly, and implementation can be time-consuming. On the other hand, flood forecasting and warning, watershed management, floodplain development guidelines, insurance, and awareness programs fall into non-structural measures [2,3]. The structural measures have a certain design capacity, and, in cases of unexpected extreme flood events, failures are possible and, thus, the damage could be disastrous [4][5][6][7]. Therefore, nonstructural approaches are very important to promote settlement and reduce property damage, ensuring the safety and well-being of the public.
Recent studies revealed that there is an increasing trend of floods [8][9][10][11]. Rapid urbanization, unplanned and uncontrolled developments on floodplains, and the blockage of drainage paths associated with the increasing population can potentially increase the flood risk [12][13][14][15]. Impervious surfaces decrease the infiltration and thus increase the surface runoff. Therefore, the usual shape of the hydrograph is changed and flash floods are frequent. In addition, the changing climate has increased the number of floods and their magnitude in some areas [16][17][18][19]. Therefore, more adverse impacts can be witnessed and expected in the future. Thus, as it was stated, non-structural solutions are much needed for society.
The hydrodynamic simulations to forecast the inundation extent is considered as one of the widely used methods among non-structural measures. However, decision makers and researchers have diverse experiences regarding the selection of hydrodynamic models. A wide variety of numerical models developed by many researchers/decision makers to estimate/predict floods and flood-induced risk can be found in the literature. However, the selection based on the dimension of modelling (either one-dimensional (1D), two-dimensional (2D), or combined model (1D/2D)) are still problematic and questionable [20][21][22].
The Hydraulic Engineering Centre's River Analysis System (HEC-RAS) is a public domain software that is widely used in hydraulic-related applications throughout the world. Many researchers have conducted different studies on comparisons of HEC-RAS (i.e., 1D, 2D, or combined (1D/2D)) models with other open-source or perpetual licence hydrodynamic models such as TELEMAC-2D, LISFLOOD-FP, MIKE-11, MIKE-21, etc. [23,24]. Furthermore, it is reported that the prediction capability of models remains the same but performance varies with changes in the friction parameters. However, due to the limitations of 1D models in capturing the properties of the floodplain, coupled 1D/2D or 2D modelling approaches were suggested by some of the researchers [23,24].
Timbadiya et al. [25] suggested that the selection of numerical model type also plays a crucial role in the accuracy of hydraulic simulations similar to precise topography and channel geometry. In most cases, 1D models are now replaced by the 2D or combined 1D/2D models for higher performances. Hydraulic simulations for complex flood plain scenarios were usually carried out using the combined 1D/2D models [20,21]. The main channel is modelled as a 1D case, while the flood plain is modelled as a 2D case under the combined 1D/2D simulations. The literature showcases many studies based on the combined 1D/2D hydraulic modelling [20,[26][27][28]. However, limited applications can be found for comparison studies of 1D to 2D and combined 1D/2D. Vozinaki et al. [21] attempted to compare 1D to combined 1D/2D hydraulic simulations; however, it was not that comprehensive. In addition, there is no literature on the selection of 1D, 2D, and combined 1D/2D hydraulic models in the context of Sri Lanka, even though the country has major flood events annually.
Therefore, this study, for the first time, presents a comprehensive comparison study of 1D to 2D and combined 1D/2D hydraulic models using a case study in Sri Lanka. The lower stretch of the Kelani River with the influence of tides and the use of high-resolution data were integrated for the modelling purposes of the HEC-RAS hydraulic model. The calibration and validations were conducted for both flows and inundation extents with measured stages and satellite-observed inundation extents.

Hydraulic Modelling of River Flow
The schematic diagram of the river flow is expressed in Figure 1. The 2D flow conditions can be seen from the plan view ( Figure 1a) and the longitudinal view ( Figure 1c).

Governing Equations for 2D HEC-RAS Hydrodynamic Model
The HEC-RAS 2D model is developed to model and simulate complex floodplains in situations where the 1D flow is no longer valid enough for expected outcomes. The 2D unsteady flow varies with time and along two spatial dimensions. The governing laws of 2D unsteady flow are the conservation of mass (continuity) and conservation of momentum. The two-dimensional unsteady continuity equation is mathematically expressed in Equation (1).
∂H ∂t where H is the water surface elevation (m), h is the water depth (m), u and v are the depth-averaged velocities in the x and y direction (m/s). The conservation of momentum is calculated by Newton's second law of motion, which states that the sum of forces acting on an element equals the rate of change of momentum, which is written by considering the gravitational force, eddy viscosity, friction, and Coriolis effect [29]. The Coriolis effect and the eddy viscosity terms are usually neglected in the 2D model due to the size of the river basin and to maintain uniformity in the comparisons with the 1D flow and the unavailability of the required parameters for the calibration of the eddy viscosity coefficient [30,31]. Therefore, the modified full momentum equations can be written as Equations (2) Momentum balance in y direction where g is the gravitational acceleration (m/s 2 ), n is the Manning's coefficient, and R is the wetted perimeter (m). In order to reduce the computational time and the numerical instabilities, the HEC-RAS 2D unsteady flow Saint-Venant equations (shallow water equation) are often simplified using diffusive wave approximation. However, those simplifications are only valid for certain flow conditions. Nevertheless, the rivers that are influenced by tides are advised to incorporate the full momentum equations [31]. The Kelani River, which was used as an example in this study, has tidal influences [32,33] and, therefore, comprehensive analysis is recommended.

Governing Equations for 1D/2D HEC-RAS Hydrodynamic Model
The main river stretch is modelled as a 1D model, while the flood plain is a model with a 2D modelling approach in the coupled 1D/2D HEC-RAS model. The HEC-RAS uses a tight coupling technique. In other terms, both the 1D domain and the 2D domain are coupled on a time step basis. A lateral structure is used in order to establish the connection between the 1D and 2D domains [34].
The main river stretch solves the 1D Saint-Venant (shallow water) equation, which is comprised of continuity and momentum equations as shown in Equations (4) and (5).

∂A ∂t
where, A, Q, S 0 , and S f are the cross section area (m 2 ), the water flowrate (m 3 /s), the slope of the riverbed, and the energy slope, respectively. In addition, flow in the flood plain is calculated using 2D continuity equations and momentum equations as given in Equations (1)-(3).

Case Study Application
The Kelani River basin in Sri Lanka is the second-largest river basin based on the catchment area. River Kelani extends from the central hills of the country to Colombo, at a distance of about 145 km when measured along the channel centreline. The river basin is located between Northern latitude 6 • 47 to 7 • 05 and Eastern longitudes 79 • 52 to 80 • 13 , with a basin area of nearly 2230 km 2 as shown in Figure 2. Broadly, the river basin can be categorized into upper and lower basins. The upper basin features mountainous terrain, whereas the lower basin is generally flat terrain. The lower basin lies below the Hanwella River gauging station, which has an approximate area of 500 km 2 (refer to Figure 2, light brown area).
The upper catchment is predominantly covered with vegetation, whereas the lower catchment is heavily urbanized. The river basin receives an average annual rainfall of nearly 2400 mm and carries a peak discharge of 800-1500 m 3 /s during the monsoonal periods (i.e., especially in the southwest monsoon period from May to September). During the southwest monsoon period, the lower reach of the basin is frequently flooded as recorded from the flood gauge located at Nagalagam Street (refer to Figure 2).
Recurring flood events happened in consecutive years, 2016, 2017, and 2018, at the Lower Kelani River basin. Out of them, the 2016 event was recorded as the most severe flood in almost 30 years (compared to the flood that happened in 1989). Rapid urbanization, unplanned developments, and reduced drainage density were identified as the major causes for these flood events. According to (Sri Lanka post-disaster needs assessment: May 2016 floods and landslides, 2016) records, more than 60% of the total population in Colombo and Gampaha Districts were affected and significant damage to the infrastructure happened due to the 2016 event.

Overall Methodology
It is interesting to observe the flood events under 1D and 2D hydraulic models and then to compare their capabilities. Therefore, Hydraulic Engineering Centre's River System (HEC-RAS) was used to analyse the extreme events that happened in 2016 and 2018 using a combined one-dimensional and two-dimensional (1D/2D) model and two-dimensional (2D) model. In addition, a high-resolution Digital Elevation Model (DEM) was used to map the inundations during the flood events along the river and floodplains.

Elevation and Modification of River Bathymetry
The Digital Elevation Model (DEM) developed using Light Detection and Ranging (LiDAR) was used to capture the topographic features. However, the LiDAR data are unable to detect the terrain features underneath the water [35]. Figure 3a explains this feature. Therefore, measured bathymetries (river cross sections) were fed using a conventional manner. Perpendicular interpolations were conducted for the sections spaced approximately at a 1 km distance, which was also used by Nandalal [3]. These interpolations were important when the change in velocity head was too large to determine the energy gradient along the river. Furthermore, the inline structures (i.e., bridges and salinity barrier) were neglected when modelling the river stretch. In developing the 1D/2D HEC-RAS model, the 2D area was connected to the main river stretch with lateral structures established next to the river over banks. A lateral structure has the same topographical features next to the river over banks. However, sudden variations in the elevation data can lead to model instabilities. As shown in Figure 3b, a smoothening was conducted to the elevations of the lateral structure by following "Savitzky-Golay" algorithm developed by Gorry [36]. However, some places were adjusted manually in order to have a good representation of the natural terrain (see Figure 3b).
The models were developed on a personal computer (HP-Elite Book i5) with an 8 core processor at a speed of 3.

Land Use Characteristics
The upper basin is predominantly covered with heavy vegetation such as forests, grasslands, scrublands, and cultivations, whereas the lower basin is heavily urbanized [37]. The preliminary manning's coefficients were assigned based on the suggestions of Chow [38,39].

Boundary Conditions
An hourly flow hydrograph and a stage hydrograph were established as the upper and lower boundary conditions, respectively, for both models. These unsteady flow hydrographs are given in Figure 4 for years 2016 and 2018. Lower boundary conditions were established at the sea outfall by giving hourly tidal fluctuation (as stage). Therefore, the stage hydrograph at Nagalagam Street was able to replicate the tidal fluctuation including backwater effect [39]. Additionally, an energy slope for the distribution of flow along the boundary line was assigned after measuring an average slope from the DEM for the 2D model.

Implicit Weighting Factor, Calculation Time Step, and Optimal Mesh Size
As was already stated, the Kelani River is a tidally influenced river, and the bed elevation at the outfall is usually lower than the mean sea level. Therefore, the backwater effect has a significant impact on the flows. As the tidal fluctuations are very dynamic, Brunner [40] suggested maintaining an implicit weighting factor around 0.6 in order to have an accurate system. Therefore, a trial and error procedure was applied to identify the optimum value for theta weighting factor. The optimal weighting factor was identified as 0.71 for the Kelani River. The weighting factor is capable of performing stable calculation.
In addition, Courant number (C) was defined to ensure the model stability selection of time step and the mesh size. This Courant number is expressed in Equation 6.
where u is the flow velocity (m/s), ∆t is computational time step (s), and ∆x is the spatial step (m), i.e., cross section spacing in 1D models and cell size in 2D models. According to the Equation (6), the time step of 60 s was selected for the calculation. Minimum of 15 m fine cell resolution in the main channel and maximum of 30 m coarse cell resolution were used for the computational mesh for suitable locations of the river basin.

Water Level Comparison
The better agreement between the simulated water level in Nagalagam Street gauge and the observed water level was ensured with the use of the Nash-Sutcliffe efficiency coefficient (NSE) [41]. The value 1 in the Nash-Sutcliffe efficiency coefficient denotes a perfect match. Equation (7) presents the Nash-Sutcliffe efficiency coefficient.
where O i is observed flow at ith time, S i is simulated flow at ith time, O is mean of observed flow, S is mean of simulated flow, and n is the number of observations. In addition to this, Pearson coefficient of determination (R 2 ) and Root Mean Square Error (RMSE) were calculated to ensure a good capture of the correlation and minimize the error between simulated and observed water levels. These are given in Equations (8) and (9).

Inundation Extent Comparison
A comprehensive evaluation with satellite-observed inundation extents and surveyed inundations were carried out to evaluate the model's performance. Flood Area Index (FAI), accuracy, Bias score, Probability of Detection, false alarm ratio, Probability of False Detection, and Success Index were extracted from Bennett et al. [42], Falter et al. [43], and Khaing et al. [44] in order to compare the inundation extents. A self-explanatory skill matrix composed of Equations (10)-(16) was developed.
Accuracy = hits + correct negatives total (11) Bias score = hits + f alse alarms hits + misses Probability o f Detection (hit rate) = hits hits + misses (13) False alarm ratio = f alse alarms hits + f alse alarms (14) Probability o f False Detection ( f alse alarm rate) = f alse alarms correct negatives + f alse alarms (15) Success Index = 1 2 hits hits + misses + correct negatives correct negatives + f alse alarms (16) where M 1 D 1 is the number of grid cells correctly predicted as flooded by the model (hits), M 1 D 0 is the number of cells flooded in the prediction but observed as dry in the observation (false alarm), M 0 D 1 is the predicted dry area but observed wet area (misses), and the number of correct negatives are the cells predicted and observed as dry.

Calibration and Validation
The calibration of the hydrodynamic model was conducted by following a trial and error procedure while changing Manning's coefficients. The initial guesses for Manning's roughness coefficients were assigned based on the land use type along with visual observations. Suggestions of Chow (1959) were further used for the justification of the values applied for the model.

Comparison of 2016 and 2018 Flood Stages
The year 2016 flood was an extreme event and the stages recorded at Nagalagam Street was observed to be almost at 2.4 m, which was the highest recorded flood after 1989. Figure 5a showcases the stage heights comparison from the observed to simulated stages, while Figure 5b   The time to the peak is also well captured by the 2D model compared to the 1D/2D model. Nevertheless, overestimations, as discussed earlier, can be observed in the 2D model. The performance of the 1D/2D and 2D models with respect to the stages for years 2016 and 2018 were evaluated with statistical parameters suggested by Moriasi et al. [45], as stated in the methodology section. The calculated statistical parameters comparing the simulated stages with the observed stages are given in Table 1. Based on the calculated statistical parameters, the overall model is capable of capturing the flow processes with both models showing a very good performance rating. According to the estimated statistical parameters for both the 1D/2D and 2D models, the 2D model shows a high prediction capability of stages (i.e., NSE = 0.98, R 2 = 0.98, and RMSE = 0.08) compared to the 1D/2D model (i.e., NSE = 0.95, R 2 = 0.91, and RMSE = 0.14). However, the 1D/2D model was comparatively able to capture the flow processes better than the 2D model for the 2018 event. Nevertheless, by comparing the visual characteristics of the hydrograph and the statistical parameters, it can be stated that the 2D model outperforms the 1D/2D model.

Inundation Extent Comparison
The microwave images captured by Sentinel-1 during the floods of 2016 and 2018 were processed and used to compare the simulated inundation extent. The comparisons were conducted by following a cell by cell approach of sub-regions suggested by De Silva et al. [46] and evaluated with the statistical parameters mentioned above. However, the inundation extents of some areas, especially the highly urbanized and mountainous areas, are partly constricted by the satellite images [43]. Figure 6 shows a comparison between the surveyed inundations with the Sentinel-1 satellite-observed inundation extent.
Locations 1-3 were areas that actually experienced floods during the 2016 event, which are marked with the surveyed inundation and further justified with Google Earth pro highresolution images. However, these were not captured by the Sentinel-1 satellite image. The incoming signals of Sentinel-1 could have been double bounced, causing the backscatter to be higher than normal due to the densely built-up areas and, therefore, the signal was unable to reach the ground. This scenario was stated by Solbø and Solheim [47]. This could be a possible reason for not capturing the three locations. Therefore, there can be some underestimations of flood cover in satellite images, even though they are considered one of the most reliable sources to make comparisons with the simulated extent [43]. Figure 7 illustrates the inundated areas from the 2016 and 2018 floods for the river basin using satellite observations and hydrological models (1D/2D and 2D).
Based on the derived inundation maps for the 2016 and 2018 events, the predicted inundation areas from the 1D/2D model were 75.24 km 2 and 34.47 km 2 for the 2016 and 2018 floods, respectively, whereas the predicted inundation areas from the 2D model were 54.19 km 2 and 22.86 km 2 , respectively. Therefore, the 1D/2D model predicts more inundation area compared to the 2D model. Moreover, the inundated areas calculated from the Sentinel-1 satellite images are lower than the inundation extents simulated from both models. These are approximately 14.28 km 2 and 3.51 km 2 , respectively. Most of the flooded areas in the Colombo and Gampaha districts were not captured by the satellite images. However, high-resolution satellite images of Google Earth Engine (refer to Figure 6) showcase the flooded areas that were not captured by the Sentinel-1 satellite images.   The inundated area observed during the 2016 flood event is comparatively larger than in 2018, as it was an extreme event. Apparently, both models have many overlapping areas. In order to have a comprehensive assessment, sub-region-based cell by cell comparison was conducted for both models with the satellite-observed flood, which is summarised in Table 2. For the event in 2016, the 1D/2D model showed a very good agreement with the satellite-observed flood, which is indicated with a probability of detection (hit rate) of 97% and an accuracy of 88%. However, the model tends to overestimate the inundation extent, which can be identified by the false alarm ratio of 48%, a false alarm rate of 13%, and a biased score of 1.97. Nevertheless, the 1D/2D model was able to successfully predict the inundation extent with a success index of 92% for the 2016 event. This is a significantly improved result when compared to past studies for the same basin [42,46].
However, the 2D model was unable to predict the extent as it was in the 1D/2D model for 2016, which is illustrated with a hit rate of 68% and an accuracy of 86%. Furthermore, a similar overestimation was observed when compared to the 1D/2D model, which is quantified with a false alarm ratio of 50%, a false alarm rate of 12%, and a biased score of 1.49. Due to that, the occurrence or non-occurrence of the event was captured with a 78% success index. Therefore, the 2D model for the 2016 event was able to reasonably capture the inundation extent.
However, the 1D/2D model predictions showed a slightly better performance for the 2018 event compared to the 2016 flood event, which is denoted by an increased accuracy of 94%, a reduced bias score (1.35), a lower false alarm ratio (34%), and a lower false alarm rate (5%), but it has a lower hit rate (83%). As a result, the prediction ends up with a success index of 89%, slightly lower when compared to the 2016 event.
The 2D model for the 2018 event showcases similar behaviour as was shown for the 2016 event. It indicated a hit rate of 41%, an accuracy of 89%, a false alarm ratio of 50%, and a false alarm rate of 6%. However, the 2D model underestimates the inundation extent, which is denoted with a bias score of 0.99. Therefore, the success index of the 2D model for the year 2018 is 73%.
After combining the results for both events, the 1D/2D model was able to successfully predict the inundation extent with a high hit rate (90%) and a high accuracy (91%), together with an overestimation denoted by statistical parameters such as false alarm ratio (41%), false alarm rate (9%), and biased score (1.66). Based on the statistical results, the 2D model was unable to predict the inundations in the same way as the 1D/2D model. Even the predictions showed an underestimation when compared with satellite-observed inundation extents.
Even the statistical matrix showed a good agreement: in visual observations, the combined 1D/2D model (refer to Figure 8a,b-pink outline highlighted with red coloured dotted box) seems to be an overestimation. However, during the flood events, those areas were inundated. However, there was no valid proof in the scientific world other than the statements of local people. Therefore, to justify the statement, "during the flood events, those areas were inundated", we conducted a comparison with the recorded maximum discharges at Hanwella Gauge, and the stages at Nagalagam Street with past studies conducted by the Disaster Management Centre Sri Lanka (DMC) and the University of Kelaniya, Sri Lanka [48].
The initial environmental examination conducted by the University of Kelaniya [48] showed some photographs captured during the 2010 flood and the recorded flood surface elevation (2.10 m above MSL), which justify that the area was inundated with the maximum recorded discharge of 907.8 m 3 /s at Hanwella gauge and a stage of 1.58 m at Nagalagam Street gauge (refer to Figure 9).

Summary and Conclusions
The present study was conducted with high-resolution spatial data by two hydraulic models, which are solved in the 1D and 2D environments (1D/2D model and 2D model). The comparisons are conducted with coupled 1D/2D and 2D models of HEC-RAS, which has not been comprehensively studied in terms of capturing the inundations with the flows in the Sri Lankan context. Therefore, the major objective of comparing both models was conducted without changing the topography, channel geometries, boundary conditions, flows, and stages. The calibration was conducted by changing Manning's coefficient and validated for similar flood scenarios that happened closer to the calibrated event. Based on the comparisons conducted, the following conclusions were made.

•
The HEC-RAS 2D model was able to successfully capture the flow process when compared to the coupled 1D/2D HEC-RAS models during high flow conditions. • The 1D/2D HEC-RAS model was a better predictor of flows when it came to low flow situations.

•
Combining the prediction capability of flows during high flows and low flows, the HEC-RAS 1D/2D model is a better predictor than the 2D model.

•
The HEC-RAS coupled 1D/2D model is a better predictor when predicting inundation extents during high flow and low flow situations.

•
Overall, the HEC-RAS 1D/2D model is a better model compared to the 2D model in predicting inundation extents and the flows under high and low flow situations.
The results would be interesting to the stakeholders such as the Disaster Management Centre of Sri Lanka to forecast the future needs based on the capacity of the hydraulic model.