Comparison of Fast Shallow-Water Schemes on Real-World Floods

: Two-dimensional shallow-water schemes on Cartesian grids are amendable for graphics processing units and thus a convenient choice for fast flood simulations. A comparison of recent schemes and validation of important use cases is essential for developers and practitioners working with flood simulation tools. In this paper, we discuss three state-of-the-art shallow-water schemes: a first-order upwind scheme, a second-order upwind scheme, and a second-order central-upwind scheme. We analyze the advantages and disadvantages of each scheme on historical Danube river floods at three regions in Austria. We study the Lobau region as a floodplain with several small channels, the Wachau region with the meandering Danube in a steep valley, and the Marchfeld region located at the river confluence of March and Danube. The validation case studies show that the second-order schemes provide better estimates of the water levels than the first-order scheme. Still, the first order scheme is useful because it offers fast simulations and reasonable results at higher resolutions. The best trade-off between accuracy and computational effort for simulating river floods is provided by the second-order upwind scheme. DOI: 10.1061/(ASCE)HY.1943-7900.0001657. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/. in water extents between the second-order BH and HW schemes is caused by different bathymetry reconstructions.


Introduction
In flood management, there is a high demand for simulation tools that are capable of providing accurate and fast flood predictions.
To deal with the risks, it is important to not only identify the regions vulnerable to inundation, but also to estimate arrival times and maximum water levels. Such detailed results of flood simulations can provide valuable insights for practitioners. Unfortunately, flood simulation models are still not in widely used by decision makers in flood management. Leskens et al. (2014) give reasons for limited use of models for decision-making. This is mainly because of the discrepancies between what information is demanded and what can actually be offered by models in terms of output and accompanying uncertainties. There are also delays and constraints in the exchange of model information through the network of participants. People who develop new models are often not the ones who use them and are unaware of the needs of the actual users. For real-world simulations, long computation times and inflexible model setup are the main bottlenecks. Shallow-water models that overcome these obstacles can be used to interactively provide real-time guidance for decision makers (Waser et al. 2014). To study what if-scenarios, dynamic changes in the model setup are necessary, e.g., removal or insertion of flood protection measures. Typical use cases include the investigation of levee failures, e.g., a breach that widens over time. Furthermore, flows over complex terrains and overtoppings of flood protection measures should be handled by the numerical scheme. Thus, important requirements for a flood simulation tool are reliability, speed, and stability of the underlying numerical solver, as well as interactivity and flexibility when setting up the model.
In most flood simulation tools [HEC-RAS (Brunner 2010), TELEMAC (Galland et al. 1991) and Visdom (Visdom n.d.;Cornel et al. 2015)] the shallow-water equations (SWEs) are the underlying model equations. They are based on the assumption that the horizontal length scale is large compared to the vertical length scale. The SWEs can be derived by depth-averaging the Navier-Stokes and continuity equations (Temam 1984). They provide plausible and reliable results of water levels and wave arrival times for events Ph.D. Student, VRVis Zentrum für Virtual Reality und Visualisierung Forschungs-GmbH, Donau-City-Strasse 11, Vienna 1220, Austria. Email: cornel@vrvis.at such as river floods and dam breaks (Audusse et al. 2004;Brodtkorb et al. 2012;de la Asunción et al. 2013;Hervouet and Petitjean 1999;Liang and Marche 2009;Russo 2005). A common numerical method for solving SWEs is the finite volume method (FVM) (Godlewski and Raviart 1996;LeVeque 1992;LeVeque 2002;Toro 2001). For the FVM, a computational grid has to be chosen as a basis for the spatial discretization, which is usually an unstructured triangular mesh or a regular grid. Unstructured triangular meshes are able to incorporate complex geometries at the expense of worse performance because of additional complexity in both the mesh and the implementation. In contrast, rectangular grids are not able to exactly represent nonaligned topographic features. However, techniques exist that allow us to overcome this limitation, e.g., cut-cells (Liang and Borthwick 2008;Meinke et al. 2013). The striking advantage of Cartesian grids is their suitability for a straightforward parallelization on GPUs (Brodtkorb et al. 2012;Horváth et al. 2016;Vacondio et al. 2016). The advent of GPU computing has reduced computation times by a factor of up to 100 compared to conventional models (Casulli and Stelling 2013;Horváth et al. 2016;Vacondio et al. 2016). This speed-up in computation time does not only allow for smaller cell sizes and more accurate simulation runs, but also decreases model uncertainties as ensemble simulation becomes a viable option.
In recent years, some properties were proven to be essential for the stability of numerical schemes. The scheme should be positivity-preserving, that is, water depths should remain nonnegative at any point and any time. Furthermore, it should be capable of handling wet-dry boundaries, which is a challenging task. To simulate real-world floods over time periods of days or even weeks, it is important for the solution to remain stable for the duration, e.g., the scheme should not produce spurious velocity oscillations at shoal zones. Problems characterized by strong discontinuities, e.g., dam breaks or overtoppings of flood protection measures, should be captured accurately. Finally, the scheme should be well-balanced, i.e., it should be capable of balancing source and numerical flux terms for simple stationary solutions such as a still water. Most of the recently developed schemes are well-balanced (Audusse et al. 2004(Audusse et al. , 2015Bollermann et al. 2013;Durran 2010;Horváth et al. 2015;Hou et al. 2014;Li et al. 2013;Noelle et al. 2006;Russo 2005;Xing and Shu 2005).
In this paper, we compare three state-of-the-art shallow-water schemes: the first-order scheme of Chen and Noelle (CN) (Chen and Noelle 2017), and the second-order schemes of Buttinger et al.  . The predecessor of the CN and BH schemes presented by Audusse et al. (2004) is a popular scheme both in the academic world and outside it, and it can also be found in the TELEMAC-MASCARET system (Hervouet and Ata 2017). The second-order scheme by Kurganov and Petrova, which is the predecessor of the HW scheme, is also a popular scheme, especially for GPU-based solvers. The CN, BH, and HW schemes presented in this paper improve their predecessors by fixing the sources of numerical instabilities due to error accumulation and make long large-scale simulations possible.
We focus on schemes that are suitable for efficient GPU implementation, since in addition to accuracy, run time is of key importance in interactive decision support. Today's GPUs are very powerful, but they are still very limited by the available resources. On GPUs, the work is subdivided among several smaller computation units, each of which has limited resources. With this information in mind, we have to select schemes and design algorithms that do not overuse the available resources to avoid performance degradations. The chosen schemes are satisfactory from this computational point of view, while providing important stability properties at the same time, i.e., they are well-balanced and positivity-preserving.
The advantages and disadvantages of the schemes are demonstrated by simulating three historical flood scenarios of the Danube in Austria. We discuss in detail the set up of each case and the uncertainties that can affect the solutions. The presented real-world validation cases have different terrain features, ranging from relatively flat lowlands to river valleys with highly varying bottom topography. We validate the schemes using the Danube river flooding of 2011 in the Lobau national park, the Danube river flooding of 2013 in the Wachau valley, and in the same flooding of 2013 in Marchfeld.
The paper is organized as follows. First, we discuss the SWEs and the numerical schemes. Then, we present the numerical results obtained with the three different schemes and compare them with measured data. Finally, we conclude by sharing our findings.
This paper goes beyond the existing literature with the following key points: • assessment of the accuracy of flood extent predictions resulting from the CN, BH, and HW schemes, • comparison of the CN, BH, and HW schemes regarding performance in large-scale scenarios, • validation of shallow-water schemes on historical floods.

Model and Numerical Methods
In this section, we summarize the underlying numerical theory of the SWEs that we compare and validate. The hyperbolic conservation law described by the two-dimensional SWEs, also called the Saint-Venant system, can be written as where h = water height; hu = discharge along the x-axis; hv = discharge along the y-axis; u and v are the average flow velocities along the respective axes; g = gravitational constant; and B = bathymetry. The Chézy friction coefficient C is defined as C ¼ h 1=6 =n, where n = Manning roughness coefficient. Subscripts represent partial derivatives, e.g., U t stands for ∂U ∂t . In vector form the system can be written as where U ¼ ½h; hu; hv T is the vector of conserved variables, F and G are flux functions. S is the bed slope source term and models the fluid's acceleration due to the gravitational forces.
To provide realistic water flow, a bed friction source term S f is included.

Discretization
The FVM is chosen for the spatial discretization on top of a uniform grid with cell size Δx × Δy. The starting point for an FVM scheme are cell averagesŪ j;k ¼ ½h j;k ;w j;k ;ū j;k ;v j;k T of the conserved variables U j;k for a finite volume cell C j;k [ Fig. 1(a)]. For the numerical solution of the SWEs [Eq.
(2)], the discretized fluxes F jþ 1 2 ;k and G jþ 1 2 ;k at the interfaces and an appropriate source term discretization S j;k need to be specified [Fig. 1(d)]. We choose the Harten-Lax-van Leer (HLL) flux (Harten et al. 1983) as a flux function. The HLL flux function F jAE 1 2 ;k ðU − jþ 1 2 ;k ; U þ jþ 1 2 ;k Þ approximately solves a Riemann problem given by the two interface states: Care has to be taken when deriving those interface values to ensure nonnegativity of the water depths. To obtain interface point values, reconstruction steps need to be applied to the cell averages and the bathymetry. Also the bed slope source term discretization S j;k needs to be tailored to the reconstruction to achieve balance. Both the reconstruction and bed slope source terms depend on the scheme and are discussed in detail in "Schemes."

Time Integration
For the first-order CN scheme, an explicit Euler time integrator is used:Ū E;nþ1 j;k ¼Ū n j;k þ ΔtðRðŪ n Þ j;k þ SðŪ n ;BÞ j;k Þ ð 4Þ where R = flux residual. For the second-order HW and BH schemes, we use the Heun method, which is a second-order strongstability-preserving Runge-Kutta time integrator (Bouchut 2007;Gottlieb et al. 2001): The Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1967) restricts the time step Δt ¼ t nþ1 − t n and is given by where a and b represent the wave speeds at the interfaces parallel to the x-and y-axes. To ensure the stability of a second-order finite volume scheme, the CFL constant is not allowed to be greater than 0.25 in the two-dimensional case (Bouchut 2007;  (a) the conserved variables U are discretized as cell averagesŪ j;k . The bathymetry function B is discretized at cell centers B j;k for the CN and BH schemes and at cell interface midpoints B j− 1 2 ;k and B jþ 1 2 ;k for the HW scheme; (b) for the BH and HW schemes, left-and right-sided point values are computed at cell interface midpoints using a second-order reconstruction; (c) for the CN scheme, a hydrostatic reconstruction is applied that computes first-order point values from cell averages. For the BH scheme, hydrostatic reconstruction updates the second-order point values from the previous step; and (d) point values are used to compute the fluxes using the HLL flux function at the cell interfaces. Buttinger-Kreuzhuber et al. 2018). For first-order schemes CFL can be increased to 0.5.
We evaluate the friction term in a semi-implicit manner (Brodtkorb et al. 2012), which in the case of a first-order time integration is given by

Schemes
In this section, we discuss the similarities and the differences between the schemes, starting with the CN scheme.
The CN scheme (Chen and Noelle 2017) is based on the hydrostatic reconstruction (HR) scheme of Audusse (Audusse et al. 2004). The original HR (Audusse et al. 2004) uses an upwind evaluation of the adjacent cell-centered bathymetry values for the interface bathymetry values. The first-order HR-based schemes (Audusse et al. 2004;Chen and Noelle 2017) use a piecewiseconstant bathymetry approximation. Thus, interactive bathymetry modifications are easily incorporated into the simulation. In the CN scheme, the bottom interface values need to be recomputed at every time step, which causes a higher workload compared to the HW scheme, e.g., where the bathymetry can be precomputed. It is known that the original HR scheme does not properly account for the acceleration of shallow-water downhill flow (Delestre et al. 2012). This results in incorrect predictions for the water depths and velocities downstream stemming from an unreliable approximation of the bed slope source term for large cell sizes (Morales de Luna et al. 2013). In the HR scheme (Audusse et al. 2004) the source term was designed for near hydrostatic situations, which is also where the name HR comes from. The authors of the CN scheme (Chen and Noelle 2017) present a remedy to this issue and improve the solution by means of a subcell reconstruction in the case of shallow downhill flow. In contrast to the original HR scheme, the authors account also for the acceleration of the flow when reconstructing the bathymetry. This is for the cases when the water level of a downstream cell is lower than the adjacent bottom value. Together with an appropriate source term approximation, comprising both the HR water depths and bathymetry values, the CN scheme is well-balanced and positivity-preserving.
A second-order extension to the original HR scheme exists and is presented by Audusse et al. (2004) and Audusse and Bristeau (2005 (Chen and Noelle 2017) to second-order accuracy. The BH scheme improves the numerical approximation of shallow flows over complex terrain by using a second-order reconstruction tailored to handle abrupt topography changes. To ensure a physically correct solution, care has to be taken when reconstructing the interface bottom values at shoal zones over rapidly changing terrain. Such situations occur in the case of shallow flow over bottom steps, e.g., if the water level on the adjacent cell is lower than the bottom value of the current cell. In this case, the BH scheme reconstructs the interface bathymetry values from the second-order slopes of the bathymetry and the water depth. Otherwise, water level and depth slopes are reconstructed to achieve good balance (Audusse et al. 2004). Moreover, the computational burden was lowered by using a simplified quadrature for the bed slope source term. After the hydrostatic reconstruction, both the CN and the BH schemes delegate the actual work of dealing with wet-dry fronts to the HLL flux, which is capable of handling dry states.
The basis for the HW  scheme is a secondorder scheme developed by Kurganov and Petrova (KP) (Kurganov and Petrova 2007). The former KP scheme is especially suitable for GPU implementation because of its simplicity. However, the KP scheme exhibits some drawbacks. The most crucial one is that it is not well-balanced. Spurious waves can emerge in shoal zones. In long simulation runs it causes the velocities to grow at the wet-dry boundaries, which leads to restricting the time steps toward very small values . Another issue is the drying behavior. If a cell gets wet, it will never get completely dry again. To tackle these issues, the HW scheme introduces a new reconstruction technique inspired by Bollermann et al. (2013), that improves the solution at wet-dry boundaries. In the reconstruction procedure, a dimension-wise separation point is generated based on the intersection of the horizontal water line and the bathymetry. In the partially flooded cells, dimension-wise waterlines are generated and used to reconstruct the water height at the cell interfaces. This technique does not resolve the problem completely. Depending on slope configurations of the bathymetry and the water surface, the scheme might develop temporal instabilities when wetting or drying. However, they do not greatly affect the overall solution and they cease to exist once the cell dries or gets flooded again.
Unfortunately, there is one more drawback that affects not only the KP but also the HW scheme. It is not straightforward to modify the cell-centered bathymetry values for these schemes. The reason for this is the dependency on the vertex values. These values are averaged to compute the bathymetry values at cell centers and interface midpoints. However, the conserved variables, e.g., water depth and discharges, are defined at cell centers. Thus, interactive modification of the bathymetry is highly complicated. It is only possible by adjusting the four vertex values, which affect the bathymetry values in the adjacent cell centers. Instead of applying a local change, one has to modify the surrounding cells as well.
To summarize, we perform a second-order reconstruction of the water levels, water depths, and the velocities using the minmod limiter [ Fig. 1(b)] in the second-order BH and HW schemes. By choosing the minmod limiter we ensure that no high velocity spots can occur inside the domain, which is important for fast and robust flood simulations. In the next step, we apply the hydrostatic reconstruction for the CN and BH scheme. For the CN scheme, the cell averages are used to compute the hydrostatic point values as it is first-order only. For the BH scheme, the already computed secondorder point values are used and updated [ Fig. 1(c)]. All three schemes delegate the actual work of dealing with dry states to the numerical flux, the HLL flux, which is capable of handling wet-dry fronts. Additionally, interface velocities are clamped to zero in the reconstruction step if the water heights are below a defined threshold. In our real-world simulations, we set the threshold to 0.1 mm. For the second-order schemes, we perform the slope reconstruction on the velocities instead of the discharges, which improves the speed estimates at wet-dry fronts. The HW scheme is capable of drying cells by limiting the time step of nearly dry cells, in contrast to the BH and CN schemes, where a thin water layer remains. The more accurate drying feature comes at a price. Because of this, the HW scheme has a more complicated reconstruction and time integration procedure that directly leads to an increased computational demand. The main difference between the second-order schemes is the bathymetry reconstruction at the cell interfaces and the reconstruction for wet-dry fronts. Since the HW scheme requires an averaging of the bathymetry for interface midpoints, sharp features of the terrain, such as dams or dikes, are smoothed [ Figs. 1(a and b)]. In the BH scheme, abrupt changes in the topography are handled better, as a result of the upwind reconstruction of interface bathymetry values and the adaptive slope reconstruction.
From an implementation point of view, the schemes differ in the size of the computational stencil, which describes the number of neighboring cells that have to be accessed when advancing a cell. For CN, the size of the stencil is four cells, i.e., one neighboring cell in each direction, for BH eight cells, and for HW 12 cells. A larger computational stencil requires more data to be processed and leads to slower reconstruction.

Boundary Treatment
For simulating river floods, boundary conditions play an essential role. Typically, for real-world floods, hydrograph data consist of time series of water levels and/or discharge values, which are then prescribed as upstream and downstream boundary conditions (BCs) at the inlet and outlet locations. Usually, hydrograph data are measured at gauging stations and are available as scalar values for the cross-sectional area at the gauge, i.e., total discharge Q g and the water level w g .
For water-level data, the application is straightforward: the same level has to be prescribed to each cell. For discharge hydrographs, we distribute discharge values on the rasterized hydrograph interfaces assuming a constant velocity vectorū ¼ ðū;vÞ T that is normal to the cross-section. For each hydrograph interface I, we prescribe a discharge q I ¼ Q g p I , where p I is the proportional share of this interface with respect to the total discharge Q g . The proportion p I depends on the cross-sectional area of flow at this interface p I ¼ A I =A, where A is the overall cross-sectional area and A I is the contribution area of the cell interface I (Fig. 2). The value of p I is then re-evaluated every time step according to the current water level. To guarantee that the mass influx equals the prescribed discharge, the inflow boundary flux has to satisfy F ð1Þ I ¼ q I in the case of the x-dimension and G ð1Þ I ¼ q I for the y-dimension. The specified boundary flux is based on the solution of a linearized Riemann problem with prescribed interface discharges ðhuÞ I ¼ q I and a continuous extension of the y-discharge in the case of the y-dimension [for details see Pankratz et al. (2007)].
In urban regions, buildings and other structures can completely block water flow. In such cases, so-called wall BCs are used to incorporate these fluid-solid interfaces. For walls and water-level hydrographs, the essential idea is to find a boundary flux FðU I Þ whose interface state U I satisfies certain BCs based on numerical and physical grounds. At a wall boundary, there is no discharge through the boundary interface, so the normal velocity vanishes at the interface: hu I ¼ 0. In a similar way, a boundary flux is deduced for a prescribed water depth h I at each interface given by h I ¼ w g − B I . In our simulations, we use the procedure of Ghidaglia and Pascal (2005) to derive these boundary fluxes for specified water heights and walls.
At the boundary of the computational domain, where the water flow has to be truncated artificially, free outflow BCs are prescribed. Usually, at these free outflows no data are available, hence we use an extension by continuity. These outflow BCs correspond to zero-order absorbing BCs, where the incoming characteristic is set to zero (Ehrhardt 2010).
We remark that the BCs take into account the regime of the flow, e.g., if the state at the boundary is subcritical or supercritical. In other words, we set boundary values only on the characteristics that are entering the computational domain.

Validation
In this section, we present three historical floods and validate the presented schemes. In the simulation, we choose a finer 3 m resolution to accurately capture the topographic features, like weirs and dikes. A coarser 12 m resolution is used to show that in some cases, there are large differences in the simulation results because of the different bathymetry reconstructions. Furthermore, we show how and why the schemes might fail when using a coarser resolution. Our implementation uses GPUs and the CUDA programming model optimized for the NVIDIA Kepler architecture (Horváth et al. 2016). All simulations were performed on NVidia GPUs, namely Geforce GTX 1080. In the following plots, measured data from the events are annotated as ME. Table 1 contains the overall simulation run-times for each case.

Lobau
The first case study involves the Lobau area, which is the alluvial backwater and floodplain of the Donau-Auen National Park in Austria. It extends on the left bank of the Danube River from river kilometer (rkm) 1,918 to rkm 1,908 downstream of the city of Vienna. It consists of floodplain forests and surface water bodies. When the water level in the Danube rises, water flows from the river into the floodplain causing regular flooding events. The Lobau can only be flooded through a small weir, the Schönauer Schlitz [ Fig. 3(a)]. One can observe that the flow through the weir changes its direction when the Danube begins to flood. Even though Lobau is a rural area and does not contain any buildings, this use case is challenging from the simulation perspective, since it has a very complex bathymetry with lots of small channels and steep slopes.
We reconstructed the flooding of January 2011 and simulated the first four days. The size of the simulation domain is approximately 8.1 × 5 km 2 and the resolution is set to 3 m. The upstream and downstream boundaries are located in the Danube. The simulation domain and the initial state along with the inflow, outflow, and gauging locations are shown in Fig. 3(a). The shaded area is marked invalid for simulation. Water level and discharge values for the upstream and downstream boundary conditions are plotted in Fig. 3(b). A nonuniform distribution of the Manning roughness coefficient is used based on the land use [ Fig. 3(c)].  Since all three schemes produce very similar visual images, Fig. 3(d) shows only the maximum water depths computed using the CN scheme. To show differences between the schemes, we picked a smaller region, highlighted in Fig. 3(c) with a rectangular box. This region is enlarged and visualized for each scheme at day 1.5 [ Fig. 3(e)]. In the enlarged images, one can observe that there are temporal differences between schemes, which are also visible in the data recorded at the gauges. The differences are more pronounced when looking at the evolution of the water levels and water extents rather than the visualized water depths. We compare the simulated and the measured water level time series at three gauging locations in Figs. 4(a-c), where day zero on the x-axis corresponds to January 13, 2011, 12 a.m. One may notice that the second-order schemes are able to predict the correct arrival times at the gauging locations PD.LP1, PD.LP16, and PD.LP18. Small discrepancies can be observed between the simulated and measured water levels. Most likely, they are caused by uncertainties in the input data, e.g., bathymetry, roughness, or boundary condition data. The roughness values are not calibrated, and the exact initial state before the flooding is unknown, which may also explain the deviations from the observations. The first-order CN scheme is unable to correctly predict either the water levels or the wave arrival times. This is caused by a faster wave energy dissipation compared to the second order schemes. The wave arrival time is delayed [Figs. 4(b and c) at day 1]. The second-order schemes produce higher water levels at gauge PD. LP16 at day 1.5. This gauge is located near a manually operable weir for which no data is available from the event. This can be responsible for the discrepancies between the simulated and the measured water levels. Overall, the BH scheme produces results that match the measured values the most closely. Fig. 4(d) shows the simulated water extents. At the beginning, the schemes produce similar extents with a minor delay. After 12 h, the water levels produced by the HW scheme start to rise more rapidly and even exceed the CN scheme. An explanation for this is the averaged bathymetry that lowers the elevation in a critical place that would otherwise hold back the water. On the second day, the CN and BH schemes produce almost identical results. However, as the CN is first-order accurate only, it does not correctly resolve the peak. At the end time, water extents produced by the CN and HW schemes are larger compared to the BH scheme. This happens because lower elevation areas get connected and flooded during the peak and they cannot drain after the flood wave retreats.
In order to assess the performance, two indicators were measured. Fig. 5(a) shows the number of time steps needed to compute 1 s simulation time. Fig. 5(b) shows the wall-clock time needed to compute 1 s simulation time. One can observe that the second-order schemes produce almost identical results regarding the time steps and they require approximately twice as many compared to the first-order CN scheme. This difference comes from the CFL condition that halves the time-step sizes for the second-order schemes and thus doubles their number. Fig. 5(b) reveals differences in computation times. First, in the second-order schemes, we have to use second-order time integration which requires two half steps in time doubling the number of computations. Second, the schemes require different computational stencils, i.e., one cell for CN, two cells for BH, and three cells for HW, making the reconstruction procedure more expensive for larger stencils. Discrepancy in water extents between the second-order BH and HW schemes is caused by different bathymetry reconstructions.

Wachau
The second case study focuses on the Wachau valley located in Lower Austria. The valley was carved out by the Danube over thousands of years. It features a riverine landscape with a settled flood plain bounded by steep slopes. This case study aims to reproduce a 100-year Danube river flood (Blöschl et al. 2013). The focus lies on the correct prediction of water levels along the river. This case also emphasizes the proper setting of inflow and outflow conditions at the upstream and downstream boundaries of the simulated reach of the Danube. We reconstructed the flooding of June 2013. Fig. 6(a) shows the initial state and the gauge locations. Measured data are available for four gauging locations along the river, namely, near Stein-Krems (rkm 2002,7), Loiben (rkm 2005,99), Dürnstein (rkm 2009,15), andKienstock (rkm 2015,21). The dataset of Kienstock and Stein-Krems is used to prescribe the upstream and the downstream BCs, respectively. Fig. 6(b) shows the corresponding hydrograph data of the BCs. At the upstream boundary, we use a discharge hydrograph with water level information, and at the downstream we use only water levels. Roughness values were set according to land use data [ Fig. 6(c)].
We simulated 12 days, starting on May 30, 2013 at 5 p.m. with a prefilled river of still water. Simulations were performed using two different resolutions (3 and 12 m). The simulation domain size is 8.1 × 2.6 km 2 . All three schemes produced visually similar results, hence in Fig. 6(d), we show the whole domain only for the CN scheme where the water depths are visualized at 3 m resolution. The rectangular box in the figure marks a region that is enlarged to allow us to spot the differences between the solutions. The results are captured for all schemes at day 3. Fig. 6(e) shows the water depths for the 3 m resolution grid, and Fig. 6(f) shows the water depths for the 12 m grid. As one can see, the biggest difference is exhibited between the first-order and the second-order schemes at 12 m resolution. More water is present in the floodplain on the southern side of the river. The second-order schemes produce visually identical results at both resolutions. Similar differences can be observed also for the 3 m resolution, but they are less significant.
The recorded water levels at the two selected gauges for the finer grid are presented in Figs. 7(a and c), and for the coarser grid in Figs. 7(b and d). For both resolutions, the second-order schemes (BH, HW) produce almost identical results. Both are able to capture the peak accurately within a range of few centimeters. However, discrepancies can be observed before and after peak. Their source can be uncertainties in the hydrograph data, bathymetry, or roughness values. The CN scheme also produces acceptable results for the 3 m resolution where the error is less than 0.5 m on average. However, for the 12 m resolution it overestimates the water levels by about 1 m at both gauges for the whole duration. We also visualize the water extents. However, we have no information regarding the actual size [Figs. 7(e and f)]. By comparing the water extents at different resolutions we see that they reduce in size significantly for the 3 m resolution and there is less discrepancy between the solutions.
We quantify these discrepancies in the water levels by the root mean square error (RMSE) in Figs. 8(a and b). By increasing the resolution, the RMSE of the water levels decreases up to a certain cell size as the schemes better approximate the SWEs and the bathymetry is better resolved. The first-order CN scheme improves throughout the shown resolutions, ranging from 24 to 3 m. Contrarily, the second-order schemes are not obtaining significantly better results below 12 m cell size. At 12 m, they already reach a point where only parameter optimization or more accurate input data would improve the model results.
Besides the water levels, performance and timing are also very important in large-scale high-resolution simulations. For this purpose, we record the number of time steps and the runtime of the simulations. Figs. 9(a and b) show the numbers of time steps needed to compute 1 ssimulation time for the 3 and 12 m resolution domain, respectively. The second-order schemes required approximately twice as many time steps compared to the first-order CN scheme. Furthermore, we remark that, for the 3 m resolution, the number is four times as high for all schemes. For both cases, the reason is the CFL condition. To preserve stability, the firstorder CN scheme uses the CFL number 0.5. For the second-order schemes the conditions are stricter and the CFL number has to be set to 0.25, which in turn halves the time-step sizes. The CFL condition also involves the resolution, which results in four times more time steps when going from 12 to 3 m resolution.
The HW scheme has temporal instabilities, as can be seen in Fig. 9. These instabilities come from high velocity spots that can temporarily limit the time-steps sizes. Even though the HW scheme produces spurious high velocity spots and hence spikes in the plots, it is able to recover from these states and stabilize.

Marchfeld
The Marchfeld is a 900 km 2 -large plain north of the Danube, one of the largest plains in Austria. It is a sedimentary basin located in Lower Austria, between Vienna and Bratislava. The Marchfeld is mainly used for agriculture and is dubbed the granary of Austria. On the east, it is bordered by the March (Morava), which is a border river between Austria and Slovakia, and on the south by the Danube and its floodplains (e.g., Lobau). Since the area is bordered by two rivers, it is prone to natural floods. To protect the inhabitants and villages, a 26 km long dike (Marchfelddamm) was built. This dike serves as the main protection measure against floods [ Fig. 10(a)].
In this case study, we investigate the flood of June 2013 and simulate the whole Marchfeld area including both border rivers, the Danube and the March. The simulation domain size is 17 × 19.5 km 2 . An overview of the area and the gauging locations are shown in Fig. 10(a). The figure also shows the still-water initial state and the protection lines, which are part of the Marchfelddamm. The initial water level is estimated from the BC. The validation and BC data are the values measured during the actual flood at six locations. Specifically, two upstream boundaries are used, namely: Wildungsmauer (Danube) and Marchegg (March). Furthermore, one outflow was prescribed at Wolfsthal (Danube) using water level information [ Fig. 10(a)]. The corresponding hydrograph data is shown in Fig. 10(b). One may notice that discharge data is provided only for one of the upstream boundaries (Wildungsmauer); the other one (Marchegg) only has water level information. We validate our simulations by comparing the results against data from three gauge locations, namely, Bad Deutsch Altenburg, Hainburg, and Thebnerstrassl [ Fig. 10(a)]. Roughness values are estimated from land use information and are not calibrated [ Fig. 1(c)]. Furthermore, the bathymetry is created by merging multiple rasters of different resolutions to complement the main raster. Since all three schemes produce visually very similar images at 3 m resolution, Fig. 10(d) shows only the maximum water depths computed using the BH scheme. To show differences between the schemes, we picked a smaller region that is highlighted in Fig. 10(d) with a rectangular box. This region is enlarged and visualized for each scheme at day 2 [ Fig. 10(e)]. We see that the solutions computed by the second-order schemes, the BH and HW, have a very similar temporal evolution, and the differences are almost unnoticeable. However, the first-order CN scheme produces a larger flood extent at the same time. Fig. 10(f) shows the results computed at 12 m resolution. In this case, the differences appear more pronounced. The dam did not fail during the event, which is also the case for the simulations computed at 3 m resolution. However, it is not the case for the simulations at 12 m resolution. In one case, the CN scheme floods almost the whole Marchfeld and produces the largest flood extent of all schemes. The HW scheme also floods a large area behind the dam, but it is far smaller then in the case of the CN scheme. The BH scheme also floods some area behind the dam, but it is the one closest to the expected results.
Overall, we conclude that with the 12 m resolution grid, none of the schemes could satisfactorily reproduce the event since the dam failed for all of them. Figs. 11(a-f) show the water level time series recorded at the three gauges. We compare the solutions computed at 3 and 12 m resolutions. At 12 m resolution, the second-order schemes produce water levels very close to the ones measured at the gauges. The CN scheme is off by more than a meter on average. When switching to the 3 m resolution the second-order schemes improve a bit and even the CN scheme produces acceptable results.
The CN matches the peaks well at Thebnerstrassl and Hainburg, however it fails to predict the water level at low flow. Furthermore, the simulated water levels are better estimated closer to the outflow boundary, where a measured water level is prescribed. This explains why in Bad Deutsch Altenburg, the simulated water levels differ the most from the observed water levels in this scenario, even for the 3 m resolution.
By inspecting the water level validation plots only, one might conclude that the second-order schemes also perform very well for the 12 m resolution. However, by looking at the water extents in Figs. 11(g and h), a significant difference becomes obvious. This difference is also visible in the visualized water extents in Fig. 10(f), where the dam fails and the water propagates toward the Marchfeld.
For this particular case, the bathymetry resolution is simply insufficient to correctly represent the dam. The reason the CN produces the worst results of all is its first-order accuracy. If comparing the second-order schemes, a significant difference between them can be still noticed. One would assume more similar water extents, since they produce almost identical results at the gauges. However, the HW scheme produces a larger flood extent at 12 m resolution. This is because of the different bathymetry reconstruction of the two schemes. The HW scheme cannot resolve correctly the approximately 20 m-wide dike at 12 m resolution and hence floods the area behind it. This is due to the fact that the bathymetry reconstruction of the HW scheme has an averaging effect on the reconstructed bathymetry so the dam is overtopped.
Finally, we compare the schemes from the performance point of view. Figs. 12(a and b) show the number of time steps needed to simulate 1 s for 3 and 12 m resolutions, respectively. Figs. 12(c and d) show the wall-clock time needed to compute 1 s simulation time for both resolutions. By examining the plots, one can notice that the CN scheme is the fastest for both resolutions and no instabilities are observed. As explained earlier, the CN scheme has to perform half the number of time steps compared to the second-order schemes. The BH and the HW schemes perform similar numbers of time steps. For the HW scheme, local instabilities can be observed at both resolutions. The BH scheme does not produce any spikes in the plot but creates a temporal velocity build-up at the downstream boundary at 3 m resolution around day 7, which vanishes later.

Conclusions
In this paper, we compare three numerical schemes for the shallowwater equations for various study cases. The three schemes (CN, BH, and HW) are well-balanced state-of-the-art schemes, chosen because of their suitability for fast flood simulations on GPUs. The comparison focuses on the prediction of flood wave arrival times and water level series, since they are especially relevant for decision makers. We reconstructed three historical flood events that happened on the Danube and compared the simulation results for different resolutions to reveal the strengths and weaknesses of the presented schemes.
Regarding the choice of the minimal cell size, at least two cells are required to resolve a topographic feature like a dam correctly for the CN and BH schemes. For the HW scheme, because of its averaging effect on the bathymetry, at least four cells are necessary. This also leads to significant differences in water extents at coarser resolutions, when comparing the two second-order schemes. All three schemes improve their solutions on grid refinement up to a certain resolution, where errors from uncertainty in the roughness parameter as well as uncertainty in the input and boundary data start to significantly affect the model results.
Regarding computational effort, the CN scheme is the most lightweight scheme, as it is first-order in space and time. Compared with second-order schemes, a CFL number two times higher can be set, allowing a time step twice as big. Furthermore, no second half-step in the time integration is required. Since the computational stencil is smaller by one cell in each direction, the CN scheme is more than four times faster than the BH scheme at the same cell size. Thus, it enables rough but quick estimations. Looking at the second-order schemes, the BH scheme outperforms the HW scheme as a result of a less complex reconstruction procedure at wet-dry fronts.
The validation case studies show that the second-order BH and HW schemes provide better estimates of the time-dependent water levels than the first-order CN scheme, which overestimates water levels in river scenarios, especially near the discharge-based boundaries. In fact, the water levels of the second-order schemes at the lower resolution (12 m) are closer to the measurements than the water levels of the first-order scheme at the higher resolution (3 m). In terms of run-time, the CN scheme at 3 m cell size is slower than the BH scheme on the 12 m grid.
Thus, considering the trade-off between accuracy and computational effort, the second-order BH scheme is our recommended choice for computing risk maps for river flood simulations. Still, the first-order CN scheme is useful for interactive decision support systems as it offers the speed that is required to run real-time simulations at higher resolutions.

Data Availability Statement
The bathymetry, roughness, and gauge data used in this study are available on request from the via donau-Österreichische Wasserstraßen-Gesellschaft mbH. The source codes of the implemented numerical solvers for the shallow water schemes and the Visdom framework are property of the VRVis Zentrum für Virtual Reality und Visualisierung Forschungs-GmbH.