- Email: [email protected]
Analysis of an open source quadtree grid shallow water flow solver for flood simulation
Quaternary International 384 (2015) 118e128
Contents lists available at ScienceDirect
Quaternary International journal homepage: www.elsevier.com/locate/quaint
Analysis of an open source quadtree grid shallow water flow solver for flood simulation Hyunuk An a, Soonyoung Yu b, *, Giha Lee c, Yeonsu Kim d a
Department of Agricultural and Rural Engineering, Chungnam National University, Daejeon 305-764, Republic of Korea Research Institute for Social Criticality, Pusan National University, Busan 609-735, Republic of Korea c Department of Constructional Disaster Prevention Engineering, Kyungpook National University, Sangju 742-711, Republic of Korea d International Water Resources Research Institute, Chungnam National University, Daejeon 305-764, Republic of Korea b
a r t i c l e i n f o
a b s t r a c t
Article history: Available online 19 February 2015
The shallow water flow solver within Gerris Flow Solver implements second-order accurate Gudunov type numerical schemes, with preserving the balance of source and flux terms on quadtree cut cell grids. The solver provides flexible grid generation, with using adaptive quadtree grids and the cut cell method, which can improve computational efficiency in 2D simulations. The solver, however, cannot implicitly discretize nonlinear friction without linearization since the current version of Gerris supports the implicit discretization only for linear friction. We compare implicit treatment for Manning's friction to the linearization. Simulation results using a benchmark test show the difference between two methods is only significant when the roughness coefficient is unrealistically high. Two real flood events, Malpasset dam break in France and Baeksan levee failure in Korea, are simulated using the quadtree grid based shallow water model, with adaptively refining meshes near water fronts and river boundaries. Comparison of simulation results with observations and measurements demonstrates that the grid adaptation can save approximately 85%e95% of the computational cost while preserving the accuracy. The research result shows that the linearized friction using the current version of Gerris can be applied for flood simulation and multiple inflow boundaries can be treated as source terms. © 2015 Elsevier Ltd and INQUA. All rights reserved.
Keywords: Gerris Flood simulation Quadtree grid Bed friction Inflow boundary
1. Introduction Numerical modeling of shallow water flow is an important research topic for flood risk management. The simulation results play a significant role in national decision-making on flood prevention and control. One-dimensional shallow water models are quite often preferred in the field of engineering (e.g., FLUCOMP; Fread, 1993; Ervine and MacLeod, 1999; HEC-RAS; ISIS; MIKE11), mostly due to the fast computation (Syme, 2011). The rapid development in computer technologies, however, makes twodimensional shallow water flow models acceptable in practice. A number of 2D models have been developed (i.e., Bates and De Roo, 2000; Yoon and Kang, 2004; Mignot et al., 2006; George, 2011; Kim and Cho, 2011). In general, the 2D models are more accurate and reliable than 1D models for complex flow simulation given that the
* Corresponding author. E-mail addresses: [email protected] (H. An), [email protected] (S. Yu), [email protected] (G. Lee), [email protected] (Y. Kim). http://dx.doi.org/10.1016/j.quaint.2015.01.032 1040-6182/© 2015 Elsevier Ltd and INQUA. All rights reserved.
water is allowed to move only in the longitudinal direction in 1D models. 2D modeling is, however, still computationally heavy to engineers despite the advances in computing technologies, particularly for flood risk assessment or real time flood forecasting, which forces engineers to choose 1D models over 2D, giving up accuracy. One of the solutions to improve computational efficiency in 2D modeling is the implementation of adaptive mesh refinement (AMR), in which grid resolution is adjusted with following flow features. Generally, the finer grid resolution is achieved along wave fronts or complex geometry than within inundated area. The optimal grid resolution is time varying given that flows change with time. The method of adaptive meshing makes numerical simulations much more efficient (Berger and Oliger, 1984; Berger and Colella, 1989; Borthwick et al., 2001), and its adoption for numerical simulations becomes practice in engineering and science, including shallow water flow modeling. George and LeVeque (2008) used block-structured AMR for tsunami modelling. Their AMR algorithms allow nesting of rectangular subgrids of multiple levels and the subgrids evolve spatially and temporally. George (2011) applied the same model to simulate dam break flow and showed
H. An et al. / Quaternary International 384 (2015) 118e128
that the shallow water model with AMR is efficient and suitable for flood simulation. Alternatively, Liang et al. (2007) and Liang (2010) simulated inundation scenarios on adaptive quadtree grids. The adaptive quadtree grid is a way to carry out AMR and uses a quadtree and hierarchical data structure. Liang (2010) compared the performance of their model with the commercial code and showed that the shallow water model on adaptive quadtree grids is highly efficient while preserving accuracy. Popinet (2011, 2012) developed a quadtree adaptive tsunami model to solve shallow water equations within the Gerris Flow Solver framework (Popinet, 2003). Lee et al. (2011) presented a dynamically adaptive quadtree grid generation system for the solution of a two-dimensional two-layer shallow water model in order to improve computational efficiency. Gerris Flow Solver has the capability of adaptive quadtree grid generation. The Gerris solver was initially designed to solve the incompressible Euler equation, and then has been extended to multiphase incompressible NaviereStokes equations (Popinet, 2009), spectral wave models (Popinet et al., 2010), and shallow water equations to simulate tsunami propagation (Popinet, 2011, 2012). Recently, An and Yu (2012) applied the cut cell method to the shallow water flow solver of Gerris, which provides the model with more flexibility. In addition, the Gerris Flow Solver is open source software, and users can see the internal structure and even modify the codes for better representation. Moreover, the Gerris Flow Solver is well equipped with auxiliary programs such as GfsView for visualization, which most open source software are lack of as Chien (2008) mentioned. Gerris supports a variety of file formats such as Cartesian Grid Data (CGD), ESRI ASCII grid, xyz files to treat initial and boundary conditions, which facilitates specifying a time varying boundary condition. Gerris can produce outputs as text, image, MPG video, or KMZ files for Google Earth as well as GfsView files. Lee et al. (2013) compared the simple pre- and postprocessing in Gerris to those in FLUMEN (FLUvial Modeling Engine), and verified the computational efficiency of quadtree grids over triangular grids. The objective of this paper is to analyze the performance of the shallow water flow solver within Gerris Flow Solver for flood simulation in 2D through real flood cases. The motivation of this paper is the fact that there are few researches with applying the quadtree grid shallow water model to real flood events. Liang et al. (2007) and Liang (2010) simulated flood scenarios using a quadtree model for flood risk assessment. Their simulation results are in good agreement with those obtained by using a commercial model, called TUFLOW. In addition, Kesserwani and Liang (2012) simulated the Malpasset dam break event using the dynamically adaptive grid based shallow water model. The shallow water flow solver of Gerris has been intensively verified through real tsunami cases (Popinet, 2011, 2012) as well as a number of benchmark tests (Popinet et al., 2010; An and Yu, 2012). However, it is still informative to analyze the model performance on real flood events for potential Gerris users because Gerris has several favorable features as stated above. In addition, Gerris has several options to treat inflow boundaries as well as friction, and allows choosing mass conservation or lake-at-rest equilibrium when the dynamic adaptive mesh refinement is applied for flood simulations. The effects of those options are investigated through a benchmark test and real flood simulations in this paper.
119
broad range of flood simulations in practice. A hyperbolic conservation form of the equations is written as follows:
vq vf vg þ þ ¼ s; vt vx vy
(1)
where t denotes time, x and y are Cartesian coordinates, and q, f, g, and s are vectors representing conserved variables, fluxes in the xand y-direction, and source terms, respectively. By neglecting the friction term, the vectors can be written as:
0
h
0
1
B C q ¼ @ hu A; 0 B g¼@
hv hv
1 hu . 2 B 2 C f ¼ @ hu þ gh 2 A; 1
huv . C A; hu þ gh2 2 2
huv 0
1 0 B C s ¼ @ hgzx A;
(2)
hgzy
where h is water depth, u and v are depth-averaged velocity components in the x and y-direction, g is the acceleration of gravity, and zx(¼ vz/vx) and zy(¼ vz/vy) are bed slopes in the x and y-direction. The friction term will be added to s in Eq. (2) for completeness, without changing the fundamental structure of the equations in Section 2.3. 2.2. Spatial and temporal discretization Recently, a number of well-balanced numerical schemes have been proposed for shallow water flow equations (e.g., Audusse et al., 2004; Audusse and Bristeanu, 2005; Caleffi et al., 2006; Kim et al., 2008; Caleffi and Valiani, 2009) given that the gradient of water depths and the bed slope terms are often imbalanced in a discretized domain in case of irregular beds. The imbalance may cause unphysical perturbations when dryewet transitions occur or flows reach a steady state condition, which significantly affects the accuracy of flood prediction. In order to solve the problem of imbalance, Gerris implements the hydrostatic reconstruction technique proposed by Audusse et al. (2004), whose satisfactory performance has been verified in Popinet (2011) and An and Yu (2012). Following An and Yu (2012), the second-order accurate discretization on a two dimensional quadtree cut cell grid is written as:
Ai;j nþ1 qi;j qni;j þ F i;j ¼ S i;j þ Sci;j Dt
(3)
where Ai,j is the area of the cell (i, j), Sci,j is an additional source term to balance in the presence of a source, and Fi,j is the numerical flux. A MUSCL-type unsplit discretisation is combined with a predictorcorrector time-stepping scheme to construct a second-order accurate Godunov type scheme (Popinet, 2011). The predicator step computes a temporary cell-center value over the half time step Dt/ 2: nþ1=2
Ai;j qi;j
¼ Ai;j qni;j
Dt n F i;j S ni;j Sc ni;j ; 2
(4)
The corrector step updates the conservative values over a time step Dt as: 2. Numerical model 2.1. Governing equation Two-dimensional depth-integrated shallow water equations are used for flood event simulation. They have performed well for a
nþ1=2 nþ1=2 nþ1=2 : Ai;j qnþ1 ¼ Ai;j qni;j Dt F i;j S i;j Sc i;j i;j
(5)
The above scheme is second-order accurate in space and time, and the details of discretization and the well-balanced scheme are described in An and Yu (2012).
120
H. An et al. / Quaternary International 384 (2015) 118e128
where Dq ¼ qnþ1 qn. Ignoring the high-order terms and substituting Eq. (10) into Eq. (9) give
2.3. Friction term discretization Frictional effects should be incorporated to the right-hand side of Eq. (1) for real flood simulation. The Manning formula is considered as:
DðDtÞ ¼
*
(12)
0 0 B vSf2 B vSf B ¼ B vh B vq @ vSf 3
vh
0 vSf2 vðhuÞ vSf3 vðhuÞ
0 1 vSf2 C C vðhvÞ C C: C vSf A
(13)
3
vðhvÞ
n
Denoting Sf ðDtÞ¼ D1 ðDtÞSf n yields the implicit scheme written as:
¼ Ai;j qni;j
Dt n F i;j S ni;j Sc ni;j Ai;j Sf ni;j ; 2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 n : ¼ Ai;j qi;j Dt F i;j S i;j Sc i;j Ai;j Sf i;j
nþ1=2
Ai;j qnþ1 i;j
n
where I is a unit matrix and
cretization methods for general source terms while implicit discretization for linear friction and Coriolis force. The explicit method is written as follows:
Ai;j qi;j
vSf vq
Dt n Dt n Dt nþ1=2 nþ1=2* F i;j S ni;j Scni;j ; qi;j ; ¼ qi;j þ Sf 2 2 2 * nþ1=2 nþ1=2 nþ1=2 nþ1=2 ; qnþ1 ¼ Ai;j qni;j Dt F i;j S i;j Sci;j ¼ qnþ1 þ DtSf ðDtÞ: i;j i;j
nþ1=2*
Ai;j qnþ1 i;j
I Dt
(6)
where n is the Manning's roughness coefficient. The addition of this friction term may cause the numerical instability in a discretized domain, in particular when a dryewet transition exists. In other words, if water depth becomes close to zero, the Manning's friction can result in an exaggerated force, which reverses the flow in a nonphysical manner. Therefore, the friction term should be properly treated. The current version of Gerris provides explicit dis-
Ai;j qi;j
(11)
where
0
1 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 1 B gn2 u u2 þ v2 C B C Sf1 B C 1=3 C; Sf ¼ @ Sf2 A ¼ B h B C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C B Sf3 @ 2 2 2 A gn v u þ v h1=3
qnþ1 ¼ qn þDtD1 ðDtÞSf n ;
¼ Ai;j qni;j
(7) However, the above equations may cause the numerical instability although it can be straightforwardly applied. Therefore, implicit treatment is applied for numerical stability in this paper. Three implicit methods are implemented and compared because only the linear friction is treated implicitly in the current version of Gerris. The first method is a recently developed implicit discretization based on splitting method (Liang and Marche, 2009). The implicit discretization treats the friction term with using the splitting method, which is equivalent to solve the following equation:
dq ¼ Sf : dt
Liang and Marche (2009) presented this method in onedimensional spaces, and we fully extend their method in twodimension in Eqs. (11)e(14). The friction force in Eq. (14) can cause an exaggerated reverse flow, which makes no physical sense, although the friction terms are discretized implicitly. Therefore, the following condition is added:
n
2qn n ; if Sf > Dt nþ1=2 qnþ1=2 qn ; if Sf ¼ > Dt Dt
Sf ¼ qn Sf
nþ1=2
2 Dt
(15)
where Sf and q are momentum components of Sf and q, respectively. With the above condition, the maximum effect of friction forces is to stop the fluid. The second method is another implicit discretization presented in Liang (2010), who applied the method of Liang and Marche (2009) on two-dimensional equations in a directional manner, which gives
(8)
The Euler implicit discretization of the above equation yields
qnþ1 qn ¼ Sf nþ1 : Dt
(9)
The right-hand side of Eq. (9) can be represented by the Taylor series expansion as:
Sf nþ1 ¼ Sf n þ
(14)
vSf vq
n
ðDqÞ þ o Dq2 ;
(10)
00 B B0 vSf B ¼B B vq @ 0
0 vSf2 vðhuÞ 0
0
1
C C C C: C vSf3 A vðhvÞ 0
(16)
Note that Eq. (16) ignores the effect of off-diagonal entries as compared with Eq. (13). The third method is the linearization of friction terms. The current version of Gerris allows the implicit discretization of the linear friction as follows:
H. An et al. / Quaternary International 384 (2015) 118e128
0
1 0 Sf ¼ @ thu A; thv
(17)
where t is the linear friction coefficient. Then the same numerical scheme as Eq. (14) can be applied with replacing Eq. (13) with
0 0 vSf ¼ @0 vq 0
0 t 0
1 0 0 A: t
(18)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The linear friction coefficient is t ¼ gn2 u2 þ v2 =h4=3 for Eq. (18). With assuming t as a constant, the Manning's friction can be treated as the linear friction in Gerris. The significance of this treatment will be investigated in the next section. 3. Numerical simulations and discussion The effect of the three different treatments for the friction term is investigated, calling Eqs. (13), (16) and (18) “nonlinear friction”, “Liang (2010)’s approach”, and “linearized friction”, respectively. In addition, the model performance and other practical issues, e.g., dynamic mesh refinement and boundary condition treatments using Gerris for real flood simulations are discussed. 3.1. Oscillations in a parabolic container with nonlinear friction A one-dimensional simple problem is numerically simulated in Fig. 1, using the three methods to handle the nonlinear friction. This benchmark test has been widely used to verify well-balanced numerical schemes (Liang and Marche, 2009; Popinet, 2011; An and Yu, 2012). The simulation domain has the topography of z ¼ h0(x/ a)2 where h0 ¼ 10 m and a ¼ 3000 m are constant parameters. The initial water profile is set to be inclined as shown in Fig. 1, with letting the water surface oscillate due to the bed topography and
121
friction. Four different Manning's coefficients of n ¼ 0.4, 0.2, 0.1, 0.05 are applied for comparison. Note that we only compare simulation results with using the three different approaches, given that there is no analytical solution with nonlinear friction, such as Manning's friction, although analytical solutions are given by Sampson et al. (2006) with linear friction. Simulation results in Fig. 1 show that Liang (2010)’s approach gives the water profiles between those computed using “linearised friction” and “nonlinear friction”, which is probably because Liang (2010)’s approach considers only the diagonal terms in Eq. (16) and ignores some of the nonlinear effects. The differences among the three treatments, however, become insignificant when the roughness coefficient gets smaller. The effect of the different treatment can hardly be seen in (c) and (d) with the roughness coefficient of 0.1 and 0.05, respectively. Given that floodplains usually do not have Manning's roughness coefficients of 0.2 or greater (Arcement and Schneider, 1989), it is expected that the three different treatments for the nonlinear friction will yield very similar results in real flood simulations. 3.2. Malpasset dam break The Malpasset dam break event is simulated to describe the model performance on real flood events. It is also verified that the difference between “nonlinear” and “linearized” friction treatment is not significant in a real flood simulation, which in general has the Manning's roughness coefficient under 0.1. This flood event has been widely used for 2D model verification (i.e., Hervouet, 2000; Valiani et al., 2002; Brufau et al., 2004; Kim and Cho, 2011; Kesserwani and Liang, 2012) given the reliable observation data. The Malpasset dam, located in the Reyran River Valley in France, was built for irrigation and drinking water storage. The dam formed a reservoir that contained over 5.5 107 m3 of water at maximum capacity. On December 2, 1959, the dam failed explosively due to intense rainfall which caused water levels in the reservoir to rise
Fig. 1. Oscillations in a parabolic container with using the different treatment of Manning's friction.
122
H. An et al. / Quaternary International 384 (2015) 118e128
rapidly. This dam break event created high waves of more than 40 m moving swiftly downstream. The dam break waves destroyed downstream villages, with causing 421 casualties and USD 68 million worth of infrastructural damages. After the flooding, a local police member surveyed maximum water level marks on the left and the right side of the Reyran River Valley. The 17 police surveyed points (P1eP17) are shown in Fig. 2. In addition, Fig. 2 describes the topography of the flow domain and the location of gauge points de France (S6eS14) in a scaled model performed by Electricite (EDF). The EDF laboratory built a non-distorted 1/400 scale model to better understand the flooding and calibrated the measurements against observations. The maximum water levels were observed at 14 gauges (S1eS14) in the physical model, of which the first five gauges (S1eS5) were located in the reservoir and omitted in Fig. 2. The Manning's friction coefficient (n) was estimated to be 0.025e0.033 sm1/3, and 0.033 sm1/3 is used in this paper following Valiani et al. (2002). For numerical calculation, the dam is treated as a straight line linking the two points, (4701.2 m, 4143.4 m) and (4655.5 m, 4392.1 m) in Fig. 2. The initial water level in the reservoir upstream from the dam is set to be 100 m and the rest of the domain is set dry. Quadtree meshes are adaptively refined around wave fronts and the river boundary with the minimum LEVEL 3(Dx z 2125 m) and the maximum LEVEL 9(Dx z 33 m). Gerris allows users to choose refinement criteria, which can be water depth, flow velocity, or others. Alternatively, users can set a complex refinement criterion given that Gerris can read C codes in a script file. For example, different ranges of LEVELs can be set depending on location, or a function of variables can be used as a refinement criterion using C codes. Based on experience, using the gradient of water level (h þ z) as the refinement criterion produces reliable results in flood simulations, which results that a cell is not refined at dried zones while refined at wetted zones when jVðh þ zÞjD > ε, where D ¼ Dx ¼ Dy is the cell size and ε is an adjustable parameter. The same flood event is simulated on the uniform grids of LEVEL 7, 8, 9, respectively, in order to compare the computational cost versus the accuracy of adaptive quadtree modeling to that of uniform grid modeling. Table 1 shows that the CPU times using the uniform grids of LEVEL 7, 8, 9 are 13, 253, and 1838 s; while the CPU time using the adaptive grid is 82 s, with saving approximately 95% of the computational time. Fig. 3 compares the simulated and the observed maximum water levels at the police survey points in (a) and at the gauging points of the EDF physical model in (b). The simulated maximum water levels using the adaptive grid or the
Fig. 3. Simulated and observed maximum water levels at (a) police survey points and de France (EDF). (b) gauging points in the physical model performed by Electricite
uniform grid of LEVEL 9 agree well with the observations, and are comparable to the results simulated by the other studies such as Valiani et al. (2002). On the other hand, the simulation results obtained using the uniform grids of LEVEL 7 and 8 are relatively poor at the both sides of the valley in Fig. 3(a), probably because of their insufficient resolution to represent steep topography. These results in Fig. 3 indicate that the simulation of this event requires LEVEL 9 (Dx z 33 m) or the finer grid resolution, which takes a lot of computational times as shown in Table 1. Fig. 3 and Table 1 demonstrate that the adaptive grid refinement with the minimum LEVEL 3 and the maximum LEVEL 9 produces the same results as the uniform grid of LEVEL 9 while saving approximately 95% of the CPU cost.
Fig. 2. Topography of the flow domain in the Malpasset dam break event. P1eP17 are police survey points and S6eS14 are measurement points in the physical model performed by de France (EDF). Electricite
H. An et al. / Quaternary International 384 (2015) 118e128 Table 1 Computational costs (Sec). Flood case
Level 7 Level 8 Level 9 Adaptive Adaptive/Finest (%)
Malpasset dam break 13 Baeksan levee failure e
253 1838 46,627 e
82a 6,717b
4.5 14.4
a With the minimum LEVEL 3(Dx z 2125 m) and the maximum LEVEL 9(Dx z 33 m). b With the minimum LEVEL 4(Dx z 390 m) and the maximum LEVEL 8(Dx z 33 m).
123
The effects of friction term discretization are compared in Fig. 4 where “linearized friction” represents Eq. (18) while “nonlinear friction” represents Eq. (13) presented in this paper. As expected from Section 3.1 and the fact that the Manning's roughness coefficient is below 0.05, there is no significant difference between two different treatments of the nonlinear friction term. This result is consistent with the result of Section 3.1. This simulation reconfirms that the linearized friction works well with the Manning's roughness coefficient for floodplains. Mass conservation is another important issue in flood simulation when the dynamical mesh refinement is applied, although in general the mass is well conserved without any numerical treatment when finite volume models are used. Gerris provides two options such as preserving mass conservation and lake-at-rest equilibrium (Popinet, 2011) for the dynamic mesh refinement. Based on our experiences, the mass conservation option produces more stable results in real flood simulations while the lake-at-rest equilibrium option may cause wrong solutions on complex topography. Thus, the mass conservation option is chosen with the dynamic mesh refinement for this simulation.
3.3. Baeksan levee failure
Fig. 4. Comparison of maximum water levels computed with “linearized” and “nonlinear” friction treatment. Adaptive grids are used.
Beaksan levee in the Nam River, Korea failed on August 10, 2002 due to intensive rainfall and consequent piping. Eighty (80) houses and 350 ha farmland were inundated. Consequently, 200 people had to be evacuated for 11 days. Although the damage of this event is not large, as opposed to the dramatically large Malpasset dam break case, this event is a general flood event caused by levee failures and suitable for the verification of flood inundation models. Fig. 5(a) shows the location of the failure and the flow simulation domain as well as three water level observation stations nearby. Topography of the flow domain is created in Fig. 5(b), with combining cross-sectional data of the river and the digital elevation map provided by National Geographic Information Institute of Korea. For the levee failure simulation, the height of the levee is added to the original topographic data at the break point. The width and the elevation of the failure are 10.3 m and 15 m, respectively, which were surveyed by the ministry of construction and transportation after the event. In Fig. 5(b), the points (G1eG6) are where simulation results are compared for model verification. The mesh resolution is adaptively adjusted following the same mesh refinement criterion as in the Malpasset dam break case, with the
Fig. 5. (a) Location and (b) topography of flow domain for the Beaksan levee failure event in Nakdong River basin, Korea.
124
H. An et al. / Quaternary International 384 (2015) 118e128
Fig. 6. (a) Upstream and (b) downstream boundary conditions calculated by using one-dimensional model HEC-RAC and water level data from three observation stations in Fig. 5(a).
Fig. 7. Land use of the flow domain in the Beaksan levee failure, which was surveyed in 1995e1999 by National Academy of Agriculture Science of Korea. Here rice paddy, dry field, orchard, and green house are categorized to farmland (A1), while grass land, forest, water, Etc are categorized to the other (A3) in Eq. (19), respectively.
Fig. 8. (a) Simulated water levels near the failure and (b) flow rates through the levee failure. G1 and G2 in (a) are located in the river and the protected lowland area, respectively. In (b), flow rates in the positive (þ) means that the water is flowing from the river into the protected lowland area.
H. An et al. / Quaternary International 384 (2015) 118e128
125
Fig. 9. (a) Simulated water depths, (b) spatial resolution of meshes, where the color blue represents level 4(Dx z 390 m) and the color red does level 8(Dx z 33 m), and (c) velocity vectors in the Beaksan levee failure simulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
minimum LEVEL 4(Dx z 390 m) and the maximum LEVEL 8(Dx z 33 m). The mass conservation option is chosen over the lake-at-rest equilibrium with the dynamic mesh refinement for this simulation as well. The initial condition is set to steady state using boundary conditions at t ¼ 0. The upstream and the downstream boundary conditions are given with the flow rate and the water level, respectively as shown in Fig. 6, which are calculated by using the widely used one-dimensional river routing model HEC-RAS and water level data from the observation stations. In addition, the Manning's coefficient is estimated calculating the area-weighted mean, following Ministry of Land, Transport and Maritime Affairs of Korea (2008) as:
n2 ¼
n21 A1 þ n22 A2 þ n23 A3 ; A1 þ A2 þ A3
(19)
where n1 ¼ 0.06, n2 ¼ 0.047, and n3 ¼ 0.05 are the Manning's coefficient of farmland, road, and the other area, and A1, A2, and A3 are the area of farmland, road, and the other, respectively. The estimated roughness coefficient of the flood plane is 0.058 sm1/3 given the land use map of the flow domain in Fig. 7, while the roughness coefficient of the river is set to be 0.023 sm1/3 referring to the river design standard of Korea (KWRA, 2005). In Fig. 7, rice paddy, dry field, orchard, and green house are categorized as farmland (A1), while grass land, forest, water, etc are categorized as other (A3) in Eq. (19), respectively. Simulated water levels and the flow rates with the levee failure are shown in Fig. 8. G1 and G2 in Fig. 8(a) are located in the river side and the protected lowland side of the failure, respectively. The flow in the positive (þ) direction in Fig. 8(b) means that the water is flowing from the river into the protected lowland area. Fig. 8(a) shows that the water level in the protected land (G2) rapidly
Fig. 10. Surveyed and simulated inundation areas.
126
H. An et al. / Quaternary International 384 (2015) 118e128
increases and reaches a peak about t ¼ 21 hr. After that, the water level in the protected land (G2) is higher than that in the river (G1), which causes the flow from the protected lowland area into the river in the negative () direction in Fig. 8(b). The peak flow rates in the positive (þ) and the negative () direction are approximately 210 m3/s and 60 m3/s, respectively. Flow in the positive (þ) direction develops rapidly while the flow in the negative () direction progresses slowly, probably due to the large difference in water levels along the river boundary at the time of the levee failure as shown in Fig. 6(a). The total water volume flowing out from the river to the protected lowland area through the levee failure is approximately 9.8 106 m3. Fig. 9(a) illustrates simulated water depths with time, or the flood propagation. The maximum water depth at the protected land is more than 5 m, and approximately 3 km2 of the land experiences inundation. The maximum water depth at t ¼ 48 hr is the lower than that at t ¼ 12 hr due to the flow in the negative () direction starting around t ¼ 21 hr. Fig. 9(b) shows the spatial resolution of the mesh. The grid is refined near wave fronts or along the river boundary, which allocates coarser grids on dry or inundated areas
during calculation. Fig. 9(c) displays simulated velocity vector fields with time, which confirms the fast flow into the protected lowland area right after the levee failure and then the slow flow to the river at later times observed in Fig. 8. The water in the river flows quickly toward the protected land at t ¼ 1,3 hr and propagates rapidly to the south at t ¼ 5 hr. At t ¼ 48 hr, the water slowly flows back to the river. After flooding, Gyeongnam Development Institute in Korea surveyed the inundated area and prepared an inundation trace map. Fig. 10 compares the inundation trace map and the simulated inundation area. Overall, the inundation area simulated using Gerris agrees well with the surveyed one, although the difference between surveyed and simulated is noticeable in mountainous area. The simulated inundation area agrees with the local topographic feature, while the surveyed inundation area apparently includes mountainous area that cannot have been inundated in Fig. 10. The inundation trace map overestimates the inundation area, given that the inundated area was estimated through site reconnaissance only without aerial photography.
Fig. 11. Comparison of water depths at G1eG6 in Fig. 11 simulated by using adaptive grids and uniform grids. In addition, the effect of inflow boundary treatments is compared. The adaptive_A represents that the inflow condition with converting the flow to the water level while the adaptive_B treats the inflow condition as a flux source.
H. An et al. / Quaternary International 384 (2015) 118e128
In order to compare the computational cost versus the accuracy of this sophisticated adaptive quadtree modeling to that of fine grid modeling, the same simulation is performed on the uniform grid of LEVEL 8. The CPU times with using the uniform and the adaptive grid are 46,627 and 6717 s, respectively in Table 1. Gerris seems to efficiently handle the grid adaption given that the adaptive quadtree grid reduces the CPU time approximately 7 times while producing very similar water level calculations at the 6 points (G1eG6) in Fig. 11. Some of the points such as G3 and G6 are initially located at the dry area and coarser grids are applied when adaptive quadtree grids are used, which makes the significant difference of water levels at initial times in Fig. 11. In addition, a negligible difference appears between the uniform grid and the adaptive quadtree grid modeling as the water recedes at later times, probably due to the adaptively refined meshes. Fig. 11 also shows the impact of boundary condition treatments. Gerris supports Dirichlet, Neumann, and reflective boundary conditions. The water level boundary such as Fig. 6(b) can be treated as the Dirichlet boundary condition in a straightforward way; while the inflow boundary such as Fig. 6(a) cannot be easily treated in shallow water type models. A way for Gerris to incorporate the inflow boundary is to search a water level given the specific inflow, with performing iterative calculations. Once the water level is determined, the Dirichlet boundary condition is applied using the estimated water level. Figs. 9e10 are the results using this approach, which converts the inflow boundary of Fig. 6(a) to the water level boundary. This method is, however, available only if there is one inflow boundary in the current version of Gerris. If there are more than two inflow boundaries in a simulation domain, this method is not available. Alternatively, users can treat the inflow boundary as a distributed flux source on specific areas near the boundary given that Gerris supports the explicit treatment of source terms as mentioned in Section 2.3. This alternative method is applicable to a variety of situations such as inflow and outflow conditions, although it is less physics-based. The two approaches
Fig. 12. Simulated velocity vectors 1 h after the Beaksan levee failure when the inflow condition of Fig. 6(a) is considered as a distributed flux source. The velocity vectors near the inflow boundary 1 h after the Beaksan levee failure in Fig. 9(c) is displayed in the red box for comparison. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
127
are compared in Fig. 11 where adaptive_A and adaptive_B represent the conversion of inflow to water level and the alternative approach, respectively. Fig. 11 shows that two methods yield the identical results. In addition, the velocity vectors 1 h after the levee failure calculated using the alternative method, with treating the inflow as a source are described in Fig. 12, which is almost identical to that shown in Fig. 9(c) except the flow near the inflow boundary. This difference is probably due to the conversion of the inflow boundary condition to a source. In conclusion, the alternative method to treat the inflow boundary can produce reasonable results in flood simulation, unless the flow behavior near the inflow boundary is important. The” linearized” friction treatment using the current version of Gerris gives the almost identical results to the “nonlinear” friction treatment of Manning's friction, as similar to Figs. 1(d) and 4, although not shown in this paper. 4. Conclusion The performance of the open source shallow water flow solver within Gerris Flow Solver for flood simulations is demonstrated and relevant practical issues are discussed. The Gerris Flow Solver has a number of favorable features for shallow water flow simulation, including adaptive mesh generation. The current version of Gerris, however, does not support implicit treatment for general friction source such as the Manning's friction, although the implicit treatment is required for numerical stability. Gerris supports the implicit discretization only for the linear friction. Thus, the applicability of linearized friction for flood simulation is investigated through the benchmark test. Simulation results indicate that the difference between the implicit treatments of the nonlinear friction and the linearization becomes insignificant when the Manning's roughness coefficient gets smaller. If the roughness coefficient is under 0.1, the linearization yields very similar results to that computed using the implicit discretization of the nonlinear friction, which means the current version of Gerris can deal with the Manning's friction in real flood simulation given that floodplains have roughness coefficients smaller than 0.1. Two real flood events, Malpasset dam break and Beaksan levee failure, were simulated in this paper. Malpasset dam break is widely used for the verification of flood models, given that the maximum water level marks were police-surveyed and a scale model was performed to better understand the event. Simulation results agree well with the observations and the measurements. The second case is a rural flood event caused by a levee failure. The simulated inundation area reasonably agrees with the inundation trace map, with noticeable mismatch observed only in mountainous areas. The mismatch between the simulation and the inundation trace map is probably from the fact that the inundation trace map overestimated the inundation areas given insufficient site survey. In addition, these two flood simulations show that the quadtree grid modeling is not only accurate but also effective and efficient. The adaptive grid modeling spent CPU times only 5%e15% of the finest uniform grid modeling while preserving the similar accuracy. The methods to treat friction terms and boundary conditions were investigated as well. Simulation results of the two real flood cases show that the linearization of the Manning's friction leads to the almost identical results to the implicit discretization of the nonlinear Manning's friction, which verify the conclusion from benchmark tests that the linearized and the nonlinear friction treatment give almost identical results with the Manning's coefficient for floodplains. Two different handlings of the inflow boundary condition were investigated. One method is to convert the inflow boundary into the water level and treat it using a Dirichlet boundary condition, while the other is to consider the inflow as a distributed flux source near the boundary. The
128
H. An et al. / Quaternary International 384 (2015) 118e128
comparison of simulation results using these two methods shows no significant differences except near the inflow boundary. In conclusion, the current version of Gerris is capable of representing flow dynamics of a broad range of floods. Acknowledgements This research was supported by a grant (11-TI-C06) from Advanced Water Management Research Program funded by Ministry of Land, Infrastructure and Transport of Korean government. This subject was also supported by Korea Ministry of Environment (MOE) as “GAIA Program-2014000540005”. References An, H., Yu, S., 2012. Well-balanced shallow water flow simulation on quadtree cut cell grids. Advances in Water Resources 39, 60e70. Arcement, G.J., Schneider, V.R., 1989. Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains. United States Geological Survey Water-Supply Paper 2339. US Government Printing Office, Washington, DC, USA. Audusse, E., Bristeanu, M.O., 2005. A well-balanced positivity preserving “secondorder” scheme for shallow water flows on unstructured meshes. Journal of Computational Physics 206 (1), 311e333. Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM Journal on Scientific Computing 25 (6), 2050e2065. Bates, P.D., De Roo, A.P.J., 2000. A simple raster-based model for flood inundation simulation. Journal of Hydrology 236 (1e2), 54e77. Berger, M., Colella, P., 1989. Local adaptive mesh refinement for shock hydrodynamics. Journal of Computational Physics 82 (1), 64e84. Berger, M., Oliger, J., 1984. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics 53 (3), 484e512. n, S., Jo zsa, J., 2001. Adaptive quadtree model of shallowBorthwick, A.G.L., Cruz Leo flow hydrodynamics. Journal of Hydraulic Research 39 (4), 413e424. zquez-Cendo n, M.E., 2004. Zero mass error using Brufau, P., García-Navarro, P., Va unsteady wetting-drying conditions in shallow flows over dry irregular topography. International Journal for Numerical Methods in Fluids 45 (10), 1047e1082. Caleffi, V., Valiani, A., 2009. Well-balanced bottom discontinuities treatment for high-order shallow water equations WENO scheme. Journal of Engineering Mechanics 135 (7), 684e696. Caleffi, V., Valiani, A., Bernini, A., 2006. Fourth-order balanced source term treatment in central WENO schemes for shallow water equations. Journal of Computational Physics 218 (1), 228e245. Chien, N.Q., 2008. Using open source software to simulate estuarine hydrodynamics. Journal of Water Resources and Environmental Engineering 23, 160e166. Ervine, D.A., Macleod, A.B., 1999. Modelling a river channel with distant floodbanks. Proceedings of the Institution of Civil Engineers: Water, Maritime and Energy 136 (1), 21e33. Fread, D.L., 1993. Flow routing. In: Maidment, D.R. (Ed.), Handbook of Applied Hydrology. McGraw-Hill, New York, pp. 10.1e10.36.
George, D.L., 2011. Adaptive finite volume methods with well balanced Riemann solvers for modeling floods in rugged terrain: Application to the Malpasset dam break flood (France, 1959). International Journal for Numerical Methods in Fluids 66 (8), 1000e1018. George, D.L., eVeque, R.J., 2008. High resolution methods and adaptive refinement for tsunami propagation and inundation. In: Hyperbolic Problems: Theory, Numerics, Applications. Springer, Berlin, Germany, pp. 541e550. Proc. 11th Int. Conf. on Hyperbolic Problems, Lyon, France, 17e21 July 2006. Hervouet, J.M., 2000. A high-resolution 2-D Dam-break model using parallelization. Hydrological Processes 14 (13), 2211e2230. Kesserwani, G., Liang, Q., 2012. Dynamically adaptive grid based discontinuous Galerkin shallow water model. Advances in Water Resources 37, 23e39. Kim, H., Cho, Y., 2011. Numerical model for flood routing with a Cartesian cut-cell domain. Journal of Hydraulic Research 49 (2), 205e212. Kim, D., Cho, Y., Kim, H., 2008. Well-balanced scheme between flux and source terms for computation of shallow-water equations over irregular bathymetry. Journal of Engineering Mechanics 134 (4), 277e290. Korea Water Resources Association(KWRA), 2005. River Design Standard of Korea. KWRA (in Korean). Lee, W.K., Borthwick, A.G.L., Taylor, P.H., 2011. A fast adaptive quadtree scheme for a two-layer shallow water model. Journal of Computational Physics 230, 4848e4870. Lee, D.E., An, H.U., Lee, G.H., Jung, K.S., 2013. Applicability evaluation of flood inundation analysis using quadtree grid-based model. Journal of Korea Water Resources Association 46 (6), 655e666 (in Korean). Liang, Q., 2010. Flood simulation using a well-balanced shallow flow model. Journal of Hydraulic Engineering, ASCE 136 (9), 669e675. Liang, Q., Marche, F., 2009. Numerical resolution of well-balanced shallow water equations with complex source terms. Advances in Water Resources 32 (6), 873e884. Liang, Q., Zang, J., Borthwick, A.G.L., Taylor, P., 2007. Shallow flow simulation on dynamically adaptive cut cell quadtree grids. International Journal for Numerical Methods in Fluids 53, 1777e1799. Mignot, E., Paquier, A., Haider, S., 2006. Modeling floods in a dense urban area using 2D shallow water equations. Journal of Hydrology 327, 186e199. Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics 190 (2), 572e600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228 (16), 5838e5866. Popinet, S., 2011. Quadtree-adaptive tsunami modelling. Ocean Dynamics 61, 1261e1285. Popinet, S., 2012. Adaptive modelling of long-distance wave propagation and finescale flooding during the Tohoku tsunami. Natural Hazards and Earth System Sciences 12 (4), 1213e1227. Popinet, S., Gorman, R.M., Rickard, G.J., Tolman, H.L., 2010. A quadtree-adaptive spectral wave model. Ocean Modelling 34, 36e49. Sampson, J., Easton, A., Singh, M., Moving boundary shallow water flow in parabolic bottom topography. ANZIAM Journal, 47, C373e387. Syme, B., 2011. Pros and Cons of 1D and 2D Modeling. In: Floodplain Management Association Conference. Sheraton San Diego Hotel & Marina, San Diego, USA. Valiani, A., Caleffi, V., Zanni, A., 2002. Case Study: Malpasset Dam-break simulation using a two-dimensional finite volume method. Journal of Hydraulic Engineering 128 (5), 460e472. Yoon, T.H., Kang, S.K., 2004. Finite volume model for two-dimensional shallow water flows on unstructured grids. Journal of Hydraulic Engineering 130 (7), 678e688.
Contents lists available at ScienceDirect
Quaternary International journal homepage: www.elsevier.com/locate/quaint
Analysis of an open source quadtree grid shallow water flow solver for flood simulation Hyunuk An a, Soonyoung Yu b, *, Giha Lee c, Yeonsu Kim d a
Department of Agricultural and Rural Engineering, Chungnam National University, Daejeon 305-764, Republic of Korea Research Institute for Social Criticality, Pusan National University, Busan 609-735, Republic of Korea c Department of Constructional Disaster Prevention Engineering, Kyungpook National University, Sangju 742-711, Republic of Korea d International Water Resources Research Institute, Chungnam National University, Daejeon 305-764, Republic of Korea b
a r t i c l e i n f o
a b s t r a c t
Article history: Available online 19 February 2015
The shallow water flow solver within Gerris Flow Solver implements second-order accurate Gudunov type numerical schemes, with preserving the balance of source and flux terms on quadtree cut cell grids. The solver provides flexible grid generation, with using adaptive quadtree grids and the cut cell method, which can improve computational efficiency in 2D simulations. The solver, however, cannot implicitly discretize nonlinear friction without linearization since the current version of Gerris supports the implicit discretization only for linear friction. We compare implicit treatment for Manning's friction to the linearization. Simulation results using a benchmark test show the difference between two methods is only significant when the roughness coefficient is unrealistically high. Two real flood events, Malpasset dam break in France and Baeksan levee failure in Korea, are simulated using the quadtree grid based shallow water model, with adaptively refining meshes near water fronts and river boundaries. Comparison of simulation results with observations and measurements demonstrates that the grid adaptation can save approximately 85%e95% of the computational cost while preserving the accuracy. The research result shows that the linearized friction using the current version of Gerris can be applied for flood simulation and multiple inflow boundaries can be treated as source terms. © 2015 Elsevier Ltd and INQUA. All rights reserved.
Keywords: Gerris Flood simulation Quadtree grid Bed friction Inflow boundary
1. Introduction Numerical modeling of shallow water flow is an important research topic for flood risk management. The simulation results play a significant role in national decision-making on flood prevention and control. One-dimensional shallow water models are quite often preferred in the field of engineering (e.g., FLUCOMP; Fread, 1993; Ervine and MacLeod, 1999; HEC-RAS; ISIS; MIKE11), mostly due to the fast computation (Syme, 2011). The rapid development in computer technologies, however, makes twodimensional shallow water flow models acceptable in practice. A number of 2D models have been developed (i.e., Bates and De Roo, 2000; Yoon and Kang, 2004; Mignot et al., 2006; George, 2011; Kim and Cho, 2011). In general, the 2D models are more accurate and reliable than 1D models for complex flow simulation given that the
* Corresponding author. E-mail addresses: [email protected] (H. An), [email protected] (S. Yu), [email protected] (G. Lee), [email protected] (Y. Kim). http://dx.doi.org/10.1016/j.quaint.2015.01.032 1040-6182/© 2015 Elsevier Ltd and INQUA. All rights reserved.
water is allowed to move only in the longitudinal direction in 1D models. 2D modeling is, however, still computationally heavy to engineers despite the advances in computing technologies, particularly for flood risk assessment or real time flood forecasting, which forces engineers to choose 1D models over 2D, giving up accuracy. One of the solutions to improve computational efficiency in 2D modeling is the implementation of adaptive mesh refinement (AMR), in which grid resolution is adjusted with following flow features. Generally, the finer grid resolution is achieved along wave fronts or complex geometry than within inundated area. The optimal grid resolution is time varying given that flows change with time. The method of adaptive meshing makes numerical simulations much more efficient (Berger and Oliger, 1984; Berger and Colella, 1989; Borthwick et al., 2001), and its adoption for numerical simulations becomes practice in engineering and science, including shallow water flow modeling. George and LeVeque (2008) used block-structured AMR for tsunami modelling. Their AMR algorithms allow nesting of rectangular subgrids of multiple levels and the subgrids evolve spatially and temporally. George (2011) applied the same model to simulate dam break flow and showed
H. An et al. / Quaternary International 384 (2015) 118e128
that the shallow water model with AMR is efficient and suitable for flood simulation. Alternatively, Liang et al. (2007) and Liang (2010) simulated inundation scenarios on adaptive quadtree grids. The adaptive quadtree grid is a way to carry out AMR and uses a quadtree and hierarchical data structure. Liang (2010) compared the performance of their model with the commercial code and showed that the shallow water model on adaptive quadtree grids is highly efficient while preserving accuracy. Popinet (2011, 2012) developed a quadtree adaptive tsunami model to solve shallow water equations within the Gerris Flow Solver framework (Popinet, 2003). Lee et al. (2011) presented a dynamically adaptive quadtree grid generation system for the solution of a two-dimensional two-layer shallow water model in order to improve computational efficiency. Gerris Flow Solver has the capability of adaptive quadtree grid generation. The Gerris solver was initially designed to solve the incompressible Euler equation, and then has been extended to multiphase incompressible NaviereStokes equations (Popinet, 2009), spectral wave models (Popinet et al., 2010), and shallow water equations to simulate tsunami propagation (Popinet, 2011, 2012). Recently, An and Yu (2012) applied the cut cell method to the shallow water flow solver of Gerris, which provides the model with more flexibility. In addition, the Gerris Flow Solver is open source software, and users can see the internal structure and even modify the codes for better representation. Moreover, the Gerris Flow Solver is well equipped with auxiliary programs such as GfsView for visualization, which most open source software are lack of as Chien (2008) mentioned. Gerris supports a variety of file formats such as Cartesian Grid Data (CGD), ESRI ASCII grid, xyz files to treat initial and boundary conditions, which facilitates specifying a time varying boundary condition. Gerris can produce outputs as text, image, MPG video, or KMZ files for Google Earth as well as GfsView files. Lee et al. (2013) compared the simple pre- and postprocessing in Gerris to those in FLUMEN (FLUvial Modeling Engine), and verified the computational efficiency of quadtree grids over triangular grids. The objective of this paper is to analyze the performance of the shallow water flow solver within Gerris Flow Solver for flood simulation in 2D through real flood cases. The motivation of this paper is the fact that there are few researches with applying the quadtree grid shallow water model to real flood events. Liang et al. (2007) and Liang (2010) simulated flood scenarios using a quadtree model for flood risk assessment. Their simulation results are in good agreement with those obtained by using a commercial model, called TUFLOW. In addition, Kesserwani and Liang (2012) simulated the Malpasset dam break event using the dynamically adaptive grid based shallow water model. The shallow water flow solver of Gerris has been intensively verified through real tsunami cases (Popinet, 2011, 2012) as well as a number of benchmark tests (Popinet et al., 2010; An and Yu, 2012). However, it is still informative to analyze the model performance on real flood events for potential Gerris users because Gerris has several favorable features as stated above. In addition, Gerris has several options to treat inflow boundaries as well as friction, and allows choosing mass conservation or lake-at-rest equilibrium when the dynamic adaptive mesh refinement is applied for flood simulations. The effects of those options are investigated through a benchmark test and real flood simulations in this paper.
119
broad range of flood simulations in practice. A hyperbolic conservation form of the equations is written as follows:
vq vf vg þ þ ¼ s; vt vx vy
(1)
where t denotes time, x and y are Cartesian coordinates, and q, f, g, and s are vectors representing conserved variables, fluxes in the xand y-direction, and source terms, respectively. By neglecting the friction term, the vectors can be written as:
0
h
0
1
B C q ¼ @ hu A; 0 B g¼@
hv hv
1 hu . 2 B 2 C f ¼ @ hu þ gh 2 A; 1
huv . C A; hu þ gh2 2 2
huv 0
1 0 B C s ¼ @ hgzx A;
(2)
hgzy
where h is water depth, u and v are depth-averaged velocity components in the x and y-direction, g is the acceleration of gravity, and zx(¼ vz/vx) and zy(¼ vz/vy) are bed slopes in the x and y-direction. The friction term will be added to s in Eq. (2) for completeness, without changing the fundamental structure of the equations in Section 2.3. 2.2. Spatial and temporal discretization Recently, a number of well-balanced numerical schemes have been proposed for shallow water flow equations (e.g., Audusse et al., 2004; Audusse and Bristeanu, 2005; Caleffi et al., 2006; Kim et al., 2008; Caleffi and Valiani, 2009) given that the gradient of water depths and the bed slope terms are often imbalanced in a discretized domain in case of irregular beds. The imbalance may cause unphysical perturbations when dryewet transitions occur or flows reach a steady state condition, which significantly affects the accuracy of flood prediction. In order to solve the problem of imbalance, Gerris implements the hydrostatic reconstruction technique proposed by Audusse et al. (2004), whose satisfactory performance has been verified in Popinet (2011) and An and Yu (2012). Following An and Yu (2012), the second-order accurate discretization on a two dimensional quadtree cut cell grid is written as:
Ai;j nþ1 qi;j qni;j þ F i;j ¼ S i;j þ Sci;j Dt
(3)
where Ai,j is the area of the cell (i, j), Sci,j is an additional source term to balance in the presence of a source, and Fi,j is the numerical flux. A MUSCL-type unsplit discretisation is combined with a predictorcorrector time-stepping scheme to construct a second-order accurate Godunov type scheme (Popinet, 2011). The predicator step computes a temporary cell-center value over the half time step Dt/ 2: nþ1=2
Ai;j qi;j
¼ Ai;j qni;j
Dt n F i;j S ni;j Sc ni;j ; 2
(4)
The corrector step updates the conservative values over a time step Dt as: 2. Numerical model 2.1. Governing equation Two-dimensional depth-integrated shallow water equations are used for flood event simulation. They have performed well for a
nþ1=2 nþ1=2 nþ1=2 : Ai;j qnþ1 ¼ Ai;j qni;j Dt F i;j S i;j Sc i;j i;j
(5)
The above scheme is second-order accurate in space and time, and the details of discretization and the well-balanced scheme are described in An and Yu (2012).
120
H. An et al. / Quaternary International 384 (2015) 118e128
where Dq ¼ qnþ1 qn. Ignoring the high-order terms and substituting Eq. (10) into Eq. (9) give
2.3. Friction term discretization Frictional effects should be incorporated to the right-hand side of Eq. (1) for real flood simulation. The Manning formula is considered as:
DðDtÞ ¼
*
(12)
0 0 B vSf2 B vSf B ¼ B vh B vq @ vSf 3
vh
0 vSf2 vðhuÞ vSf3 vðhuÞ
0 1 vSf2 C C vðhvÞ C C: C vSf A
(13)
3
vðhvÞ
n
Denoting Sf ðDtÞ¼ D1 ðDtÞSf n yields the implicit scheme written as:
¼ Ai;j qni;j
Dt n F i;j S ni;j Sc ni;j Ai;j Sf ni;j ; 2 nþ1=2 nþ1=2 nþ1=2 nþ1=2 n : ¼ Ai;j qi;j Dt F i;j S i;j Sc i;j Ai;j Sf i;j
nþ1=2
Ai;j qnþ1 i;j
n
where I is a unit matrix and
cretization methods for general source terms while implicit discretization for linear friction and Coriolis force. The explicit method is written as follows:
Ai;j qi;j
vSf vq
Dt n Dt n Dt nþ1=2 nþ1=2* F i;j S ni;j Scni;j ; qi;j ; ¼ qi;j þ Sf 2 2 2 * nþ1=2 nþ1=2 nþ1=2 nþ1=2 ; qnþ1 ¼ Ai;j qni;j Dt F i;j S i;j Sci;j ¼ qnþ1 þ DtSf ðDtÞ: i;j i;j
nþ1=2*
Ai;j qnþ1 i;j
I Dt
(6)
where n is the Manning's roughness coefficient. The addition of this friction term may cause the numerical instability in a discretized domain, in particular when a dryewet transition exists. In other words, if water depth becomes close to zero, the Manning's friction can result in an exaggerated force, which reverses the flow in a nonphysical manner. Therefore, the friction term should be properly treated. The current version of Gerris provides explicit dis-
Ai;j qi;j
(11)
where
0
1 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 1 B gn2 u u2 þ v2 C B C Sf1 B C 1=3 C; Sf ¼ @ Sf2 A ¼ B h B C qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C B Sf3 @ 2 2 2 A gn v u þ v h1=3
qnþ1 ¼ qn þDtD1 ðDtÞSf n ;
¼ Ai;j qni;j
(7) However, the above equations may cause the numerical instability although it can be straightforwardly applied. Therefore, implicit treatment is applied for numerical stability in this paper. Three implicit methods are implemented and compared because only the linear friction is treated implicitly in the current version of Gerris. The first method is a recently developed implicit discretization based on splitting method (Liang and Marche, 2009). The implicit discretization treats the friction term with using the splitting method, which is equivalent to solve the following equation:
dq ¼ Sf : dt
Liang and Marche (2009) presented this method in onedimensional spaces, and we fully extend their method in twodimension in Eqs. (11)e(14). The friction force in Eq. (14) can cause an exaggerated reverse flow, which makes no physical sense, although the friction terms are discretized implicitly. Therefore, the following condition is added:
n
2qn n ; if Sf > Dt nþ1=2 qnþ1=2 qn ; if Sf ¼ > Dt Dt
Sf ¼ qn Sf
nþ1=2
2 Dt
(15)
where Sf and q are momentum components of Sf and q, respectively. With the above condition, the maximum effect of friction forces is to stop the fluid. The second method is another implicit discretization presented in Liang (2010), who applied the method of Liang and Marche (2009) on two-dimensional equations in a directional manner, which gives
(8)
The Euler implicit discretization of the above equation yields
qnþ1 qn ¼ Sf nþ1 : Dt
(9)
The right-hand side of Eq. (9) can be represented by the Taylor series expansion as:
Sf nþ1 ¼ Sf n þ
(14)
vSf vq
n
ðDqÞ þ o Dq2 ;
(10)
00 B B0 vSf B ¼B B vq @ 0
0 vSf2 vðhuÞ 0
0
1
C C C C: C vSf3 A vðhvÞ 0
(16)
Note that Eq. (16) ignores the effect of off-diagonal entries as compared with Eq. (13). The third method is the linearization of friction terms. The current version of Gerris allows the implicit discretization of the linear friction as follows:
H. An et al. / Quaternary International 384 (2015) 118e128
0
1 0 Sf ¼ @ thu A; thv
(17)
where t is the linear friction coefficient. Then the same numerical scheme as Eq. (14) can be applied with replacing Eq. (13) with
0 0 vSf ¼ @0 vq 0
0 t 0
1 0 0 A: t
(18)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi The linear friction coefficient is t ¼ gn2 u2 þ v2 =h4=3 for Eq. (18). With assuming t as a constant, the Manning's friction can be treated as the linear friction in Gerris. The significance of this treatment will be investigated in the next section. 3. Numerical simulations and discussion The effect of the three different treatments for the friction term is investigated, calling Eqs. (13), (16) and (18) “nonlinear friction”, “Liang (2010)’s approach”, and “linearized friction”, respectively. In addition, the model performance and other practical issues, e.g., dynamic mesh refinement and boundary condition treatments using Gerris for real flood simulations are discussed. 3.1. Oscillations in a parabolic container with nonlinear friction A one-dimensional simple problem is numerically simulated in Fig. 1, using the three methods to handle the nonlinear friction. This benchmark test has been widely used to verify well-balanced numerical schemes (Liang and Marche, 2009; Popinet, 2011; An and Yu, 2012). The simulation domain has the topography of z ¼ h0(x/ a)2 where h0 ¼ 10 m and a ¼ 3000 m are constant parameters. The initial water profile is set to be inclined as shown in Fig. 1, with letting the water surface oscillate due to the bed topography and
121
friction. Four different Manning's coefficients of n ¼ 0.4, 0.2, 0.1, 0.05 are applied for comparison. Note that we only compare simulation results with using the three different approaches, given that there is no analytical solution with nonlinear friction, such as Manning's friction, although analytical solutions are given by Sampson et al. (2006) with linear friction. Simulation results in Fig. 1 show that Liang (2010)’s approach gives the water profiles between those computed using “linearised friction” and “nonlinear friction”, which is probably because Liang (2010)’s approach considers only the diagonal terms in Eq. (16) and ignores some of the nonlinear effects. The differences among the three treatments, however, become insignificant when the roughness coefficient gets smaller. The effect of the different treatment can hardly be seen in (c) and (d) with the roughness coefficient of 0.1 and 0.05, respectively. Given that floodplains usually do not have Manning's roughness coefficients of 0.2 or greater (Arcement and Schneider, 1989), it is expected that the three different treatments for the nonlinear friction will yield very similar results in real flood simulations. 3.2. Malpasset dam break The Malpasset dam break event is simulated to describe the model performance on real flood events. It is also verified that the difference between “nonlinear” and “linearized” friction treatment is not significant in a real flood simulation, which in general has the Manning's roughness coefficient under 0.1. This flood event has been widely used for 2D model verification (i.e., Hervouet, 2000; Valiani et al., 2002; Brufau et al., 2004; Kim and Cho, 2011; Kesserwani and Liang, 2012) given the reliable observation data. The Malpasset dam, located in the Reyran River Valley in France, was built for irrigation and drinking water storage. The dam formed a reservoir that contained over 5.5 107 m3 of water at maximum capacity. On December 2, 1959, the dam failed explosively due to intense rainfall which caused water levels in the reservoir to rise
Fig. 1. Oscillations in a parabolic container with using the different treatment of Manning's friction.
122
H. An et al. / Quaternary International 384 (2015) 118e128
rapidly. This dam break event created high waves of more than 40 m moving swiftly downstream. The dam break waves destroyed downstream villages, with causing 421 casualties and USD 68 million worth of infrastructural damages. After the flooding, a local police member surveyed maximum water level marks on the left and the right side of the Reyran River Valley. The 17 police surveyed points (P1eP17) are shown in Fig. 2. In addition, Fig. 2 describes the topography of the flow domain and the location of gauge points de France (S6eS14) in a scaled model performed by Electricite (EDF). The EDF laboratory built a non-distorted 1/400 scale model to better understand the flooding and calibrated the measurements against observations. The maximum water levels were observed at 14 gauges (S1eS14) in the physical model, of which the first five gauges (S1eS5) were located in the reservoir and omitted in Fig. 2. The Manning's friction coefficient (n) was estimated to be 0.025e0.033 sm1/3, and 0.033 sm1/3 is used in this paper following Valiani et al. (2002). For numerical calculation, the dam is treated as a straight line linking the two points, (4701.2 m, 4143.4 m) and (4655.5 m, 4392.1 m) in Fig. 2. The initial water level in the reservoir upstream from the dam is set to be 100 m and the rest of the domain is set dry. Quadtree meshes are adaptively refined around wave fronts and the river boundary with the minimum LEVEL 3(Dx z 2125 m) and the maximum LEVEL 9(Dx z 33 m). Gerris allows users to choose refinement criteria, which can be water depth, flow velocity, or others. Alternatively, users can set a complex refinement criterion given that Gerris can read C codes in a script file. For example, different ranges of LEVELs can be set depending on location, or a function of variables can be used as a refinement criterion using C codes. Based on experience, using the gradient of water level (h þ z) as the refinement criterion produces reliable results in flood simulations, which results that a cell is not refined at dried zones while refined at wetted zones when jVðh þ zÞjD > ε, where D ¼ Dx ¼ Dy is the cell size and ε is an adjustable parameter. The same flood event is simulated on the uniform grids of LEVEL 7, 8, 9, respectively, in order to compare the computational cost versus the accuracy of adaptive quadtree modeling to that of uniform grid modeling. Table 1 shows that the CPU times using the uniform grids of LEVEL 7, 8, 9 are 13, 253, and 1838 s; while the CPU time using the adaptive grid is 82 s, with saving approximately 95% of the computational time. Fig. 3 compares the simulated and the observed maximum water levels at the police survey points in (a) and at the gauging points of the EDF physical model in (b). The simulated maximum water levels using the adaptive grid or the
Fig. 3. Simulated and observed maximum water levels at (a) police survey points and de France (EDF). (b) gauging points in the physical model performed by Electricite
uniform grid of LEVEL 9 agree well with the observations, and are comparable to the results simulated by the other studies such as Valiani et al. (2002). On the other hand, the simulation results obtained using the uniform grids of LEVEL 7 and 8 are relatively poor at the both sides of the valley in Fig. 3(a), probably because of their insufficient resolution to represent steep topography. These results in Fig. 3 indicate that the simulation of this event requires LEVEL 9 (Dx z 33 m) or the finer grid resolution, which takes a lot of computational times as shown in Table 1. Fig. 3 and Table 1 demonstrate that the adaptive grid refinement with the minimum LEVEL 3 and the maximum LEVEL 9 produces the same results as the uniform grid of LEVEL 9 while saving approximately 95% of the CPU cost.
Fig. 2. Topography of the flow domain in the Malpasset dam break event. P1eP17 are police survey points and S6eS14 are measurement points in the physical model performed by de France (EDF). Electricite
H. An et al. / Quaternary International 384 (2015) 118e128 Table 1 Computational costs (Sec). Flood case
Level 7 Level 8 Level 9 Adaptive Adaptive/Finest (%)
Malpasset dam break 13 Baeksan levee failure e
253 1838 46,627 e
82a 6,717b
4.5 14.4
a With the minimum LEVEL 3(Dx z 2125 m) and the maximum LEVEL 9(Dx z 33 m). b With the minimum LEVEL 4(Dx z 390 m) and the maximum LEVEL 8(Dx z 33 m).
123
The effects of friction term discretization are compared in Fig. 4 where “linearized friction” represents Eq. (18) while “nonlinear friction” represents Eq. (13) presented in this paper. As expected from Section 3.1 and the fact that the Manning's roughness coefficient is below 0.05, there is no significant difference between two different treatments of the nonlinear friction term. This result is consistent with the result of Section 3.1. This simulation reconfirms that the linearized friction works well with the Manning's roughness coefficient for floodplains. Mass conservation is another important issue in flood simulation when the dynamical mesh refinement is applied, although in general the mass is well conserved without any numerical treatment when finite volume models are used. Gerris provides two options such as preserving mass conservation and lake-at-rest equilibrium (Popinet, 2011) for the dynamic mesh refinement. Based on our experiences, the mass conservation option produces more stable results in real flood simulations while the lake-at-rest equilibrium option may cause wrong solutions on complex topography. Thus, the mass conservation option is chosen with the dynamic mesh refinement for this simulation.
3.3. Baeksan levee failure
Fig. 4. Comparison of maximum water levels computed with “linearized” and “nonlinear” friction treatment. Adaptive grids are used.
Beaksan levee in the Nam River, Korea failed on August 10, 2002 due to intensive rainfall and consequent piping. Eighty (80) houses and 350 ha farmland were inundated. Consequently, 200 people had to be evacuated for 11 days. Although the damage of this event is not large, as opposed to the dramatically large Malpasset dam break case, this event is a general flood event caused by levee failures and suitable for the verification of flood inundation models. Fig. 5(a) shows the location of the failure and the flow simulation domain as well as three water level observation stations nearby. Topography of the flow domain is created in Fig. 5(b), with combining cross-sectional data of the river and the digital elevation map provided by National Geographic Information Institute of Korea. For the levee failure simulation, the height of the levee is added to the original topographic data at the break point. The width and the elevation of the failure are 10.3 m and 15 m, respectively, which were surveyed by the ministry of construction and transportation after the event. In Fig. 5(b), the points (G1eG6) are where simulation results are compared for model verification. The mesh resolution is adaptively adjusted following the same mesh refinement criterion as in the Malpasset dam break case, with the
Fig. 5. (a) Location and (b) topography of flow domain for the Beaksan levee failure event in Nakdong River basin, Korea.
124
H. An et al. / Quaternary International 384 (2015) 118e128
Fig. 6. (a) Upstream and (b) downstream boundary conditions calculated by using one-dimensional model HEC-RAC and water level data from three observation stations in Fig. 5(a).
Fig. 7. Land use of the flow domain in the Beaksan levee failure, which was surveyed in 1995e1999 by National Academy of Agriculture Science of Korea. Here rice paddy, dry field, orchard, and green house are categorized to farmland (A1), while grass land, forest, water, Etc are categorized to the other (A3) in Eq. (19), respectively.
Fig. 8. (a) Simulated water levels near the failure and (b) flow rates through the levee failure. G1 and G2 in (a) are located in the river and the protected lowland area, respectively. In (b), flow rates in the positive (þ) means that the water is flowing from the river into the protected lowland area.
H. An et al. / Quaternary International 384 (2015) 118e128
125
Fig. 9. (a) Simulated water depths, (b) spatial resolution of meshes, where the color blue represents level 4(Dx z 390 m) and the color red does level 8(Dx z 33 m), and (c) velocity vectors in the Beaksan levee failure simulation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
minimum LEVEL 4(Dx z 390 m) and the maximum LEVEL 8(Dx z 33 m). The mass conservation option is chosen over the lake-at-rest equilibrium with the dynamic mesh refinement for this simulation as well. The initial condition is set to steady state using boundary conditions at t ¼ 0. The upstream and the downstream boundary conditions are given with the flow rate and the water level, respectively as shown in Fig. 6, which are calculated by using the widely used one-dimensional river routing model HEC-RAS and water level data from the observation stations. In addition, the Manning's coefficient is estimated calculating the area-weighted mean, following Ministry of Land, Transport and Maritime Affairs of Korea (2008) as:
n2 ¼
n21 A1 þ n22 A2 þ n23 A3 ; A1 þ A2 þ A3
(19)
where n1 ¼ 0.06, n2 ¼ 0.047, and n3 ¼ 0.05 are the Manning's coefficient of farmland, road, and the other area, and A1, A2, and A3 are the area of farmland, road, and the other, respectively. The estimated roughness coefficient of the flood plane is 0.058 sm1/3 given the land use map of the flow domain in Fig. 7, while the roughness coefficient of the river is set to be 0.023 sm1/3 referring to the river design standard of Korea (KWRA, 2005). In Fig. 7, rice paddy, dry field, orchard, and green house are categorized as farmland (A1), while grass land, forest, water, etc are categorized as other (A3) in Eq. (19), respectively. Simulated water levels and the flow rates with the levee failure are shown in Fig. 8. G1 and G2 in Fig. 8(a) are located in the river side and the protected lowland side of the failure, respectively. The flow in the positive (þ) direction in Fig. 8(b) means that the water is flowing from the river into the protected lowland area. Fig. 8(a) shows that the water level in the protected land (G2) rapidly
Fig. 10. Surveyed and simulated inundation areas.
126
H. An et al. / Quaternary International 384 (2015) 118e128
increases and reaches a peak about t ¼ 21 hr. After that, the water level in the protected land (G2) is higher than that in the river (G1), which causes the flow from the protected lowland area into the river in the negative () direction in Fig. 8(b). The peak flow rates in the positive (þ) and the negative () direction are approximately 210 m3/s and 60 m3/s, respectively. Flow in the positive (þ) direction develops rapidly while the flow in the negative () direction progresses slowly, probably due to the large difference in water levels along the river boundary at the time of the levee failure as shown in Fig. 6(a). The total water volume flowing out from the river to the protected lowland area through the levee failure is approximately 9.8 106 m3. Fig. 9(a) illustrates simulated water depths with time, or the flood propagation. The maximum water depth at the protected land is more than 5 m, and approximately 3 km2 of the land experiences inundation. The maximum water depth at t ¼ 48 hr is the lower than that at t ¼ 12 hr due to the flow in the negative () direction starting around t ¼ 21 hr. Fig. 9(b) shows the spatial resolution of the mesh. The grid is refined near wave fronts or along the river boundary, which allocates coarser grids on dry or inundated areas
during calculation. Fig. 9(c) displays simulated velocity vector fields with time, which confirms the fast flow into the protected lowland area right after the levee failure and then the slow flow to the river at later times observed in Fig. 8. The water in the river flows quickly toward the protected land at t ¼ 1,3 hr and propagates rapidly to the south at t ¼ 5 hr. At t ¼ 48 hr, the water slowly flows back to the river. After flooding, Gyeongnam Development Institute in Korea surveyed the inundated area and prepared an inundation trace map. Fig. 10 compares the inundation trace map and the simulated inundation area. Overall, the inundation area simulated using Gerris agrees well with the surveyed one, although the difference between surveyed and simulated is noticeable in mountainous area. The simulated inundation area agrees with the local topographic feature, while the surveyed inundation area apparently includes mountainous area that cannot have been inundated in Fig. 10. The inundation trace map overestimates the inundation area, given that the inundated area was estimated through site reconnaissance only without aerial photography.
Fig. 11. Comparison of water depths at G1eG6 in Fig. 11 simulated by using adaptive grids and uniform grids. In addition, the effect of inflow boundary treatments is compared. The adaptive_A represents that the inflow condition with converting the flow to the water level while the adaptive_B treats the inflow condition as a flux source.
H. An et al. / Quaternary International 384 (2015) 118e128
In order to compare the computational cost versus the accuracy of this sophisticated adaptive quadtree modeling to that of fine grid modeling, the same simulation is performed on the uniform grid of LEVEL 8. The CPU times with using the uniform and the adaptive grid are 46,627 and 6717 s, respectively in Table 1. Gerris seems to efficiently handle the grid adaption given that the adaptive quadtree grid reduces the CPU time approximately 7 times while producing very similar water level calculations at the 6 points (G1eG6) in Fig. 11. Some of the points such as G3 and G6 are initially located at the dry area and coarser grids are applied when adaptive quadtree grids are used, which makes the significant difference of water levels at initial times in Fig. 11. In addition, a negligible difference appears between the uniform grid and the adaptive quadtree grid modeling as the water recedes at later times, probably due to the adaptively refined meshes. Fig. 11 also shows the impact of boundary condition treatments. Gerris supports Dirichlet, Neumann, and reflective boundary conditions. The water level boundary such as Fig. 6(b) can be treated as the Dirichlet boundary condition in a straightforward way; while the inflow boundary such as Fig. 6(a) cannot be easily treated in shallow water type models. A way for Gerris to incorporate the inflow boundary is to search a water level given the specific inflow, with performing iterative calculations. Once the water level is determined, the Dirichlet boundary condition is applied using the estimated water level. Figs. 9e10 are the results using this approach, which converts the inflow boundary of Fig. 6(a) to the water level boundary. This method is, however, available only if there is one inflow boundary in the current version of Gerris. If there are more than two inflow boundaries in a simulation domain, this method is not available. Alternatively, users can treat the inflow boundary as a distributed flux source on specific areas near the boundary given that Gerris supports the explicit treatment of source terms as mentioned in Section 2.3. This alternative method is applicable to a variety of situations such as inflow and outflow conditions, although it is less physics-based. The two approaches
Fig. 12. Simulated velocity vectors 1 h after the Beaksan levee failure when the inflow condition of Fig. 6(a) is considered as a distributed flux source. The velocity vectors near the inflow boundary 1 h after the Beaksan levee failure in Fig. 9(c) is displayed in the red box for comparison. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
127
are compared in Fig. 11 where adaptive_A and adaptive_B represent the conversion of inflow to water level and the alternative approach, respectively. Fig. 11 shows that two methods yield the identical results. In addition, the velocity vectors 1 h after the levee failure calculated using the alternative method, with treating the inflow as a source are described in Fig. 12, which is almost identical to that shown in Fig. 9(c) except the flow near the inflow boundary. This difference is probably due to the conversion of the inflow boundary condition to a source. In conclusion, the alternative method to treat the inflow boundary can produce reasonable results in flood simulation, unless the flow behavior near the inflow boundary is important. The” linearized” friction treatment using the current version of Gerris gives the almost identical results to the “nonlinear” friction treatment of Manning's friction, as similar to Figs. 1(d) and 4, although not shown in this paper. 4. Conclusion The performance of the open source shallow water flow solver within Gerris Flow Solver for flood simulations is demonstrated and relevant practical issues are discussed. The Gerris Flow Solver has a number of favorable features for shallow water flow simulation, including adaptive mesh generation. The current version of Gerris, however, does not support implicit treatment for general friction source such as the Manning's friction, although the implicit treatment is required for numerical stability. Gerris supports the implicit discretization only for the linear friction. Thus, the applicability of linearized friction for flood simulation is investigated through the benchmark test. Simulation results indicate that the difference between the implicit treatments of the nonlinear friction and the linearization becomes insignificant when the Manning's roughness coefficient gets smaller. If the roughness coefficient is under 0.1, the linearization yields very similar results to that computed using the implicit discretization of the nonlinear friction, which means the current version of Gerris can deal with the Manning's friction in real flood simulation given that floodplains have roughness coefficients smaller than 0.1. Two real flood events, Malpasset dam break and Beaksan levee failure, were simulated in this paper. Malpasset dam break is widely used for the verification of flood models, given that the maximum water level marks were police-surveyed and a scale model was performed to better understand the event. Simulation results agree well with the observations and the measurements. The second case is a rural flood event caused by a levee failure. The simulated inundation area reasonably agrees with the inundation trace map, with noticeable mismatch observed only in mountainous areas. The mismatch between the simulation and the inundation trace map is probably from the fact that the inundation trace map overestimated the inundation areas given insufficient site survey. In addition, these two flood simulations show that the quadtree grid modeling is not only accurate but also effective and efficient. The adaptive grid modeling spent CPU times only 5%e15% of the finest uniform grid modeling while preserving the similar accuracy. The methods to treat friction terms and boundary conditions were investigated as well. Simulation results of the two real flood cases show that the linearization of the Manning's friction leads to the almost identical results to the implicit discretization of the nonlinear Manning's friction, which verify the conclusion from benchmark tests that the linearized and the nonlinear friction treatment give almost identical results with the Manning's coefficient for floodplains. Two different handlings of the inflow boundary condition were investigated. One method is to convert the inflow boundary into the water level and treat it using a Dirichlet boundary condition, while the other is to consider the inflow as a distributed flux source near the boundary. The
128
H. An et al. / Quaternary International 384 (2015) 118e128
comparison of simulation results using these two methods shows no significant differences except near the inflow boundary. In conclusion, the current version of Gerris is capable of representing flow dynamics of a broad range of floods. Acknowledgements This research was supported by a grant (11-TI-C06) from Advanced Water Management Research Program funded by Ministry of Land, Infrastructure and Transport of Korean government. This subject was also supported by Korea Ministry of Environment (MOE) as “GAIA Program-2014000540005”. References An, H., Yu, S., 2012. Well-balanced shallow water flow simulation on quadtree cut cell grids. Advances in Water Resources 39, 60e70. Arcement, G.J., Schneider, V.R., 1989. Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains. United States Geological Survey Water-Supply Paper 2339. US Government Printing Office, Washington, DC, USA. Audusse, E., Bristeanu, M.O., 2005. A well-balanced positivity preserving “secondorder” scheme for shallow water flows on unstructured meshes. Journal of Computational Physics 206 (1), 311e333. Audusse, E., Bouchut, F., Bristeau, M.O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM Journal on Scientific Computing 25 (6), 2050e2065. Bates, P.D., De Roo, A.P.J., 2000. A simple raster-based model for flood inundation simulation. Journal of Hydrology 236 (1e2), 54e77. Berger, M., Colella, P., 1989. Local adaptive mesh refinement for shock hydrodynamics. Journal of Computational Physics 82 (1), 64e84. Berger, M., Oliger, J., 1984. Adaptive mesh refinement for hyperbolic partial differential equations. Journal of Computational Physics 53 (3), 484e512. n, S., Jo zsa, J., 2001. Adaptive quadtree model of shallowBorthwick, A.G.L., Cruz Leo flow hydrodynamics. Journal of Hydraulic Research 39 (4), 413e424. zquez-Cendo n, M.E., 2004. Zero mass error using Brufau, P., García-Navarro, P., Va unsteady wetting-drying conditions in shallow flows over dry irregular topography. International Journal for Numerical Methods in Fluids 45 (10), 1047e1082. Caleffi, V., Valiani, A., 2009. Well-balanced bottom discontinuities treatment for high-order shallow water equations WENO scheme. Journal of Engineering Mechanics 135 (7), 684e696. Caleffi, V., Valiani, A., Bernini, A., 2006. Fourth-order balanced source term treatment in central WENO schemes for shallow water equations. Journal of Computational Physics 218 (1), 228e245. Chien, N.Q., 2008. Using open source software to simulate estuarine hydrodynamics. Journal of Water Resources and Environmental Engineering 23, 160e166. Ervine, D.A., Macleod, A.B., 1999. Modelling a river channel with distant floodbanks. Proceedings of the Institution of Civil Engineers: Water, Maritime and Energy 136 (1), 21e33. Fread, D.L., 1993. Flow routing. In: Maidment, D.R. (Ed.), Handbook of Applied Hydrology. McGraw-Hill, New York, pp. 10.1e10.36.
George, D.L., 2011. Adaptive finite volume methods with well balanced Riemann solvers for modeling floods in rugged terrain: Application to the Malpasset dam break flood (France, 1959). International Journal for Numerical Methods in Fluids 66 (8), 1000e1018. George, D.L., eVeque, R.J., 2008. High resolution methods and adaptive refinement for tsunami propagation and inundation. In: Hyperbolic Problems: Theory, Numerics, Applications. Springer, Berlin, Germany, pp. 541e550. Proc. 11th Int. Conf. on Hyperbolic Problems, Lyon, France, 17e21 July 2006. Hervouet, J.M., 2000. A high-resolution 2-D Dam-break model using parallelization. Hydrological Processes 14 (13), 2211e2230. Kesserwani, G., Liang, Q., 2012. Dynamically adaptive grid based discontinuous Galerkin shallow water model. Advances in Water Resources 37, 23e39. Kim, H., Cho, Y., 2011. Numerical model for flood routing with a Cartesian cut-cell domain. Journal of Hydraulic Research 49 (2), 205e212. Kim, D., Cho, Y., Kim, H., 2008. Well-balanced scheme between flux and source terms for computation of shallow-water equations over irregular bathymetry. Journal of Engineering Mechanics 134 (4), 277e290. Korea Water Resources Association(KWRA), 2005. River Design Standard of Korea. KWRA (in Korean). Lee, W.K., Borthwick, A.G.L., Taylor, P.H., 2011. A fast adaptive quadtree scheme for a two-layer shallow water model. Journal of Computational Physics 230, 4848e4870. Lee, D.E., An, H.U., Lee, G.H., Jung, K.S., 2013. Applicability evaluation of flood inundation analysis using quadtree grid-based model. Journal of Korea Water Resources Association 46 (6), 655e666 (in Korean). Liang, Q., 2010. Flood simulation using a well-balanced shallow flow model. Journal of Hydraulic Engineering, ASCE 136 (9), 669e675. Liang, Q., Marche, F., 2009. Numerical resolution of well-balanced shallow water equations with complex source terms. Advances in Water Resources 32 (6), 873e884. Liang, Q., Zang, J., Borthwick, A.G.L., Taylor, P., 2007. Shallow flow simulation on dynamically adaptive cut cell quadtree grids. International Journal for Numerical Methods in Fluids 53, 1777e1799. Mignot, E., Paquier, A., Haider, S., 2006. Modeling floods in a dense urban area using 2D shallow water equations. Journal of Hydrology 327, 186e199. Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. Journal of Computational Physics 190 (2), 572e600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. Journal of Computational Physics 228 (16), 5838e5866. Popinet, S., 2011. Quadtree-adaptive tsunami modelling. Ocean Dynamics 61, 1261e1285. Popinet, S., 2012. Adaptive modelling of long-distance wave propagation and finescale flooding during the Tohoku tsunami. Natural Hazards and Earth System Sciences 12 (4), 1213e1227. Popinet, S., Gorman, R.M., Rickard, G.J., Tolman, H.L., 2010. A quadtree-adaptive spectral wave model. Ocean Modelling 34, 36e49. Sampson, J., Easton, A., Singh, M., Moving boundary shallow water flow in parabolic bottom topography. ANZIAM Journal, 47, C373e387. Syme, B., 2011. Pros and Cons of 1D and 2D Modeling. In: Floodplain Management Association Conference. Sheraton San Diego Hotel & Marina, San Diego, USA. Valiani, A., Caleffi, V., Zanni, A., 2002. Case Study: Malpasset Dam-break simulation using a two-dimensional finite volume method. Journal of Hydraulic Engineering 128 (5), 460e472. Yoon, T.H., Kang, S.K., 2004. Finite volume model for two-dimensional shallow water flows on unstructured grids. Journal of Hydraulic Engineering 130 (7), 678e688.