Hydrodynamic simulation of overland flooding over low-lying flat lands: A case study of the severe 2011 flood in Sam-Khok and Klong Luang districts, Thailand

The unprecedented flood event in the central Chao Phraya basin in 2011 was a catastrophic natural disaster affecting Thai society and economic sectors. During this event massive amounts of overland flow played a major role in damaging properties and environmental systems, with typical drainage systems rather ineffective in low-lying flat lands. The capability to effectively handle flood problem under these conditions is indeed a great challenge. To achieve such a goal, hydrological management schemes need to be improved based on insightful knowledge of flood hydrodynamics. In this study, we aim to gain insight into the aforementioned issue by performing and analyzing simulation of overland flow in low-lying flat lands. The flood plain of Sam-Khok and Klong Luang districts was selected as the study area and the simulation performed with an unstructured grid shallow water model together with a high-fidelity topographic data (LIDAR). The simulation results have been validated with field data (water trace height). The major hydrodynamic mechanisms such as flow patterns, flow magnitudes are characterized and the effective discharge quantified.


INTRODUCTION
After years of debate regarding the trigger for the 2011 flood in the Chao Phraya Watershed (CPYW), hydrological scientists and water resources experts have reached agreement that the most significant factor influencing the propagation of the flood over this region was an unusual synoptic condition characterized by the frequent Tropical Cyclone (TCs) events.Specifically, the TCs in early 2011 started from the onset of the monsoon season (mid-May) instead of the middle period of the season (July-August), (Thai Meteorological Department, 2011).As a result, wet season over the CPYW was longer in 2011 than in years with normal synoptic conditions.Hence, more rainy days were encountered with a heightened probability of high intensity storm events of days or weeks duration.These characteristics reflect both the amount and timing of rainfall which needed to be regulated by the watershed.In the case of the 2011 flood, the watershed was ineffectual because the total rainfall amount during the first few months of the TC season was already beyond the normal annual amount for the region.The total rainfall during the 2011 rainy season exceeded 140% of the average rainfall in the rainy season during the past two decades (Komori et al., 2012;Wongsa, 2013;Takebayashi et al., 2013).The response of the watershed is severely sensitive to triggering by TCs.As a result, triggering TCs were very significant because the input rainfall increased over a short drainage time, and thus, overflow could easily occur.Besides, influenced by severe weather with intense precipitation, the watershed lost its capacity to hold water in the soil; hence, excess rainfall flowed down slope further reducing the drainage capacity of the major river and reservoirs in the watershed.The other significant factor was the geomorphology of the watershed.The Chao Phraya delta consists mostly of gentle slope topography, and typically has poor drainage capabilities.Most of the 2011 water mass moved along the flood plain of the Chao Phraya River, exhibiting as a wide open channel flow.To precisely determine the aforementioned flow phenomena, the drainage system and primary flow behaviors should be well defined.In the present study, the term 'floodplain flow' is used to represent the hydrological events occurring across adjacent areas of the floodplain, rather than just in areas near the river.Here, the drainage system considered can be separated from the river channel based on natural or manmade banks.Hence, the source of flood waters flowing in the system considered originated from the levee overtopping and/or incoming flow through breached channel banks or dikes along the river.However, overland propagation with hydraulic discontinuities underlines the physics of the phenomena, hence, the term overland floods represents the flow characteristic in the case that the flood runs across the crop fields and urban areas.
The impact of overland flow was massive during the 2011 flood event; however, the capacity to effectively handle flood conditions is indeed a great challenge.To achieve such a goal, hydrological management schemes need to be improved based on insightful knowledge of flood hydrodynamics.In this study, we aim to gain insight into the aforementioned issue by simulation of overland flow over low-lying flat lands.The flood plain in Sam-Khok and Klong Luang districts was selected as the study area.The simulation was performed with an unstructured-grid shallow water model together with a high-fidelity topographic data (LIDAR).The simulation results have been validated with field data (water trace height).The major hydrodynamic mechanisms such as flow patterns and flow magnitudes are characterized and the effective discharge quantified.

ENVIRONMENTAL SETTING
The Chao Phraya flood plain consists of numerous sub drainage systems, catchment areas within low-lying land as well as gentle slope topography.These, however, are mostly non-recording gauge areas.This implies that there is a significant lack of supported hydrological data, which thus imposes even more challenges during flood event back analysis.Hence, an alternative approach for flow assumption validation and model calibration needs to be derived; this issue will be discussed later.Here, the flood plain in Sam-Khok and Klong Luang districts of Pathum Thani was selected as the study area (Figure 1).The size of this area is approximately 126 km 2 covering the irrigation units of Pathum Thani, Northern Rangsit and Ayutthya.The study area is located between latitudes of 14°2'31.88''N and 14°20'24.74''Nand longitudes of 100°30'41.12''Eand 100°37'5.17''E.It is mainly low-lying land and stretches in a gently-inclined plain to the Chao Phraya River mouth with a slope less than 0.1%.This drainage path is separated from the Chao Phraya River and the adjacent catchment system by the crest of the 3309 Road (Rd) and the Phahonyothin Rd as the western and eastern ridges, respectively.The elevation of both divides is approximately 3.50 meters above the mean sea level (m asl).
The drainage system can also be further divided into two sub catchment areas: the western sub-catchment area (east of the Chao Phraya River) covering most of the Bangsai-Sam-Khok district and the eastern sub-catchment area covering Chiang Rak Noi-Klong Nueng district.These sub-catchment areas are separated by a railway as a topographic divide in a north-south direction with an elevation of approximately 3.50 m asl (see Figure 1).Moreover, the two large polder systems have been constructed in the eastern sub-catchment area for flood protection of the Navanakorn Industrial Estate and Thailand Science Park (NSTDA)-Asian Institute of Technology (AIT)-Thammasat University (TU).In terms of major land-use, paddy fields cover almost the entire western sub-catchment area while the eastern sub-catchment area represents a more complicated pattern with mixed land-use consisting of paddy fields, urban areas and industrial estates.In terms of drainage system infrastructure, the two major canals (Prem Pracha Korn and Chiang Rak Noi canals) are located along the north-south and east-west axis of the system.Many water gates are also installed along the canal for regulating the water flow inside the system and from the adjacent sources.

Hydrological data survey
Using high precision measurement devices (i.e., Differential GPS), a field survey of the highest elevation of the 2011 flood event, as evidenced by the water trace and land surface level with respect to the mean sea level, can be performed in the study area with high confidence.The spatial distribution of the corresponding flood height data is shown in Figure 1.The southerly flow profile of the 2011 flood can be determined based on the trends in the survey data together with the available Light Detecting and Ranging Digital Elevation Models (LIDAR DEMs, http://thaisdi.gistda.or.th) as presented in Figure 2. Here, the flow profile can be identified as either M1 or M2 e.g.flow decelerates or accelerates, respectively, in the downstream direction.With the additional information that the flood height in the north was higher than flood heights in all other areas, we then conclude that the flow profile was a gradually varying flow over a mild slope channel with M2 type.Overall, the flow assumption for the flood event considered may be easily understood with only one-dimensional open channel approaches but reconstructing related hydrodynamics and processes together with hydraulic properties need to be performed with more comprehensive two or three dimensional models.Several localities influencing hydraulic head upstream and downstream of the drainage system must be determined properly.The proper hydrodynamic solver and associated boundary conditions need to be justified.

Hydrodynamics model
For capturing flood propagation mechanisms such as during the 2011 flood event in Thailand with adequate accuracy, and that reflects the hydraulic phenomena well, a proper hydrodynamic solver needs to be selected.In addition, at least two major capabilities need to be considered.First, the tool must be capable of fully resolving fluid motion conservation laws because the flow variable imbalance can significantly induce large discrepancy in flood depth and current speed.Second, a capability for handling high gradient solutions and jump phenomena is needed.This will allow one to perform the meshing process with ease because a specific meshing treatment such as the high density grid can be relaxed.Thus, a finite volume based solver fits the requirements perfectly.This is because the finite volume method can be applied to the conservation laws directly.Moreover, the specific nature of the finite volume approach, e.g.flux communication, allows us to handle problems with high gradient efficiently.Hence, the finite volume coastal ocean model, FVCOM (Chen et al., 2003) was used to investigate the present flood.The wet/dry method in this solver allows us to determine the region of effective flood damage accurately.
The governing equations used in the present study are Boussinesq, hydrostatic approximations of the primitive equations of mass and momentum conservation, which can be expressed in σ coordinates as follows: where u, v and ω are velocity components (m s -1 ) of the x (m), y (m) and σ directions.D = H + η is the total water depth (m) where H is the still water depth (m) measured from the mean sea level (m asl) and η is the fluctuation elevation of free surface (m), respectively.g is the acceleration due to gravity (m s -2 ); ρ 0 is the reference density (kg m -3 ); ρ is the perturbed density (kg m -3 ); K m is the vertical eddy viscosity (m 2 s -1 ); and F u and F v are the horizontal momentum diffusion terms.The kinematic boundary condition at surface is 0 The bottom momentum and kinematic boundary conditions are where (τ bx , τ by ) is bottom stress of x and y components.Specifically, the bottom stress is specified as The drag coefficient, C d , is determined by matching a logarithmic bottom layer to the model at a height, z b , above the bottom, thus where k = 0.4 is the von Karman constant and z 0 is the bottom roughness parameter.In this study, we set the bottom roughness to 1 cm.

Topographic data
Using the LIDAR DEMs together with canal and road geometry from the Royal Irrigation Department (RID) and high accuracy spot elevation, the model input can cover all essential details of the drainage system structure, including drainage divisions, relief of the flood plain, and existing infrastructure.

Computational domain and boundary condition
The study area was decomposed into 200,000 sub-domains represented by triangular grids.High density grids were located along the main road and flood prevention structures (Figure 3).The details of the boundary condition (BC) are summarized in Tables I and II and fully explained below.
For the baseline case or reconstruction of the 2011 flood event, the boundary condition representing the influence of any existing regulators (e.g., the topographic divides or dikes) must be imposed.Here, a wall condition was applied at both the western (3309 Rd with dike on the crest) and eastern boundaries (Phahonyothin Rd).An open boundary with fixed water elevation was also set at the southern side of the computational domain.For realistic simulation, two major polder systems and important infrastructure (such as railway and local roads) must be also included in the domain.To handle this, high density meshing was applied over those regions.The failure of the dike around the polder boundary was managed by partial modification of the polder boundary elevation.The inflow boundary condition was applied along the northern side of the domain.However, the flood discharge rate or other means of measuring current speed over the study area were unavailable.Hence, determining accurate incoming flood discharges to be used in the simulation was indeed a great challenge.The only clue that we have in this regard is the survey flow profile as mentioned earlier.Thus, to determine the inflow rate at the upstream boundary, we needed to find the discharge rate (Q) that provided the best fit between the model and survey flow profile.For this purpose, an experimental Q was selected in the range of 500 to 1000 cubic meter per second (m 3 s -1 ), which is approximately 15% and 30% of the carrying capacity of the Chao Phraya River.
For scenario-based studies, various boundary conditions were applied on the western boundary under two different schemes.
The S1 scheme was based on assumption that the capacity of Chao Phraya River had been exceeded.Hence, conventional practices such dike and polder systems may be applied along the western boundary.In order to find a solution for flood management guidelines, we also assumed that the drainage capacity of the downstream area was enhanced, which is the only reasonable assumption.Thus, the downstream boundary is now set to be radiative while all other boundaries are set the same as in the baseline case.
For the S2 scheme, the boundary condition was similar to that of the S1 scheme except for the western boundary where radiative boundary conditions were imposed.The conditions were set to transfer fluxes out of the computational domain.However, this is quite an extreme assumption meaning that we allowed flood flow to adjacent areas.For both scenarios, the two polder systems were included in the domain by applying wall boundary conditions along their boundaries.Furthermore, a new dike of 4.5 m asl along the Chiang Rak Noi canal with a narrow gap was also applied by modification of the elevation data.The incoming flood rate for the inflow boundary was identical or twice that of the effective discharge determined from the baseline case.
For the sake of clarity, we have also included the simulation schematic for both the baseline and scenario cases in Figure 4.

Model validation and analysis for the baseline case
Using the FVCOM model with sufficient topographical data together with proper flow assumptions, the propagation of the flood during the peak flow period and its corresponding flow profile over the study area can be reconstructed (Figure 2).The flow profile resolved by the model was in good agreement with that of the field survey with small errors.This indicates that the simulation can capture the behaviors of the 2011 flood event under proper hydrodynamic conditions and hydraulic properties.Hence, we finally estimate the appropriate value of the incoming flood discharge to be within 940 m 3 s -1 .Furthermore, based on the present simulation results, the underlying mechanism of the flood event can now be more accurately characterized as follows.
As presented in Figure 5, the southward direction characterizes the propagation pattern of the flood, and also suggests the source of incoming flow.It is known that most of the inundated water during the peak flow was slowly drained From Figure 6, we can see that the drainage divides and infrastructure inside the drainage system impede the flow movement.The aforementioned mechanism can be seen from the motion of flood along the western part.As the hydraulic head along the western drainage divide was preserved by constructing a dike on the 3309 Rd, the flow direction conformed to the shape of the drainage divide.This suggests that the flow direction could turn westward as head loss may occur along the 3309 Rd.Theoretically, this would enhance the drainage capability and reduce flood levels over the western part of the study area.This assumption will be applied in another simulation in order to gain more understanding of flood mechanism sand assess how flood levels and drainage capability would be different under different flood prevention practices and incoming flow rates.The computational set up and results are presented next.

Scenario-based assessment
Probable flood magnitude and associated hydrodynamic processes are presented and discussed in the previous section.Also discussed in the background section is the reason behind the 2011 flood, which relates to the triggering of extraordinary synoptic conditions.Considering the related issue of such climate dynamics over the tropical region requires a high precision prediction which is the subject of ongoing research (Veerakachen et al., 2014).Though, in the present study, we have not yet obtained any reliable information as to whether or not the reoccurrence of the severe climatic conditions (such as triggering by stronger and more frequent TCs) would take place again in the near future.This encourages many related agencies and communities to pay more attention to flood mechanisms for a specific drainage system (Kwak et al., 2012) as well as the whole system like CPYW.This is mainly for improvement of preparedness and prediction skills.In the current study, the behavior of the 2011 flood was further determined under different prevention/protection schemes as well as various hypothetical incoming flow magnitudes.The setup of the computational domain and the details of each scenarios were already described in the methodology, with the results of each scheme as follows.
S1 scheme: In this scheme, all components of the catchment were almost similar to that of the baseline case, except that permanent dikes (sheet piles) were installed along an east-west direction along the Chiang Rak Noi canal leaving a narrow gap in the middle.We also assumed that the drainage capacity of the downstream area was enhanced.No leakage or no failure of the dike bounding the polder systems was allowed.We assumed that the incoming flood discharge of around 940 m 3 s -1 or more has returned with the inundation map for the scenario shown in Figure 7 (left panel).From Figure 7, the inundation area can be seen to be reduced in the middle to southern side of the study area.This would mostly be due to the effect of the permanent dike preventing the area beneath it from severe flooding and creating a detention area above the dike.Hence, the flood level over the latter area was much less than that of the baseline case.However, as the magnitude of the incoming flood discharge increases greatly, this practice is would be rather ineffective as the inundation area beneath the dike increases.
S2 scheme: For this scheme, the model configuration was exactly the same as in the previous case.However, to resolve the increasing drainage rate in the study area, radiative boundary conditions were applied to the western side of the domain.The condition of the incoming flood discharge was the same as the S1 scheme.The results of the scheme are presented in Figure 7 (right panel).Overall, in terms of inundation, the spreading of flood area given by setup of S2 is relatively lower than that of S1.However, the flood level was still high compared to the baseline case.One of the interesting points is that there was a positive head gradient through the western boundary (as the dike on the top of the 3309 Rd had been removed) even though the incoming flood discharge was rather large in magnitude.This tends to be the most effective practice in terms of enhancing the drainage rate of the western part of the drainage area.This, however, is subject to the carrying capacity of Chao Phraya River.Theoretically, the western part can alternate between a retention area or drainage area depending on the aforementioned condition.
Other practices such as constructing a dike along the Chiang Rak Noi canal leaving a narrow gap may yield less efficient flood reduction in downstream areas and might result in even higher water levels than that of the 2011 flood in the upper area of the dike.

CONCLUSIONS
The development of numerical simulation models for flood management and field surveys of highest water surface levels during the Thailand 2011 flood has been achieved.The flood plain in Sam-Khok and Khong Luang districts including the Navanakorn Industrial Estate, Thailand Science Park, Asian Institute of Technology and Thammasat University was selected as the study area.Precise surveys of land surface level and the highest elevation of the 2011 flood, with respect to the mean sea level, were performed.Furthermore, a numerical hydrodynamic model was developed to simulate the behaviors of a massive flood over a low slope flood plain and applied to investigate flood management practice in the study area and its vicinity.Using the hydrodynamic model FVCOM together with high accuracy elevation data (LIDAR DEMs) and flow assumptions guided by survey data, the flow characteristics of the 2011 massive flood can be correctly modeled.
The numerical results were in good agreement with the observed data.The incoming flow rate from upstream was estimated to be 940 m 3 s -1 .The numerical model was employed to investigate the effects of various flood prevention measures.
It was found that increasing the water draining through the western boundary (no dike on the top of the 3309 Rd) would be the most effective scheme in reducing inundation in the area (subject to the carrying capacity of Chao Phraya River).

Figure 1 .
Figure 1.Overview of the study area and details of its parts (map drawn according to UTM (Zone 47)).The figure also includes LIDAR DEM and overlaid spot elevation of 2011 flood level traces (m asl).The red polygon encloses the Navanakorn Industrial Estate and the orange polygon encloses Thailand Science Park, Asian Institute of Technology and Thammasat University

Figure 3 .
Figure 3. Computational domains for the baseline case (upper) and scenario-based study (lower)

Figure 4 .
Figure 4. Simulation schematic for baseline case (a) and scenario cases (b)

Figure 7 .
Figure 7. Maximum flood level in meters given from the S1 (left) and S2 (right) scenarios.Upper belongs to the 940 m 3 s -1 inflow while lower is driven by the inflow magnitude of 1880 m 3 s -1