Abstract
Flood simulation of a region in southern Thailand during January 2017 is presented in this work. The study area covers the Tapi river, the longest river in southern Thailand. The simulation is performed by applying the two-dimensional shallow water model in the presence of strong source terms to the local bottom topography. The model is solved numerically by our finite volume method with well-balanced property and linear reconstruction technique. This technique is accurate and efficient at solving for complex flows in the wet/dry interface problem. Measurements of flows are collected from two gauging stations in the area. The initial conditions are prepared to match the simulated flow to the measurements recorded at the gauging stations. The accuracy of the numerical simulations is demonstrated by comparing the simulated flood area to satellite images from the same period. The results are in good agreement, indicating the suitability of the shallow water model and the presented numerical method for simulating floodplain inundation.
1. Introduction
To simulate flooding over an affected area of terrain, the two-dimensional shallow water model is one of the most efficient models. Since the nonlinear shallow water model is complicated, an efficient and accurate numerical method is required to find approximate solutions in terms of water depth and velocity. The finite volume method is an accurate numerical method that can be developed to solve the problem (for more details and reviews, see [1, 2]). This scheme requires an accurate numerical flux scheme for approximating the flux at cell interfaces in the shallow water equations. One extensively used scheme is Harten-Lax-van Leer (HLL) [3–5]. The modified version for two-dimensional problems is HLLC [5]. To obtain a more accurate approximation, the weighted average flux (WAF) has been introduced [4, 5–10, 11]. The WAF approximation is widely applied to the finite volume method. It can solve various types of problems (see, for instance, [5, 7, 9, 11, 12]).
In this work, we will apply the finite volume method with the WAF approximation for simulating a flood. The accuracy of numerical scheme depends on the method for approximating the bottom slope as discussed in [5]. Here, we improve the accuracy of numerical results by applying a linear reconstruction as described in [13]. The numerical scheme is second-order accurate in space for smooth flows with smooth bottom [13]. Instead of approximating bottom slope as presented in [11], a well-balanced scheme with bottom slope approximation is developed. The structured rectangular meshes are used due to its simple structure to develop a well-balanced WAF finite volume scheme [7, 10, 14, 15]. In addition, this kind of discretization can be applied directly to simulate real flood using the digital elevation model (DEM) from [16]. In application, we have to solve the nonlinear model that interacts with the nonlinear source term from the bottom topography. In this case, we improve the accuracy of numerical results by applying a linear reconstruction [1–3] for both water depth and bottom profile. The developed scheme is second-order accurate in space not only for smooth flow problem but also high-gradient water depth flow (see numerical experiments in [13]). To ensure second-order accuracy in time integration, we apply the second-order Runge–Kutta (RK2) method. When dealing with another source term of friction slope, the splitting implicit method proposed by Kesserwani and Liang [17, 18] is applied in our scheme. By combining all of these techniques, our numerical scheme is accurate and efficient; this will be demonstrated by the numerical experiments described below before applying the scheme to simulate a natural flood event recorded in Thailand in January 2017. The study area is the Tapi river basin located in southern Thailand. We propose a method to prepare the initial conditions necessary to conduct the flood simulation. All data utilized in our work for demonstrating the capability and performance of our numerical scheme are provided on the internet: the bottom topology, the satellite imagery, and the flow data collected at two gauging stations.
The paper is organized as follows. We describe the finite volume method with the weighted average flux and linear reconstruction for two-dimensional shallow water equations in Sections 2.1–2.3. The well-balanced finite volume method is presented in Section 2.4. Numerical tests are shown in Sections 3.1–3.3. Flood simulations are shown in Section 4. The conclusions are finally given in Section 5.
2. Numerical Scheme for Shallow Water Equations
Flood over a large area can be simulated by considering the two-dimensional shallow water equations. The governing equations are given bywhere is the water depth, and are the flow velocities in the x- and y-direction, respectively, is the acceleration due to gravity, and z is the bottom elevation. and are the friction terms in the in x- and y-direction, respectively. Here, with n denoting Manning’s Roughness coefficient.
The conservative form of equations (1)–(3) iswhere , , , and .
Next, we will apply our developed scheme to numerically approximate the water level and velocity at each location and time of our studied area.
2.1. Finite Volume Method
A discretized form of (4) iswhere is the approximation of , over cell , given bywith and for and . is the source term approximation at cell . and are the numerical fluxes in the x- and y-direction, respectively. We will apply the weighted average flux (WAF) to approximate these terms. Details are provided in Section 2.2.
The discretization in time is performed by the second-order Runge-Kutta (RK2) method to ensure the second-order accuracy in time of our method.
2.2. Weighted Average Flux (WAF)
We first consider the approximation of numerical flux in the x-direction at interface () denoted by as follows:where is the solution of the Riemann problem from constant data and .
The weighted average flux in two dimensions is proposed by [5, 11]. It is composed of three flux components as follows:where p denotes the component of the numerical flux vector at the interface and is the value of flux in the region k of the solution of the Riemann problem at component p. The first region , the third region , and the flux in the intermediate region are approximated by the Harten-Lax-van Leer (HLL) approach [11], where and are the solutions from the left and the right limits at the interface and are the velocities in y-direction from the left and the right limits at the interface. The weighted values are calculated from the wave speeds in the left, the right, and the intermediate regions, respectively. WAF approximation in the y-direction, and at other interfaces can be obtained similarly.
To avoid unexpected oscillations near a discontinuous water level profile, the WAF can be applied while enforcing the total variation diminishing (TVD); more details can be found in [13].
2.3. Linear Reconstruction
Second-order accuracy in space from constant data can be obtained by applying linear reconstruction [1–3, 13]. For instance, in the x-direction, the unknown variables are reconstructed before calculating numerical fluxes bywhere is a that there are various forms. Here, we applied the minmod slope limiter given bywhere
By the same concept, the linear reconstruction in the y-direction can be obtained by applying equations (9)–(12).
2.4. Well-Balanced Scheme
A well-balanced concept is included in our developing scheme for preserving the stationary solution at the steady state. The concept is obtained from considering just the one-dimensional flow at the steady state where the stationary solution must satisfy the following conditions:
Following Bermudez and Vazquez [19], a numerical scheme is called a well-balanced scheme if it satisfies the exact C-property, namely,
Similarly, for the two-dimensional flow problem, the conditions are
To obtain the well-balanced scheme, we follow the reconstruction approach proposed by Audusse et al. [20]. We reconstruct h at the interfaces in the x- and y-direction bywhere and .
The advantage of these reconstructions is that it can preserve the non-negativity of the water depth [20].
In this work, we propose a technique to modify the conservative variables to be
The finite volume scheme is then expressed bywhere the numerical flux in the x-direction and the bottom slope terms are modified to
Similarly, the numerical flux in the y-direction, and can be obtained by the decomposition of the bottom slopes in the third component. By applying these reconstructions, the finite volume method with the WAF approximation becomes a well-balanced scheme that preserves the exact C-property in two dimensions at the steady state. Some numerical tests are shown in the next section to confirm this property.
Moreover, to overcome the difficulties in calculating the source term for wet/dry problem, the friction term is approximated by the splitting implicit technique (see more details in [13, 17, 18]).
3. Numerical Tests
In this work, the scheme without linear reconstruction is referred to as scheme I, and the scheme with linear reconstruction is referred to as scheme II. Scheme II is second-order accurate in space for a smooth flow over a smooth bottom (see numerical experiment in [13]). Before applying our numerical method to simulate the observed flood event in Thailand, the accuracy of numerical scheme is checked by performing three test cases: still water stationary state, convergence of flow to still water stationary state, and partial dam-break flow. The problems and simulation setting are given in the following sections.
3.1. Still Water Stationary State
This experiment is performed to check the well-balanced property of the present scheme. The scheme is a well-balanced scheme if it satisfies the exact C-property; that is, the numerical solution should preserve the still water stationary solution at steady state. In this experiment, we consider a rectangular domain of 1500 m long with the discontinuous bottom defined by
The initial water depth is and velocity is zero in an entire domain. Simulation is run on uniform 1000 cells. The numerical result of water depth and velocity by scheme II at final time 100 s preserves still water stationary solution as shown in Figure 1.
(a)
(b)
To test the ability of the scheme to handle wet/dry still water stationary state, an additional experiment is performed on the same domain with initial conditions and initial velocity zero in the entire domain. The numerical result of water depth and velocity by scheme II with 1000 uniform cells at final time 100 s also preserves still water stationary solution as shown in Figure 2.
(a)
(b)
The numerical results from this experiment show that the present scheme satisfies the exact C-property for both wet and wet/dry problems over a discontinuous bottom; hence, the present scheme is a well-balanced scheme.
3.2. Convergence of Flow to Still Water Stationary State
To test the convergence of flow whether it reaches the still water stationary state, we perform a numerical experiment by considering a rectangular domain of 1500 m long with the discontinuous bottom defined in (20). Numerical schemes without special technique are usually unstable when dealing with wet/dry and discontinuous bottom. The initial water depth is defined by
Initial velocity is zero in the entire domain. Simulation is run on 1000 cells. The numerical results obtained by scheme II at time 0 s, 150 s, 400 s, and 1000 s are shown in Figure 3. The water depth and velocity at the final time 4000 s are shown in Figure 4. It shows that the numerical scheme can be used to simulate dam-break flow over a discontinuous bottom. It preserves well-balanced property without unexpected oscillations over the bottom.
(a)
(b)
(c)
(d)
(a)
(b)
3.3. Partial Dam-Break
This problem is considered on a 200 m × 200 m rectangular domain as shown in Figure 5. The initial water level is 10 m on the upstream side and 5 m on the downstream side. The bottom profile is assumed to be flat and frictionless. The partial dam-break is set at the middle of the domain.
The simulation is performed using scheme II on uniform grid cells. The final simulation time is 10 s. The water level and its contour plot at 7.2 s and the final time 10 s are shown in Figures 6 and 7, respectively. Plots of vector fields at 7.2 s and the final time 10 s are shown in Figures 8 and 9, respectively. Since there is no exact solution to this problem, we have checked the accuracy by comparing the water level profile at the same time with the results in [21, 22]. At the same simulation time, our results are very close to their results.
(a)
(b)
(a)
(b)
For the simulation of dry bed case, the initial water level is set to be 10 m on the upstream side and zero on the downstream side. The water level and its contour plot at the final time 7.2 s is shown in Figure 10. The plot of vector field at 7.2 s is shown in Figure 11. This shows the applicability of the present scheme for simulating wet/dry flow with moving shock on downstream.
(a)
(b)
3.4. Dam-Break Flows over Three Humps
This numerical experiment is considered to show the applicability of the present scheme for solving strong interaction between high-gradient water depth and friction bottoms in wet and dry case. The problem is dam-break flow over three humps defined by
The computational domain is rectangular with 75 m × 30 m. The dam is located at 16 m from the upstream boundary with initial water depth = 1.875 m. Downstream is initially dry. We perform the simulations using scheme II with two cases of resolutions, , and uniform grid cells with Manning’s coefficient 0.018.
The water depth profile, the contour plot, and the velocity fields at using and are shown in Figures 12–14, respectively. The obtained results agree well with the previous results presented by [18, 23]. As shown in Figures 12 and 13, the strong shock front that attacks the largest hump is detected correctly. The water depth profile using grid cells is slightly smoother than that using grid cells. This demonstrates that the developed scheme is capable of simulating flows that have both wet and dry beds, as well as the effect from a large bottom slope even using rectangular cells.
(a)
(b)
(a)
(b)
(a)
(b)
The comparison between scheme I and scheme II in terms of accuracy is shown in Figure 15. The contour plots at simulation time are slightly different. These show that both schemes I and II can solve complex wet/dry flow interacting with bottom effects that are usually represented in a real-world problem. When we consider the computational time, scheme I uses CPU time 39 s while scheme II uses 71 s performing on Intel(R) Core i7, CPU 3.6 GHz and RAM 8 GB. It is almost nearly double in this case due to the application of linear construction for every cell in each time step. Hence, for practical propose in the next section, we will apply just scheme I to simulate the real problem for saving computational time.
(a)
(b)
Generally, the numerical schemes without conservative property may suffer from mass-lost problem during wet/dry simulations. Thus, we have checked this issue by performing the next simulation using grid cells with initial mass 893.3824 in a close system (all reflected boundaries). When the flow reaches steady state at , the water mass remains 893.3824. The present scheme can preserve mass during time integration. Flow velocity is also zero. This result is shown in Figure 16.
(a)
(b)
It should be noted that the rectangular grid cells are used in this simulation and the water depth profile is relatively fitted to the circular dry bottom domain at steady state. More accurate solution can be obtained by applying more mesh refinements at the discretization step.
4. Case Study: Tapi Flood Simulation in Thailand, 2017
In this section, we will apply the developed scheme to simulate the flooding in the Tapi basin (Figure 17). This basin is located in the south of Thailand. The studied area covers an area of approximately 13,454 km2 with 8 tributaries. Most of the area is high and used for agriculture, especially fruit trees and rubber plantations. Since the bottom topography is not smooth, the numerical results obtained from schemes I and II are slightly different as discussed in Section 3.4. Thus, only scheme I will be applied to simulate the real flow for saving CPU time.
(a)
(b)
In our simulation, we consider the smaller area shown in Figure 17(b), which is located in Phrasaeng, Ban Na San, Wiang Sa, and Khian Sa districts. The latitude is from N to N, and the longitude is from E to E. It covers an area of 1,620 km2. This area includes a 96 km stretch of the Tapi river. There are two gauging stations, X37A in Phrasaeng district and X217 in Khian Sa district, separated by about 96 kilometers. The input topography obtained from the NASA Shuttle Radar Topographic Mission (SRTM) is in the digital elevation model (DEM) format (Figure 18). The resolution corresponds to a grid cell size of 90 × 90 m.
(a)
(b)
The cause of flood event is prolonged due to heavy rainfall within the basin. The amount of water entering the basin can be assessed from the measured discharge.
The observed discharge rates at both gauging stations on January 9–13, 2017, are shown in Table 1 [25]. The upstream flow is at station X37A, while the downstream flow is at station X217. We performed numerical simulation from January 9 to 13, 2017. Since the situation is heavy flood, the initial flood area and initial water depth over the Tapi river are unavailable. We propose a method to prepare the initial conditions of water depth and velocity to yield an outflow corresponding to the observed outflow given and assert an inflow corresponding to the observed inflow conditions on January 8. The inflow was set to that observed, and the simulation was run until the outflow reached the observed level. The obtained initial state was used as the initial condition for the simulation of the period from January 9 to 13. The contour plot of simulated water depth on January 8 is shown in Figure 19.
From the initial conditions established for January 8, we run the first simulation to predict the flooded area on January 11. The discharge observed at station X37A is set as an input or upstream condition throughout the simulation. We assume a constant discharge value throughout each day and update the value daily for each of the 3 days simulated. The result is displayed on the Google Earth along the Tapi river. The satellite image of the area on the same day can be obtained from [26] (Figure 20(a)). The difference between the simulation result and satellite image is shown in Figure 20(b). Inundation occurs along the Tapi river due to heavy rainfall that is reflected in the large values of discharge at both gauging stations. The amount of water is very high, exceeding the maximum capacity of the river channel. The simulation and satellite image are in good agreement, suggesting the correctness of our prepared initial conditions and the ability of our numerical scheme to solve complex flow problems.
(a)
(b)
Next, we continue the simulation to predict the flooded area on January 13. Heavy rainfall over the Tapi river basin continued during this period. This can be observed from the discharge value which reaches the maximum value 876.8 at station X217 on January 13. Severity of flooding is expected to increase as a direct result. The real satellite image data and the difference between simulation and real data are shown in Figures 21(a) and 21(b), respectively. As expected, the flooded area is larger. The water depth level is relatively high which is approximately . This location is in Chaiburi district, a populated area. Many people were affected by this severe flood event; the water height of the observation data in this area is also in the range when compared to the surrounding environment. These observations confirm the agreement between our simulations and the actual flood for both the spread of flooded area in two-dimensional plan and water height in the vertical direction.
(a)
(b)
5. Conclusions
Flood simulation during January 2017 in Thailand is presented in this work. The study area is the Tapi river basin which covers many provinces in southern Thailand. We apply the well-balanced finite volume method to solve the two-dimensional shallow water model with strong source term from an irregular bottom profile in DEM format. The simulated period is from January 9 to 13. Discharge data are collected from two gauging stations. The initial conditions are difficult to obtain, and here, we use the data from January 8 and then run numerical simulations until the numerical results are close to the observed data. The simulation results show flooded area on January 11 and 13 that agree well with the satellite images displayed by the Google Earth program. The range of predicted water depth at some locations is in the same range as that indicated by news photos. All of these simulation results show the capability and the performance of our numerical scheme to solve complex shallow water flows in real situations that can be applied to study other areas further.
Data Availability
The image and time series data supporting this manuscript are from previously reported studies and datasets, which have been cited by references [11, 13, 15, 25, 26] in our manuscript. The processed data are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to thank the Hydro and Agro Informatics Institute (HAII) for supporting data. They would also like to thank Mr. Brian Kubera at the Faculty of Science, Kasetsart University, for proof reading throughout the paper.