Two-Dimensional Hydrodynamic Modelling of Flood Inundation for a Part of the Mekong River with TELEMAC-2D

Aims: This paper presents a study on the development of a 2-dimensional (2D) hydrodynamic model based on TELEMAC-2D for the flood simulation of the river from Kratie to Kampong Cham in Cambodia, a part of the Mekong River. The motivation behind the research was to study the feasibility of TELEMAC-2D in flood forecasting, and specifically to determine its adequacy in flood simulations with a focus on the reduction in model run-time through parallelization. Place and Duration of Study: DHI-NTU Centre, Nanyang Environment and Water Research Institute (NEWRI), Nanyang Technological University, between November 2013 and March 2014. Methodology: We simulated an actual flood event which occurred between June to November in 2001 for the stretch of the Mekong River from Kratie to Kampong Cham and compared the model simulations with MODIS satellite Images for specific days in the pre-, peakand post-flood period. Results: It was found that during the peak-flood period, there was high percentage (> 90%) match between the simulation results and observation obtained from satellite images while the match was below 50% for the preand postflood periods. Conclusion: The 2D simulation results were consistent with observations from satellite imaging. Original Research Article Vu et al.; BJECC, 5(2): 162-175, 2015; Article no.BJECC.2015.013 163 The discrepancy at preand post-flood may be due to the fact that (i) the model takes into account only hydrodynamic processes of flows in the river and flood plain, it does not consider other hydrological processes such as infiltration or evaporation which may be important during the preand postflood periods, and (ii) the resolution of MODIS satellite image at 500m x 500m may be too coarse and therefore not sufficient to identify flooded areas when the area is small or water depth low. Finally, it was found that the computing time can be reduced significantly with parallelization using multi-core processors, albeit with lesser advantage in speedup when the number of cores increased beyond 4.


INTRODUCTION
The Mekong River is the 7 th longest river in Asia, defined by its length and discharge [1]. It has a huge impact to the livelihoods of the people of Cambodia, Laos and Southern Vietnam, which rely heavily on it for socio-economic benefits arising from fishery and agricultural farming activities. Flooding in the Mekong River generally occurs during the summer flood season (i.e. from May to October) due to overflows from the river branches and channels that spread into flood plains. Additionally, the low elevation of the floodplain allows the river to be more prone to overflow from the increased water levels. Tidal motion of the sea may also cause exacerbate the flooding problem. In 2000, it was discovered that tidal backwater obstructed the draining of high water levels in the Mekong Delta, inevitably resulting in floods [2]. Annual floods with normal magnitude are known to bring benefits to the agriculture activities and provide water resources and nutrient rich sediments. Floods therefore provide abundant benefits to the riparian countries of Lao-PDR, Cambodia, Thailand and Vietnam, in terms of agriculture, fisheries, and the creation and maintenance of wetlands. However, when extreme floods occur, whereby water levels exceed the critical beneficial level, and flooding is prolonged, this results in damage to land and property, as well as loss of lives. One of the worst flood events in the Mekong Basin occurred in the year 2000, when flood levels exceeded more than 3 meters. Over 600 communities throughout the river basins were inundated impacted, with 11 million people affected by the flood [3]. The flood also resulted in more than 800 human casualties [4]. Around 160 thousand hectares (ha) of rice crops were inundated and had to be harvested immediately and more than 55 thousand ha of rice crops were completely destroyed [5]. Therefore, studies on flood inundation with numerical models can provide useful information that allows efficient flood management for the mitigation of flood damage.
In principle, overland flows caused by over flows from river branches and canals and spread widely on flood plains are usually considered as shallow water waves which can be described mathematically by the 2-dimentional (2D) Saint-Venant equations. In this study the Saint-Venant equations (also Shallow Water Equations) were solved by the Finite Element Method in TELEMAC-2D [6], an open source hydrodynamic modeling software. The implementation of a 2D hydrodynamic numerical model for a large river system and flood plains based on the shallow water equations is data intensive and requires extensive computing facilities. Normally, in order to provide an accurate assessment of the development of a flood, the numerical models have to handle an enormous amount of information, typically involving the spatialtemporal characteristics of flood levels, regional geography, as well as infrastructure such as buildings and the drainage system. The volume of information and the computational complexity also increase when a large area is to be modelled or when a higher resolution is needed. Indeed, a flood modelling system which is deployed by sequential simulations is unlikely to be performed completely in a reasonable timeframe. In the present study, the TELEMAC-2D software was employed using a parallel processing approach based on the domain decomposition method in combination with the message passing interface protocol. This decreased the simulation time on multi-core processor computers [7].
In this paper, we first provide the background information about the model area as well as the TELEMAC-2D model, and then present the methodology for developing the 2D hydrodynamic model to investigate the flood event in 2001 for the model area. We also offer the comparison of the simulation results with satellite images of the inundation. Finally, we discuss the efficiency and stability for running TELEMAC-2D model.

Model Area
The Mekong River can be divided into two parts -1). The Upper Basin, also known as Lancang Jiang, located in China and 2) the Lower Mekong Basin (LMB) from Yunnan in China to the South China Sea. Fig. 1 shows the Mekong Basin and its main tributaries. The Upper Basin i.e. the yellow area in Fig. 1, only contributes 15 to 20 percent of the total flow in the Mekong River, with a majority of the flow contribution arising from the LMB, shown in blue in Fig. 1. The LMB contains two mains areas i.e. the middle Mekong (between China and Northern Cambodia) and the Mekong delta (the region of the river before discharge into the South China Sea).

Fig. 1. Hydrographic map of the Mekong basin (Source: MRC [2])
The study area is in Cambodia, from Kratie to Kampong Cham and consists of a single river route between the two locations. According to Mekong River Commission [2], the Cambodian floodplain takes up about 20% of the entire Mekong Basin, approximately 155,000 km 2 in basin area, and flow in this region contributes about 18% of the total flow in the Mekong River. As the Mekong River reaches Kratie, the topography of the LMB system changes abruptly. Upstream from Kratie, the river flows through mountainous or highland areas with a clearly identifiable mainstream channel. The hydrodynamics of the river system from Kratie to the river mouth is influenced by three main factors: the inflow from regions upstream from Kratie, the storage of the Great Lake (Tonle Sap in Cambodia) which connects via the Tonle Sap River to the Mekong at its confluence at Phnom Penh, and the tide level at the river mouth in Vietnam. The Great Lake functions as a natural storage reservoir for this region. It stores a large amount of flood water during the wet season (June to November) which is slowly released during the dry season (December to May), and acts as a crucial source of water supply to the delta. In the Mekong Delta, the river branches are linked together by manmade channels forming an extremely complex network. The tide effect is significant at stations near the river mouth in Vietnam to Phnom Penh in Cambodia. The details can be found in Kummu et al. [8].
This study was carried out as part of a larger study for real-time flood management for the region. The area from Kratie to Kampong Cham, with a river segment of 113 km in length and flood plain of 2,800 km 2 in size was selected. The model was used to simulate a flood event in 2001 which was comparable with another major flood event occurring in 2000. The maximum flood extent in 2000 was monitored by the Mekong River Commission [9] and is shown in Fig. 2, where the flooded area was approximately 45,000 km 2 .

TELEMAC-2D
The TELEMAC system which was introduced by the National Hydraulics and Environmental Laboratory -a part of the R&D group of Electricite de France, is a numerical modelling system developed to study environmental processes in free surface transient flows. The primary purpose of the TELEMAC-2D model is to simulate the dynamics of flow in a water-body, including rivers, estuaries or coastal areas via the solution of Shallow Water Equations (SWEs). TELEMAC-2D model solves the depth-integrated SWEs where the horizontal length scale is much greater than the vertical, using the finite element method on an irregular mesh that covers the floodplain. At each computational node of the mesh, TELEMAC estimates the depth of flow and two horizontal velocity components [7,10,11]. Effects of bottom friction, diffusivity and Coriolis and hydraulic structures such as weirs are considered in the model. The governing equations considered are: where h is the depth of water (m); u, v are velocity components (m/s); g is the gravitational acceleration (m/s 2 ); Z is the free surface elevation (m); t is time (s); x, y are horizontal space coordinates (m); S h is the source or sink of fluid (m/s) and S x , S y are source or sink terms in dynamic equation (m/s 2 ) and ∇ ⃗ is divergence.
The TELEMAC-2D model can be used for both steady and non-steady simulations taking into account non-linear effects including propagation of long waves. TELEMAC-2D supports parallel processing which can be performed on a multicore workstation to improve computational speed, using the technique of domaindecomposition where the computational grid or mesh is split into a series of smaller subdomains. Each sub-domain is run by individual processors and the architecture therefore requires information exchange activity between computational nodes i.e. processors operate independently and merely pass information between one another, such as fluxes at the boundary of each sub-domain. Parallel virtual machine is used as the messages-passing interface (MPI) to account for the communication between computation cores.

Model Developments
The topography data used in this study consists of the Digital Elevaltion Model (DEM) with resolution of 200 m for the flood plain and 42 cross-sections for the river segment from Kratie to Kampong Cham. In the pre-processing stage, the cross-section data were interpolated and combined with the DEM of the flood plain to obtain an overall topographical grid for the area, with resolution of 200 m and each cell consisting of information of the x-and y-coodinates and eleveation z. The model domain boundary was estimated based on the maximum flood extent of historical flood events. The entire modelling domain was discretized using a system of irregular triangular elements combining the mesh for the channel and flood plains. The mesh for the channel was created at intervals of 250 m and the mesh for the flood plain was created with a default edge length of 450 m. The elevation z was assigned to nodes on the mesh using the nearest neighborhood method. A 3D view of the mesh adopted for this study is shown in Fig. 3.

MODIS satellite image processing
Satellite images of the case study area were obtained for seven days selected through the simulation period to depict pre-, peak-and postflooding. The Moderate-Resolution Imaging Spectroradiometer (MODIS) images from National Aeronautics and Space Administration (NASA), which are available online, were used. The MODIS images used for our project are MOD09A1 -MODIS Terra/Aqua Surface Reflectance 8-Day L3 Global 500 m images. Each image is a composite of eight consecutive days of surface reflectance images at 500 m resolution. Cloud coverage is reduced through using these 8 days of images, with the 'best' observation for every cell in the image retained to produce the MODIS image. The 'best' observation is selected based on criteria such as high observation coverage, low view angle, absence of clouds or cloud shadow, and aerosol loading (MODIS Surface Reflectance User's Guide, Version 1.3, 2011). Each MODIS satellite image contains seven spectral bands of data, with wavelengths ranging from 620nm to 2155 nm.
The Normalized Difference Water Index (NDWI) is applicable for remote sensing of waterbodies from space. It has been used for flood detection [12,13]. The NDWI is given by: where Green is the Green Band, Band 2, while NIR is a Near-Infra Red Band, which can be taken as Band 4 from the MODIS satellite spectral bands. Both of these Bands are able to sense similar depths through vegetation canopies. While liquid absorption by the Green Band is negligible, there is weak liquid absorption by the NIR Band, which is also enhanced by canopy scattering. Hence, the NDWI is able to clearly show changes in liquid water content of vegetation canopies [12]. The converted Band 2 and Band 4 from MODIS images were used in the computation of NDWI to compare with our model results. When the value of NDWI is more than 0, it indicates a flooded area. Fig. 4 shows an example of the processed satellite image and its segregation between areas with NDWI values more than or less than 0. The lighter areas have values less than 0, indicating low or no flooding. On the other hand, the darker areas have NDWI values more than 0, indicating inundation. The case study area of Kratie to Kampong Cham is located within the boxed area.

ArcGIS comparison
The TELEMAC simulation results were postprocessed in Blue Kenue (Canadian Hydraulics Centre). Water Depth of the simulation results for each of the same seven days were saved as shape files that can be read by ArcGIS. Since eight results were printed out for each day of the simulation, we used the results extracted for Midday (12:00:00) for comparison.
Several threshold values were used to determine the percentage of agreement between the NDWI and simulated water depth. For the NDWI values, we checked if NDWI exceeded 0, -0.1, -0.2, and -0.3; whereas for the model water depth results, we identified cells that exceeded 0, 0.3 and 0.5 in value. By matching the cells that give both NDWI and water depth exceeding the threshold values, we can then conclude whether or not the model simulation and satellite images give good agreement on flood occurrence. Figs. 5 and 6 show an example of the model result and satellite data for the flood peak on 21 Aug 2001, respectively. Fig. 7 shows the superimposed layers when NDWI > -0.3 and water depth > 0.5.

Simulation Results
The

Comparison between Simulation Results and MODIS Satellite Images
The model results and satellite information were separately analyzed and the flooded cells, identified by both the model and satellite image were compared. A match between model results and satellite image is determined in ArcGIS when a cell is identified to be flooded by both the model and satellite image. The matched cells are shown in dark blue in Fig. 9 for areas that are inundated.
In order to provide a quantitative evaluation of the agreement between results from these two approaches, an index of percentage of agreement was used: where N is the number of matched cells with NDWI > -0. 3  The threshold values were varied to determine the best conditions for optimal match between NDWI and simulated water depth results. For NDWI values obtained from MODIS data, we used threshold values of 0 to -0.3, with values exceeding -0.3 giving more areas with inundation, compared to NDWI > 0. Similarly, for modelled water depth results, water depth more than 0 showed greater inundated areas with water depths > 0.5 m. With the different combinations of flood threshold values therefore, we are able to determine the optimum combination for agreement between the pre-, peak-, and post-flood periods. Table 1 shows that the highest percentage of matched results for the three periods considered occurred when values of NDWI > -0.3 and model water depth > 0.5 m. With the lowest threshold of -0.3 for MODIS data, we extracted the maximum wetted area possible. The threshold for model water depth was set at 0.5, to extract the extreme flooded areas. Hence, this combination would give the highest agreement between both the model and satellite image.
The optimal matching combination can also be observed in Figs With this threshold, we explore the water depth threshold, which would yield the highest percentage match.
It can be seen that the percentage match for the pre-and post-flood are < 50%. However, for the peak period, agreement between NDWI and modelled water depth results are high and are about 90%. The low match during the pre-and post-flood period may be attributed to the low water levels in the river that could have been undetected by MODIS satellite images. As the water levels were too low in the river when there was still minimal flooding, the areas with water did not show up clearly on the satellite imagery and hence could not be extracted on ArcGIS. Another possible reason for the observed discrepancy is that the model takes into account only hydrodynamic processes of flows in the river and flood plain and it does not consider other hydrological processes such as infiltration or evaporation. Thus, at times before and after the flood, the model tends to provide larger flood inundation area than that the actual case since the abstractions represented by these processes are not taken into account.  August, the percentage comparison between the MODIS satellite image and modelled water depth was high. Once again, the errors could be attributed to losses not considered in the model or the shallow depths at the fringes of the inundated zone. The high percentage match between model simulation results and MODIS data indicates that TELEMAC-2D is indeed suitable for flood studies; however, issues related to the poor match between model and observations during the pre-and post-flood periods will have to be further investigated. Nevertheless, it can be concluded that for the period around the peak flood, the 2D hydrodynamic model from TELEMAC-2D is able to produce results that are in close agreement with observations.

Model Computing time with Multi-Core Processor
The simulation for the flood event from 4 th June to 30 th November 2001 was implemented on a multi-core workstation. The computing time was 177, 87, 51, 44 minutes corresponding to 1, 2, 4, and 6 cores used as shown in Fig. 13.
As expected, the computing time was lower when the number of cores used increased; however, the improvement in computational speed was less significant when the number of cores was larger than 4. The decrease in computing time was about 50.8 % when 2 cores were used. The reduction in runtime was 41.3% when 4 cores were used. Finally, it was only 13.7% when 6 cores were used. This trend reveals that the computing time approaches a limitation when the number of cores is more than 6. This limitation is presumably due to the MPI operations as reported in other studies [14]. HPC applications often require numerous information exchange activities between computational nodes. In the case of the message-passing model, the nodes do not share the physical memory; therefore they communicate by sending information back and forth many times per second, which likely raises the possibility of latency when the number of node increases beyond a certain limit.

CONCLUSION
The present study demonstrated an application of TELEMAC-2D for flood inundation modeling. The 2D simulation results were consistent with observations from satellite imaging. The best agreement in flood extent between the modelled results and observations from MODIS satellite image analysis was obtained for the flood peak condition, while the agreement was worse for pre-and post-flood peak periods. The present study also showed that computing time can be reduced significantly with parallelization using multi-core processors. However, the speedup was modest when the number of cores increased beyond 4.