Towards Coupling of 1D and 2D Models for Flood Simulation—A Case Study of Nilwala River Basin, Sri Lanka

Abstract

The Nilwala river basin is prone to frequent flooding during the southwest monsoon and second intermonsoon periods. Several studies have recommended coupling 1D and 2D models for flood modelling as they provide sufficient descriptive information of floodplains with greater computational efficiency. This study aims to couple a 1D hydrological model (HEC-HMS) with a 2D hydraulic model (iRIC) to simulate flooding in the Nilwala river basin. Hourly rainfall and streamflow data of three flood events were used for calibration and validation of HEC-HMS. The model performed exceptionally well considering the Nash–Sutcliffe coefficient, percent bias, and root mean square error. The flood event of May 2017 was simulated on iRIC using the streamflow hydrographs modelled by HEC-HMS. An overall accuracy of 81.5% was attained when the simulated extent was compared with the surveyed flood extent. The accuracy of the simulated flood depth was assessed using the observed water level at Tudawa gauging station, which yielded an NSE of 0.94, PBIAS of −4.28, RMSE of 0.18 and R² of 0.95. Thus, the coupled model provided an accurate estimate of the flood extent and depth and can be further developed for hydrological flood forecasting on a regional scale.

1. Introduction

A significant consequence of climate change is the increased frequency and intensity of extreme weather events [1]. Germanwatch described Sri Lanka as being among the top ten most affected countries due to climate change in 2018, with a global climate risk index of 19.0 [2]. Regional studies conducted on rainfall trends in Sri Lanka have indicated both increasing and decreasing patterns [3,4,5]. It was found that there was a higher tendency for short-duration high-intensity rainfall events to occur within the island [4], which often resulted in water-related natural hazards. Floods are destructive and occur more frequently than other natural disasters [6], with Asia being one of the most affected regions [7]. The Fifth Assessment Report of the IPCC states with medium confidence that the magnitude and frequency of recent floods are comparable to or surpass historical floods [1].

Sri Lanka was recently devastated by floods and landslides due to heavy southwest monsoonal rainfall in May 2017. This natural disaster led to 219 deaths, affecting approximately 230,000 families [8]. The most adversely affected districts by this disaster were Kalutara, Galle, Matara, Hambantota, and Ratnapura. Flooding occurred in towns, villages, and agricultural land bordering the Kalu, Nilwala, and Gin rivers. Matara city is located in the floodplain of the Nilwala river basin, an area prone to flooding at the onset of the southwest monsoonal rains. Matara is the administrative capital of the Matara district and is considered a commercial hub. The risk of flooding is significantly greater in urban areas, such as Matara city where impervious surfaces prevent infiltration and cause runoff to generate quickly [9].

A past study identified precipitation within the basin increases from the south to the north with an increasing elevation [10]. Thus, a significant proportion of the rainfall that contributes to flooding in the lower coastal plains occurs in the basin’s upper region. There is a considerable period of time for extreme precipitation in upper Nilwala to accumulate into floods in the lower coastal floodplains. With the recent advances in information and communication technology, this study aims to assess the feasibility of applying a 1D–2D coupled model to simulate flooding in real-time that could provide sufficient lead time for evacuation and disaster preparation.

The Post Disaster Needs Assessment (PDNA) was conducted following the disaster in 2017 and identified several vital interventions: strengthening of the country’s early warning system, the conduction of hazard risks, and vulnerability assessments were points focused on. Moreover, studies have shown that mortalities caused by floods depend on the scale or severity of the flood, the possibility of an early warning mechanism, and evacuation [11]. Computational models can be used to understand the nature of the flood and simulate hydrological responses to various control measures [12]. They can also generate flood risk maps for planning and developing the resiliency of high-risk regions [13]. The existing early warning system of the Nilwala river basin consists of a network of water gauging stations with extreme water levels categorized as alert, minor flood, and major flood. Flood warnings are issued based on this water level and its trend. There remains room for improvement of the system through flood modeling.

Hydraulic models can be classified dimensionally as 1D, 2D, and 3D models. Three-dimensional models were developed relatively recently, requiring more descriptive parameters that remain unavailable for most basins [14]. The simplest representation of flow in a stream is as a one-dimensional flow in a longitudinal direction. Although these models are computationally efficient, they demonstrate certain limitations; namely the topological discretization of topography as cross-sections of the floodplain [15], the subjectivity of the location and orientation of cross-sections, and the incapability of simulating lateral diffusion of the flood wave [16]. A 2D model solves for water flow equations, both in longitudinal and lateral planes by using either finite element or finite volume analysis methods [17]. They can provide additional detailed information on the floodplain and allow the visualization of flood extent, but are computationally more demanding, requiring greater processing power and storage capacity.

In this study, we coupled a 1D hydrological model to a 2D hydraulic model. The upper portion of the basin is modelled one-dimensionally to carry the flood wave to the two-dimensionally modelled coastal floodplains, which is an urban area with a high risk of flooding. The coupled model will generate comprehensive information of the floodplain with greater computational efficiency [13]. The coupled model allows one to alter the roughness coefficient of the simulation extent, provide a more detailed description of the topography, simulate overland flow, and can be easily visualized as two-dimensional maps [12]. Teng et al., [18] describes several methods to couple a 1D model with a 2D model. These include; modelling a section of the stream in 1D and the rest in 2D, coupling a 1D drainage model with a 2D flood model [19], and coupling a 1D river model with a 2D floodplain model [14,20,21]. This study uses the second approach, where Hydrologic Engineering Center—Hydrologic Modelling System (HEC-HMS) is used to simulate drainage in upper Nilwala and is externally coupled with the 2D hydraulic model, the Nays2DFlood solver, in the International River Interface Cooperative (iRIC).

This approach has been successfully used for flood forecasting and urban flood management. A study conducted in the Brahmani and Baitarani river delta in India by coupling SWAT and SWMM with iRIC yielded satisfactory results and was able to identify the depth and extent of inundation [19]. A similar study using a coupled model conducted in China was similarly successful and capable of assessing the effects of water conservancy structures within the simulation extent [22].

HEC-HMS is a rainfall-runoff hydrological model developed by the United States Army Corps of Engineers used extensively to simulate the hydrological response of catchments [23,24,25]. This model is used to estimate the runoff generated upstream of the 2D flood modeling extent and to provide the boundary conditions required for flood simulation. HEC-HMS was selected in this study due to its versatility in modeling hydrological processes and ability to model storm events. Nays2DFlood solver in iRIC is a 2D hydraulic model used for flood inundation mapping [19,26,27,28,29]. The iRIC is a recent hydraulic model that can model large flood events [30]. Owing to the nature of its recency, it has not been applied widely, but was able to simulate flooding in ungauged catchments where data are scarce. It is a free open-source software developed by Yasuyuki Shimizu and Toshiki Iwasaki of Hokkaido University, Japan.

2. Materials and Methods

2.1. Study Area

Figure 1 depicts the Nilwala river basin located along the southern coast of Sri Lanka. The river is 72 km in length and originates at Panilkanda, north-east of Deniyaya, at 1050 m above mean sea level [10]. The total drainage area of the river basin is 1025 km² [31]. Steep longitudinal slopes characterize the basin’s upper region, while the lower basin consists of gentle slopes [31,32]. It passes several small towns such as Deniyaya, Morawaka, and Akuressa and flows into the ocean Thotamuna close to Matara city. The catchment is located in the wet zone, with most of the upper catchment covered by rainforests. A significant portion of the basin is used for cultivating agricultural crops such as paddy, coconut, tea, rubber, cinnamon, and citronella [32,33].

2.2. Data

Hourly rainfall and streamflow data required for simulation were obtained from the Irrigation Department of Sri Lanka. The Nilwala river has three streamflow gauging stations; Urawa, Pitabeddara, and Panadugama, from which hourly streamflow data were collected. Continuous data of the 2017 flood event were not available for Urawa and Panadugama stations. Thus, calibration and validation of HEC-HMS was undertaken solely using data obtained from Pitabeddara. Hourly precipitation data from four stations in upper Nilwala; Deniyaya, Panadugama, Pitabeddara, and Urawa were additionally used. The data used in this study are summarized in

Table 1. Summary of data used in study.

Data Type	Source	Frequency	Period/Scale
Precipitation	Irrigation Department	Hourly	2017–2019
Streamflow	Irrigation Department	Hourly	2017–2019
Water level	Irrigation Department	Hourly	May–June 2017
Digital Elevation Model (DEM)	United States Geological Survey (USGS)	-	30 m

.

2.3. Rainfall-Runoff Modelling Using HEC-HMS

The Hydrologic Engineering Center—Hydrologic Modelling System (HEC-HMS) model was designed by the US Army Corps of Engineers as a rainfall-runoff model to determine the hydrological response of watersheds [34]. It was developed for dendritic watersheds and has a wide array of applications, including urban drainage, flood-loss reduction, flood early warning system planning, reservoir spillway design, stream restoration, flow forecasting, surface erosion, sediment routing, and systems operation [35].

The Nilwala basin was delineated and divided into 20 subbasins using the HEC-GeoHMS extension in ArcGIS. Terrain preprocessing and basin processing tools were used to generate the basin model file containing the drainage network and delineated basin. Next, precipitation was estimated using the Thiessen polygon method [36], and suitable weights for each subbasin were defined to create a meteorological model. The basin and meteorological models were imported into HEC-HMS. Hourly precipitation data from the Deniyaya, Panadugama, Pitabeddara, and Urawa gauging stations, and observed discharge data for the selected events were imported into the HEC-HMS interface.

HEC-HMS allows for the separate modelling of hydrological processes; loss, transformation, baseflow, and routing with several models for each process. The selection of the model for each process should be based on the catchment characteristics, data availability, and whether the simulation is event-based or continuous [25].

The Snyder unit hydrograph method [37] was selected to simulate the transform in conjunction with the deficit and constant loss method, and the recession baseflow model [38]. The lag time, $t_p$ to be used in the Snyder unit hydrograph was defined as the time period between the centroid of excess rainfall and the peak discharge. It was calculated for each subbasin using the following relationship;

$$t_p=C{C}_t{\left(L{L}_c\right)}^{0.3},$$

(1)

where $C_t$ = basin coefficient, $L$ = length of the mainstream from the outlet to the basin divide, $L_c$ = length along the mainstream from the outlet nearest the watershed centroid and $C$ = a conversion factor (0.75 for SI units).

The Muskingum–Cunge method [39] was used to model routing in all reaches with a Manning’s n coefficient of 0.040 [40,41]. Parameters such as reach length and slope were estimated using topographical maps. A deterministic optimization approach was followed with NSE set as the objective function. Optimization trials were carried out for the deficit and constant loss model, and the recession baseflow model to determine proper parameters.

Table 2. Parameters optimized in HEC-HMS model.

Parameter	Units	Optimized Value	Model
Initial deficit	mm	3.62	Deficit and constant loss
Maximum storage	mm	67.575	Deficit and constant loss
Constant rate	mm/h	6.135	Deficit and constant loss
Impervious	%	30	Deficit and constant loss
Initial discharge	m3/s/km2	0.1	Recession baseflow
Recession constant	dimensionless	0.45	Recession baseflow
Ratio to peak	dimensionless	0.25	Recession baseflow

lists the parameters optimized in HEC-HMS.

The model was initially applied to the basin’s upper region reaching up to Pitabeddara due to the lack of reliable streamflow data below this point. Once the model was calibrated and validated, it was then extended to the entire basin. For the simulation, three flood events were selected: May 2017, May 2018, and September 2019. The May 2018 and September 2019 events were used for calibration, while the May 2017 event was used to validate the model. The two selected events for calibration had both reported flooding, although not at the magnitude of the May 2017 flooding event.

2.4. Hydrodynamic Modelling Using Nays2DFlood Solver in iRIC

The International River Interface Cooperative (iRIC) is a free numerical simulation software consisting of various computational solvers to decipher hydrological problems. Nays2DFlood is one such solver in iRIC, which is used for flood simulation. The iRIC software platform was initiated by Professor Yasuyuki Shimizu (Hokkaido University) and Dr. Jonathan Nelson (USGS). Nays2DFlood relies on unsteady 2-dimensional plane flow simulations using boundary-fitted coordinates as the general curvilinear coordinates [42]. This solver can be applied to river basins with sparse data, especially in developing regions as it only requires the DEM without specific topographical dimensions of the river.

The equation of continuity (Equation (2)) and equations of motion (Equations (3) and (4)) of two-dimensional unsteady flows in a rectangular coordinate system are expressed as;

$$\frac{\partial h}{\partial t}+\frac{\partial \left(hu\right)}{\partial x}+\frac{\partial \left(hv\right)}{dy}=q+r,$$

(2)

$$\frac{\partial \left(hu\right)}{\partial t}+\frac{\partial \left(h{u}^2\right)}{\partial x}+\frac{d\left(huv\right)}{\partial y}=-hg\frac{\partial H}{\partial x}-\frac{{\tau }_x}{\rho }+D^x,$$

(3)

$$\frac{\partial \left(hv\right)}{\partial t}+\frac{\partial \left(huv\right)}{\partial x}+\frac{d\left(h{v}^2\right)}{\partial y}=-hg\frac{\partial H}{\partial x}-\frac{{\tau }_y}{\rho }+D^y,$$

(4)

where,

$$\frac{{\tau }_x}{\rho }=C_{f}u\sqrt{u^2+v^2},$$

(5)

$$\frac{{\tau }_y}{\rho }=C_{f}v\sqrt{u^2+v^2},$$

(6)

$$D^x=\frac{\partial }{\partial x}\left[v_t\frac{\partial \left(hu\right)}{\partial x}\right]+\frac{\partial }{\partial y}\left[v_t\frac{\partial \left(hu\right)}{\partial y}\right],$$

(7)

$$D^y=\frac{\partial }{\partial x}\left[v_t\frac{\partial \left(hv\right)}{\partial x}\right]+\frac{\partial }{\partial y}\left[v_t\frac{\partial \left(hv\right)}{\partial y}\right],$$

(8)

where, $h$ = water depth, $t$ = time, $u$ = flow velocity in $x$ direction, $v$ = flow velocity in $y$ direction, $g$ = gravitational acceleration, $H$ = water surface elevation, ${\tau }_x$ and ${\tau }_y$ are riverbed shear stress in $x$ and $y$ directions, respectively, $C_f$ = riverbed friction coefficient, $v_t$ = eddy viscosity coefficient, $\rho$ = density of water, $q$ = inflow through a box culvert, a sluice pipe or a pump per unit area, and $r$ = rainfall [42]. These basic equations, which present in a rectangular coordinate system, are mapped to a generalized curvilinear coordinate system. The shape of the simulation mesh can be changed by accomplishing this depending on the boundary conditions [42].

The inputs required for iRIC are the DEM, satellite images, and streamflow hydrographs at the inlet points. The 30 m SRTM DEM from the United States Geological Survey (USGS) was used for the simulation, while the required background satellite images were obtained from Google Earth. Next, a calculation grid with 80 m grid size was created by defining a grid centerline representing the river’s flow. Values for n_I, n_J, and W were defined, where n_I = number of divisions in the longitudinal direction, n_J = number of divisions in the transverse division, and W = grid width in the transverse direction. Subsequently, calculation conditions were set, including boundary conditions, initial conditions, roughness conditions, and an appropriate time step for computation. The inflow streams identified earlier were then located on the grid. The discharge hydrographs obtained from the HEC-HMS model were then imported into iRIC manually and set as boundary conditions.

Impermeable objects such as roads, banks, embankments, and other flood control structures were located on the grid as obstacles. Matara city was outlined as an area of the grid occupied by buildings, and its areal fraction occupied by buildings was defined. Most calculation conditions of iRIC were left as default with a numerical upwind scheme used as the finite differential method. The northern, eastern and western boundaries were set as inflow boundaries while the southern boundary was modelled as having free outflow. The resulting inundation map’s accuracy was determined by comparison with a survey of the flood extent carried out by the Irrigation Department. The depth of flooding was validated using observed water level data at Tudawa gauging station.

2.5. Statistical Evaluation

Both models were statistically evaluated using Nash–Sutcliffe efficiency (NSE) [43], percent bias (PBIAS) [44], and root mean square error (RMSE).

The NSE indicates the goodness of fit of simulated data to observed data. It is calculated using the following equation, given by Nash and Sutcliffe (1970).

$$\mathrm{NSE}=1-\frac{{\sum }_{i=1}^n{\left(O_i-S_i\right)}^2}{{\sum }_{i=1}^J{(O_i-\overline{O})}^2},$$

(9)

where $O_i$ = value of $i$th observation, $S_i$ = value of $i$th simulation, $\overline{O}$ = mean of observations and $n$ = total observations. NSE values range from −∞ to 1, with an NSE of 1 indicating a perfect fit.

The PBIAS measures a simulated data’s tendency to be greater than or lesser than the observed values [44]. A negative value of PBIAS indicates a model overestimation bias, while a positive value indicates a model underestimation bias. A PBIAS value of 0 is optimal.

$$\mathrm{PBIAS}=\frac{{\sum }_{i=1}^{n}100\left(O_i-S_i\right)}{{\sum }_{i=1}^J{O}_i},$$

(10)

The RMSE was performed for quantifying the prediction error in terms of the units of the variable calculated by the model [45].

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^n{(O_i-S_i)}^2}{n}},$$

(11)

3. Results and Discussion

3.1. Calibration and Validation of HEC-HMS

The Irrigation Department has only recently begun collecting hourly streamflow and rainfall data from gauging stations at Panadugama, Pitabeddara, and Urawa. Thus, only three storm events were found with significant precipitation that resulted in downstream flooding. Of the three selected events, the May 2018 and September 2019 events were selected for model calibration. The flooding event of May 2017, considerably more intense among the events, was selected for validation and subsequently modelled using iRIC.

The Snyder unit hydrograph method for transformation produced excellent results during calibration and validation processes. A past study conducted within the Nilwala river basin confirmed that the Snyder transformation method provided the most accurate results [46]. Another study conducted in the Attanagalu Oya basin using HEC-HMS concluded that the Snyder unit hydrograph simulated flow with greater reliability than the Clark unit hydrograph and SCS Curve Number methods [23]. Moreover, the study recommends applying this model in conjunction with the deficit and constant loss model to simulate river flow in the country’s wet zone [23].

The recession baseflow model in HEC-HMS simulates a catchment behavior when flow recedes exponentially and was designed primarily for event simulations [34]. Recession baseflow and the Muskingum–Cunge routing method were successfully used to predict flow in the Nilwala river basin [46]. The downstream regions of the basin display a characteristic gentle slope, modelled more accurately using the Muskingum–Cunge method [47]. Detailed river cross-sectional data could not be obtained for the study, and thus the parameters for the routing model were estimated using the basin model generated by HEC-GeoHMS and satellite imagery. Moreover, obtaining 8-point cross-sectional data could significantly improve model performance. Alternatively, a trapezoidal shape design was specified for the routing channels.

During optimization, it was understood that the loss parameters influenced the peak of discharge, while the baseflow parameters affected the initial flow and recession curve. The transformation and routing processes influenced the time of the peak. After several optimization trials, the model was adequately calibrated for selected flood events.

Figure 2a depicts the observed and simulated hydrographs for the May 2018 event, while Figure 2b shows the September 2019 flood event. The May 2018 event provided the values of 0.93 NSE and 0.3 RMSE (

Table 3. Statistical evaluation HEC-HMS calibration and validation steps.

Title 1	Event	NSE	PBIAS	RMSE
Calibration	May 2018	0.93	−10.95%	0.3
	September 2019	0.746	−4.82%	0.5
Validation	May 2017	0.927	−8.33%	0.3

), categorized as being excellent according to the criteria suggested by Moriasi et al., [48] for monthly hydrological simulations. The September 2019 event produced an NSE value of 0.746 after calibration and can be categorized as useful while having an RMSE value of 0.5. The PBIAS for both events was negative, −10.95% and −4.82%, respectively, indicating an overestimation bias.

The observed and simulated hydrograph at Pitabeddara gauging stations for the May 2017 event is shown in Figure 3. The calibrated parameters were averaged and validated using the May 2017 event. As mentioned earlier, this event is a comparatively more severe flood event, and rainfall frequency analysis has classified this event as having a return interval of 75–100 years [49]. Due to the severity, the gauge station at Pitabeddara was completely flooded during the event’s peak. The peak discharge of the flood was estimated by flood marks. Thus, there remains considerable uncertainty in the peak discharge of this event. However, the simulated flow using averaged parameters showed the goodness of fit with the observed flow with an NSE value of 0.927, a PBIAS of −8.33%, and an RMSE value of 0.3. As with calibration, the PBIAS indicates a model overestimation bias.

3.2. Hydraulic Model Simulation

Initial simulation runs identified that the model could not correctly identify the flow path of the mainstream in specific regions due to the coarser 30 m DEM used. As a result, the elevation file was modified to indicate the river’s correct flow path. A study conducted for flood simulation in Kabul using iRIC, Shokory et al., [28] modified the DEM similarly as the flood would not pass over the original DEM. Initially, the model was applied to a smaller region, including Matara and the region immediately north of the city. This overestimated the city’s flood as the model failed to account for flooding that occurred above this region. Thus, the simulation extent was extended further upstream up to Akuressa for more accurate simulation.

Six inlets were located within the simulation extent (Figure 4a). The hydrographs in Figure 4b, which were obtained using the HEC-HMS model were set as the boundary conditions for Nays2DFlood at each inlet. The highest contribution was seen from inlet 1, the main river, and included flow generated from a larger proportion of the upper catchment. Inlet 2, the second-largest contributor, represents a tributary of the Nilwala river, Kirama Oya.

Paddy cultivation represents the primary land use of the region surrounding Matara city. The 30 m DEM could not recognize the complicated network of interconnected minor canals that characterize this region. The roughness coefficient was taken as 0.15 for forested areas and 0.09 [46] for other areas based on land cover and by reviewing past studies and literature [41]. Moreover, a section of Matara city was incorrectly simulated as being flooded when in fact flooding did not occur. We attribute this to faults in the DEM, that was unable to detect flood control structures such as bunds on either side of the river due to topographical changes that might have occurred owing to its production and low resolution.

The extent of flooding was verified using a survey of the flood extent carried out by the Department of Irrigation for the 2017 May flood event (Figure 5). The survey was carried out by considering flood marks and accounts of the residents in the months shortly after the flood event. An accuracy assessment (

Table 4. Accuracy assessment of flood extent simulation using iRIC.

Type	Inundated/Unaffected	Accuracy (%)
Producer’s Accuracy	Inundated	70.10
	Unaffected	88.72
User’s Accuracy	Inundated	79.76
	Unaffected	82.40
Overall Accuracy		81.50

) using 10,000 random points was carried out to compare the simulated flood extent with the surveyed extent, with a maximum simulated depth greater than 0.5 m being considered as a flood, based on the classification of Hamdan et al., [50]. The producer’s accuracy details the number of reference points (i.e., surveyed points) classified correctly to the total number of reference points for that particular class. The user’s accuracy represents the correctly simulated points to the total number of points simulated for its respective class. The results of the assessment indicated an overall accuracy of 81.5%. The relatively lower producer’s accuracy for inundated class (70.1%) when compared with the higher user’s accuracy (81.5%) indicates a higher error of observed inundated points being simulated as unaffected.

The simulated flood depth accuracy was verified using observed water levels at Tudawa station, located north of Matara city. This resulted in; NSE = 0.94, PBIAS = −4.28, RMSE = 0.18 and R² = 0.95 (Figure 6).

In Figure 7a, depths over 7 m were observed towards the northern aspect of the simulation extent close to the main river’s inflow points and its tributary Kirama Oya (i.e., inlets 1 and 2, from Figure 4). There remains uncertainty in flood depth in the upper regions as the model fails to consider any inundation that may have occurred upstream of the simulation extent which is modelled as purely one-dimensional. Furthermore, the conservation of momentum flow and the return flow of water at the linkage sites were neglected [18]. However, when considering the region closer to Matara, a maximum depth of 5 m was simulated. Comparing the velocity magnitude map (Figure 7b) with the background image affirms that the model could correctly identify the river’s flow path.

Shokory et al., [28] studied the influence of rainfall and streamflow on flooding. The study found that rainfall did not significantly impact flooding in the Kabul case study, with the difference between the two resulting in the formation of “puddles”. The Nilwala basin showed similar characteristics as the flood occurs mainly due to heavy precipitation in the upper catchment [10]. This can be clearly understood when comparing the input hydrographs generated using HEC-HMS for the May 2017 event (Figure 4). In addition, a rainfall gauging station with the capability to measure hourly precipitation within the simulated extent was unavailable. Due to these reasons, the rainfall within the basin was omitted from the simulation.

As this study seeks to develop a tool that could be applied for real-time flood forecasting, the time required to simulate needs to be considered. The grid size of the generated mesh will affect the simulation time where a grid with higher resolution will produce more accurate results but will require a significantly longer time to complete the simulation. Thus, an 80 m grid size was selected to produce results faster with acceptable accuracy.

During the May 2017 flood event, it took approximately 36 h for the flood peak to reach Tudawa. The average time taken to run the hydraulic model was close to 5 h, but this could be further reduced by using a machine with greater computational power. However, the model did simulate within an acceptable period, with sufficient lead time for taking necessary precautions based on the generated flood risk maps.

A flood risk map (Figure 8) was generated based on the May 2017 flood event. Regions were classified as high, medium and low risk according to the flood depth [50]. Such maps generated for different return periods can aid in the design and planning of infrastructure, flood control structures, and provide guidance in community resiliency development projects. In addition, maps generated in real-time during a storm event can be used to issue flood warnings and provide specific instructions for evacuation where needed.

4. Conclusions

The Nilwala river basin’s upper region was modelled in HEC-HMS using the Snyder unit hydrograph method to simulate transformation, while the Muskingum–Cunge routing model, deficit and constant loss model, and recession baseflow model were used to model routing, loss and baseflow, respectively. The model was calibrated and validated using three flood events: May 2018 and September 2019 events for calibration and the May 2017 event for validation. The model performed well considering the Nash–Sutcliffe efficiency, percent bias, and root mean square error. Simulated stream hydrographs from HEC-HMS were externally coupled with iRIC to simulate inundation in Matara for the 2017 May flood event. The neglection of upstream flooding, momentum transfer and return flow at the linkage site presented major limitations of the model. The accuracy of the DEM used will greatly impact the performance of iRIC. However, the extent of flooding was estimated, with an overall accuracy of 81.50%. The accuracy of simulated flood depth was verified using water levels at the Tudawa gauging station which resulted in an NSE of 0.94, a PBIAS of −4.28, an RMSE of 0.18 and an R² of 0.95. Based on these values, it can be concluded that the 1D–2D coupled model can accurately identify the extent of flooding in addition to the flood depth. Nevertheless, the model’s accuracy should be further improved by simulating more events with available observed flood data. Moreover, the model should be tested in the field using real-time flood events to evaluate its feasibility as a tool for flood forecasting and early warning systems. Accordingly, a flood risk map was developed for the May 2017 flood event. Further studies should be undertaken to develop such maps for storms with different return periods and for various climate change scenarios.