Research Article  Open Access
Sutatip Vichiantong, Thida Pongsanguansin, Montri Maleewong, "Flood Simulation by a WellBalanced Finite Volume Method in Tapi River Basin, Thailand, 2017", Modelling and Simulation in Engineering, vol. 2019, Article ID 7053131, 13 pages, 2019. https://doi.org/10.1155/2019/7053131
Flood Simulation by a WellBalanced Finite Volume Method in Tapi River Basin, Thailand, 2017
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 twodimensional 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 wellbalanced 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 twodimensional 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 HartenLaxvan Leer (HLL) [3–5]. The modified version for twodimensional 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 secondorder accurate in space for smooth flows with smooth bottom [13]. Instead of approximating bottom slope as presented in [11], a wellbalanced scheme with bottom slope approximation is developed. The structured rectangular meshes are used due to its simple structure to develop a wellbalanced 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 secondorder accurate in space not only for smooth flow problem but also highgradient water depth flow (see numerical experiments in [13]). To ensure secondorder accuracy in time integration, we apply the secondorder 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 twodimensional shallow water equations in Sections 2.1–2.3. The wellbalanced 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 twodimensional shallow water equations. The governing equations are given bywhere is the water depth, and are the flow velocities in the x and ydirection, respectively, is the acceleration due to gravity, and z is the bottom elevation. and are the friction terms in the in x and ydirection, 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 ydirection, 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 secondorder RungeKutta (RK2) method to ensure the secondorder accuracy in time of our method.
2.2. Weighted Average Flux (WAF)
We first consider the approximation of numerical flux in the xdirection 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 HartenLaxvan Leer (HLL) approach [11], where and are the solutions from the left and the right limits at the interface and are the velocities in ydirection 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 ydirection, 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
Secondorder accuracy in space from constant data can be obtained by applying linear reconstruction [1–3, 13]. For instance, in the xdirection, 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 ydirection can be obtained by applying equations (9)–(12).
2.4. WellBalanced Scheme
A wellbalanced concept is included in our developing scheme for preserving the stationary solution at the steady state. The concept is obtained from considering just the onedimensional 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 wellbalanced scheme if it satisfies the exact Cproperty, namely,
Similarly, for the twodimensional flow problem, the conditions are
To obtain the wellbalanced scheme, we follow the reconstruction approach proposed by Audusse et al. [20]. We reconstruct h at the interfaces in the x and ydirection bywhere and .
The advantage of these reconstructions is that it can preserve the nonnegativity 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 xdirection and the bottom slope terms are modified to
Similarly, the numerical flux in the ydirection, 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 wellbalanced scheme that preserves the exact Cproperty 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 secondorder 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 dambreak 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 wellbalanced property of the present scheme. The scheme is a wellbalanced scheme if it satisfies the exact Cproperty; 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 Cproperty for both wet and wet/dry problems over a discontinuous bottom; hence, the present scheme is a wellbalanced 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 dambreak flow over a discontinuous bottom. It preserves wellbalanced property without unexpected oscillations over the bottom.
(a)
(b)
(c)
(d)
(a)
(b)
3.3. Partial DamBreak
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 dambreak 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. DamBreak Flows over Three Humps
This numerical experiment is considered to show the applicability of the present scheme for solving strong interaction between highgradient water depth and friction bottoms in wet and dry case. The problem is dambreak 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 realworld 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 masslost 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 km^{2} 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 km^{2}. 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 twodimensional 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 wellbalanced finite volume method to solve the twodimensional 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.
References
 R. J. Leveque, Finite Volume Method for Hyperbolic Problem, Cambridge University Press, Cambridge, UK, 2005.
 R. J. Leveque, Numerical Methods for Conservation Laws, Birkhäuser Verlag Press, Basel, Switzerland, 1992.
 E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, Berlin, Germany, 1999.
 E. F. Toro, ShockCapturing Methods for FreeSurface Shallow Flows, John Willy and Sons Ltd., Hoboken, NJ, USA, 2001.
 R. Ata, S. Pavan, S. Khelladi, and E. F. Toro, “A Weighted Average Flux (WAF) scheme applied to shallow water equations for reallife applications,” Advances in Water Resources, vol. 62, pp. 155–172, 2013. View at: Publisher Site  Google Scholar
 E. F. Toro, “A weighted average flux method for hyperbolic conservation laws,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 423, no. 1865, pp. 401–418, 1989. View at: Publisher Site  Google Scholar
 E. F. Toro, “Riemann problems and the WAF method for solving the twodimensional shallow water equations,” Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences, vol. 338, no. 1649, pp. 43–68, 1992. View at: Publisher Site  Google Scholar
 E. F. Toro and P. L. Roe, “A hybridised highorder random choice method for quasilinear hyperbolic systems,” in Proceedings of 16th International Symposium on Shock Tubes and Waves, H. Gronig, Ed., p. 701, Aachen, Germany, July 1987. View at: Google Scholar
 M. Maleewong, “Modified predictorcorrector WAF method for the shallow water equations with source terms,” Mathematical Problems in Engineering, vol. 2011, Article ID 178491, 17 pages, 2011. View at: Publisher Site  Google Scholar
 S. J. Billett and E. F. Toro, “On WAFtype schemes for multidimensional hyperbolic conservation laws,” Journal of Computational Physics, vol. 130, no. 1, pp. 1–24, 1997. View at: Publisher Site  Google Scholar
 Y. Loukili and A. W. Soulaimani, “Numerical tracking of shallow water waves by the unstructured finite volume WAF approximation,” International Journal for Computational Methods in Engineering Science and Mechanics, vol. 8, no. 2, pp. 75–88, 2007. View at: Publisher Site  Google Scholar
 E. D. FernándezNieto and G. NarbonaReina, “Extension of WAF type methods to nonhomogeneous shallow water equations with pollutant,” Journal of Scientific Computing, vol. 36, no. 2, pp. 193–127, 2008. View at: Publisher Site  Google Scholar
 T. Pongsanguansin, M. Maleewong, and K. Mekchay, “Shallowwater simulations by a wellbalanced WAF finite volume method: a case study to the great flood in 2011, Thailand,” Computational Geosciences, vol. 20, no. 6, pp. 1269–1285, 2016. View at: Publisher Site  Google Scholar
 R. J. Leveque, “Balancing source terms and flux gradients in highresolution Godunov methods: the quasisteady wavepropagation algorithm,” Journal of Computational Physics, vol. 146, no. 1, pp. 346–365, 1998. View at: Publisher Site  Google Scholar
 S. M. Amiri, N. Talebbeydokhti, and A. Baghlani, “A twodimensional wellbalanced numerical model for shallow water equations,” Scientia Iranica, vol. 20, no. 1, pp. 97–107, 2013. View at: Publisher Site  Google Scholar
 http://srtm.csi.cgiar.org.
 G. Kesserwani and Q. Liang, “Wellbalanced RKDG2 solutions to the shallow water equations over irregular domains with wetting and drying,” Computers & Fluids, vol. 39, no. 10, pp. 2040–2050, 2010. View at: Publisher Site  Google Scholar
 G. Kesserwani and Q. Liang, “A discontinuous Galerkin algorithm for the twodimensional shallow water equations,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 49–52, pp. 3356–3368, 2010. View at: Publisher Site  Google Scholar
 A. Bermudez and M. E. Vazquez, “Upwind methods for hyperbolic conservation laws with source terms,” Computers & Fluids, vol. 23, no. 8, pp. 1049–1071, 1994. View at: Publisher Site  Google Scholar
 E. Audusse, F. Bouchut, M. O. Bristeau, R. Klein, and B. Perthame, “A fast and stable wellbalanced scheme with hydrostatic reconstruction for shallow water flows,” SIAM Journal on Scientific Computing, vol. 25, no. 6, pp. 2050–2065, 2004. View at: Publisher Site  Google Scholar
 H.R. Vosoughifar, A. Dolatshah, and S.K. S. Shokouhi, “Discretization of multidimensional mathematical equations of dam break phenomena using a novel approach of finite volume method,” Journal of Applied Mathematics, vol. 2013, Article ID 642485, 12 pages, 2013. View at: Publisher Site  Google Scholar
 K. Anastasiou and C. T. Chan, “Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshs,” International Journal for Numerical Methods in Fluids, vol. 24, no. 11, pp. 1225–1245, 1997. View at: Publisher Site  Google Scholar
 Q. Liang and A. G. L. Borthwick, “Adaptive quadtree simulation of shallow flows with wetdry fronts over complex topography,” Computers & Fluids, vol. 38, no. 2, pp. 221–234, 2009. View at: Publisher Site  Google Scholar
 http://www.rid.go.th/lproject/const/water25/25river/22tapee.
 http://wwwhydro8.com.
 http://www.gistda.or.th.
Copyright
Copyright © 2019 Sutatip Vichiantong et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.