Flood Simulation by a Well-Balanced Finite Volume Method in Tapi River Basin, Thailand, 2017

Flood simulation of a region in southern +ailand during January 2017 is presented in this work. +e study area covers the Tapi river, the longest river in southern+ailand.+e simulation is performed by applying the two-dimensional shallow water model in the presence of strong source terms to the local bottom topography.+emodel is solved numerically by our finite volume method with well-balanced property and linear reconstruction technique. +is 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. +e initial conditions are prepared to match the simulated flow to the measurements recorded at the gauging stations. +e accuracy of the numerical simulations is demonstrated by comparing the simulated flood area to satellite images from the same period.+e results are in good agreement, indicating the suitability of the shallow water model and the presented numerical method for simulating floodplain inundation.


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. e finite volume method is an accurate numerical method that can be developed to solve the problem (for more details and reviews, see [1,2]). is 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][4][5]. e 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][6][7][8][9][10]11]. e 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. e 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]. e 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. e 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. e developed scheme is second-order accurate in space not only for smooth flow problem but also highgradient 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 ailand in January 2017. e study area is the Tapi river basin located in southern ailand. 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. e 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. e 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. e conclusions are finally given in Section 5.

Numerical Scheme for Shallow Water Equations
Flood over a large area can be simulated by considering the two-dimensional shallow water equations. e governing equations are given by where h is the water depth, u and v are the flow velocities in the x-and y-direction, respectively, g is the acceleration due to gravity, and z is the bottom elevation. S fx � −Cu are the friction terms in the in xand y-direction, respectively. Here, C � gn 2 /h 1/3 with n denoting Manning's Roughness coefficient. e conservative form of equations (1)- (3) is where U � (h, hu, hv) T , F � (hu, hu 2 + (1/2)gh 2 , huv) T , G � (hv, huv, hv 2 + (1/2)gh 2 ) T , and S � S b + S f � (0, −ghz x , −ghz y ) T + (0, ghS fx , ghS fy ) T . Next, we will apply our developed scheme to numerically approximate the water level and velocity at each location and time of our studied area.

Finite Volume Method.
A discretized form of (4) is where U ij is the approximation of U, over cell I ij � (x i−1/2 , x i+1/2 ) × (y j−1/2 , y j+1/2 ), given by with Δx � x i+1/2 − x i−1/2 and Δy � y j+1/2 − y j−1/2 for i � 1, 2, . . . , Kx and j � 1, 2, . . . , Ky. S ij is the source term approximation at cell I ij . F WAF and G WAF 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. e discretization in time is performed by the secondorder Runge-Kutta (RK2) method to ensure the secondorder accuracy in time of our method.

Weighted Average Flux (WAF).
We first consider the approximation of numerical flux in the x-direction at interface (x i+1/2 , y j ) denoted by F WAF i+1/2,j as follows: where U i+1/2,j is the solution of the Riemann problem from constant data U − i+1/2,j and U + i+1/2,j . e weighted average flux in two dimensions is proposed by [5,11]. It is composed of three flux components as follows: i+(1/2),j p , (p � 1, 2), where p denotes the component of the numerical flux vector at the interface and (F (k) i+(1/2),j ) p is the value of flux in the region k of the solution of the Riemann problem at component p. e first region (F (1) i+(1/2),j ) p � (F(U − i+(1/2),j )) p , the third region (F (3) i+(1/2),j ) p � (F(U + i+(1/2),j )) p , and the flux in the intermediate region (F (2) i+(1/2),j ) p are approximated by the Harten-Lax-van Leer (HLL) approach [11], where U − i+(1/2),j and U + i+(1/2),j are the solutions from the left and the right limits at the interface v − i+(1/2),j and v + i+(1/2),j are the velocities in y-direction from the left and the right limits at the interface. e weighted values ω k are calculated from the wave speeds in the left, the right, and the intermediate regions, respectively. WAF approximation in the y-direction, G WAF i,j+1/2 , 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].

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 by where σ ij is a that there are various forms. Here, we applied the minmod slope limiter given by where By the same concept, the linear reconstruction in the ydirection can be obtained by applying equations (9)-(12).

Well-Balanced Scheme.
A well-balanced concept is included in our developing scheme for preserving the stationary solution at the steady state. e 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 v � 0, 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 by ). e 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 e finite volume scheme is then expressed by where the numerical flux in the x-direction and the bottom slope terms are modified to Similarly, the numerical flux in the y-direction, G r i,j − 1/2 and G l i,j+1/2 , 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]).

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 Modelling and Simulation in Engineering in ailand, 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. e problems and simulation setting are given in the following sections.

Still Water Stationary State.
is experiment is performed to check the well-balanced property of the present scheme. e 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 e initial water depth is h + z � 16 m and velocity is zero in an entire domain. Simulation is run on uniform 1000 cells. e 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.
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 h + z � max(z, 6) m and initial velocity zero in the entire domain. e 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.
e 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.

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. e initial water depth is defined by Initial velocity is zero in the entire domain. Simulation is run on 1000 cells. e numerical results obtained by scheme II at time 0 s, 150 s, 400 s, and 1000 s are shown in Figure 3. e 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.

Partial Dam-Break.
is problem is considered on a 200 m × 200 m rectangular domain as shown in Figure 5. e initial water level is 10 m on the upstream side and 5 m on the downstream side. e bottom profile is assumed to be flat and frictionless. e partial dam-break is set at the middle of the domain. e simulation is performed using scheme II on 100 × 100 uniform grid cells. e final simulation time is 10 s. e 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.
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. e water level and its contour plot at the final time 7.2 s is shown in Figure 10. e plot of vector field at 7.2 s is shown in Figure 11. is shows the applicability of the present scheme for simulating wet/dry flow with moving shock on downstream.

Dam-Break Flows over ree Humps.
is 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. e problem is dam-break flow over three humps defined by e computational domain is rectangular with 75 m × 30 m. e dam is located at 16 m from the upstream boundary with initial water depth h + z � 1.875 m. Downstream is initially dry. We perform the simulations using scheme II with two cases of resolutions, 85 × 85, and 170 × 170 uniform grid cells with Manning's coefficient 0.018. e water depth profile, the contour plot, and the velocity fields at t � 12 s using 85 × 85 and 170 × 170 are shown in Figures 12-14, respectively. e 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. e water depth profile using 170 × 170 grid cells is slightly smoother than that using 85 × 85 grid cells.
is 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.
e comparison between scheme I and scheme II in terms of accuracy is shown in Figure 15. e contour plots at simulation time t � 12 s are slightly different. ese show that both schemes I and II can solve complex wet/dry ow interacting with bottom e ects 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.
Generally, the numerical schemes without conservative property may su er from mass-lost problem during wet/dry simulations. us, we have checked this issue by performing the next simulation using 85 × 85 grid cells with initial mass 893.3824 in a close system (all re ected boundaries). When the ow reaches steady state at 5000 s, the water mass remains 893.3824. e present scheme can preserve mass during time integration. Flow velocity is also zero. is result is shown in Figure 16.
It should be noted that the rectangular grid cells are used in this simulation and the water depth pro le is relatively tted to the circular dry bottom domain at steady state. More accurate solution can be obtained by applying more mesh re nements at the discretization step.

Case Study: Tapi Flood Simulation in
Thailand, 2017 In this section, we will apply the developed scheme to simulate the ooding in the Tapi basin ( Figure 17). is basin is located in the south of ailand. e 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 di erent as discussed in Section 3.4. us, only scheme I will be applied to simulate the real ow for saving CPU time.
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. e latitude is from 8.500400 ∘ N to 8.917500 ∘ N, and the longitude is from 99.084100 ∘ E to 99.417100 ∘ E. It covers an area of 1,620 km 2 .  Figure 18). e resolution corresponds to a grid cell size of 90 × 90 m. e cause of ood event is prolonged due to heavy rainfall within the basin. e amount of water entering the basin can be assessed from the measured discharge.    Table 1 [25]. e upstream ow is at station X37A, while the downstream ow is at station X217. We performed numerical simulation from January 9 to 13, 2017. Since the situation is heavy ood, the initial ood 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 out ow corresponding to the observed out ow given and assert an in ow corresponding to the observed in ow conditions on January 8. e in ow was set to that observed, and the simulation was run until the out ow reached the observed level. e obtained initial state was used as the initial condition for the simulation of the period from January 9 to 13. e contour plot of simulated water depth on January 8 is shown in Figure 19.

Modelling and Simulation in Engineering
From the initial conditions established for January 8, we run the first simulation to predict the flooded area on January 11. e 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. e result is displayed on the Google Earth along the Tapi river. e satellite image of the area on the same day can be obtained from [26] (Figure 20(a)). e 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. e amount of water is very high, exceeding the maximum capacity of the river channel. e 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.
Next, we continue the simulation to predict the flooded area on January 13. Heavy rainfall over the Tapi river basin continued during this period. is 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. e 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. e water depth level is relatively high which is approximately 0.5−1.5 m.
is 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 0.5−1.5 m when compared to the surrounding environment. ese 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.

Conclusions
Flood simulation during January 2017 in ailand is presented in this work. e study area is the Tapi river basin which covers many provinces in southern ailand. We apply the well-balanced finite volume method to solve the twodimensional shallow water model with strong source term from an irregular bottom profile in DEM format. e simulated period is from January 9 to 13. Discharge data are collected from two gauging stations. e 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. e simulation results show flooded area on January 11 and 13 that agree well with  the satellite images displayed by the Google Earth program. e 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 e 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.
e processed data are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest. also like to thank Mr. Brian Kubera at the Faculty of Science, Kasetsart University, for proof reading throughout the paper.