A fast second-order shallow water scheme on two-dimensional structured grids over abrupt topography

This paper presents a finite volume scheme on structured grids to simulate shallow flows over complex terrain. The situation of shallow downhill flow over a step is particularly challenging for most shallow water schemes. We study this situation in detail and devise a novel second-order reconstruction strategy, which gives superior results over former hydrostatic reconstruction (HR) schemes. The reconstruction step is based on a recent first-order hydrostatic reconstruction HR method, which improves shallow flows over steps. The proposed second-order scheme is well-balanced, positivity-preserving, and handles dry cells. When compared with the original HR, we lower the computational burden by using a simplified quadrature for the bed slope source term. We test the scheme on various benchmark setups to assess accuracy and robustness, where the method produces comparable results to other HR-based schemes in most cases and superior results in the case of shallow downhill flow over steps. The novel second-order scheme is capable of simulating large-scale real-world flood scenarios fast and accurately.


Introduction
The shallow water equations (SWEs) describe the motion of an incompressible fluid under the gravitational force. They provide plausible and reliable results of water levels for tsunamis, river floods, dam breaks, and levee breaches ( de la Asunción et al., 2013;Audusse et al., 2004;Brodtkorb et al., 2012;Hervouet and Petitjean, 1999;Liang and Marche, 2009;Russo, 2005 ). 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 equations ( Temam, 1984;Whitham, 1999 ).
The finite volume method (FVM) is a common numerical method for solving the SWEs. For the spatial discretization, a computational grid has to be chosen. Unstructured triangular meshes are able to incorporate complex geometries, however they require time-consuming mesh generation. In contrast, rectangular grids lack the pre-processing step at the expense of poor resolution of topographic features not aligned with the grid. This issue can be overcome with the cut-cell technique ( An and Yu, 2012;Ingram et al., 2003 ). Furthermore, in the context of second-order with an error term coming from the topography ( Audusse et al., 2016 ). Thus, convergence of this scheme can be expected for Lipschitz continuous bathymetry.
In first-order and second-order finite volume (FV) schemes the bottom topography is approximated by piecewise constant and piecewise linear functions, respectively, thus giving rise to discontinuities in the bottom at the discrete level. In practice, large discontinuous bottom steps may occur at coarse spatial resolutions or in stormwater scenarios, where the water layer is often thinner than the stepsize of the bottom jumps. As bottom steps appear at the discrete level of FVMs, there is a need to stretch the applicability of the SWEs also to cases involving discontinuous bathymetries. This leads to mathematical and numerical problems, since then the product of the depth and the bottom gradient cannot be understood in a distributional sense. The mathematical theory of nonconservative products in the source term is an active field of research ( Dal Maso et al., 1995 ). Even if discontinuous bathymetries are not within the theoretical assumptions of the SWEs, Morales de Luna et al. (2013) note that the SWEs still give reasonable results in the case of small enough bottom jumps. Altogether, this motivates the study of shallow water flows over bottom steps and their numerical approximation.
In the case of shallow downhill flow, the original first-order HR scheme does not properly account for the acceleration due to a sloped bottom ( Delestre et al., 2012 ). This effect can be mitigated by switching from a first-order to a second-order approximation. Morales de Luna et al. (2013) improve the original first-order HR scheme in the case of partially wet interfaces. Recently, Chen and Noelle (2017) proposed a new reconstruction which features an even better approximation of the source term in case of shallow downhill flows, leading to a new firstorder scheme, called the CN scheme in this work. Also they present a way to investigate and derive the two existing HR schemes by means of subcell reconstructions. In Xia et al. (2017) , Xia et al. present a Surface Reconstruction Method (SRM) to overcome the problem of partially wet interfaces, which they describe by the term "waterfall effect ". In their first-order scheme, a second-order approximation of the bottom is used in all cells to reconstruct the water surface and the bottoms at the interfaces for the flux and source term computation. Other strategies to further improve the flow over abrupt topography include considering the conservation of the total head instead of the conservation of the hydrostatic equilibrium. Such schemes are also called energybalanced methods and typically require additional waves in the approximate Riemann solver (ARS) to resolve the stationary bottom discontinuity at the interface. This leads to additional complexity in the solver, both implementation-wise and performance-wise ( Goutal et al., 2017;LeFloch and Thanh, 2011;Murillo and García-Navarro, 2010;Murillo and Navas-Montilla, 2016 ).
Another numerical difficulty arises at wet-dry zones, characterized by interfaces between dry and wet cells. A robust numerical scheme should be able to maintain nonnegativity of water depth, but should also avoid unphysically high velocities in these sensitive regions. Typically, at wet-dry zones, the wet cells only feature a thin layer of water thus giving rise to large velocities when the discharge is divided by a small depth. Different methods are tailored to tackle this problem ( Horváth et al., 2015;Hou et al., 2013b ). In Hou et al. (2013b) , a novel source treatment, which is slightly faster than the original source terms as in Audusse et al. (2004) , for unstructured grids is introduced.
In this paper, we present a new two-dimensional scheme, which is second-order accurate. It is based on the hydrostatic reconstruction procedure of Chen and Noelle (2017) . The second-order accuracy allows us to reduce the discretization error and perform more accurate simulation runs. We apply a simple source treatment, which is computationally efficient and leads to a minor loss of accuracy in typical use cases. Furthermore, our proposed reconstruction is adapted to limit the velocities, which reduces the occurrence of unphysically high velocities, and to correctly reconstruct the solution variables in the vicinity of abrupt changes in the bottom topography and in the water levels. In particular, our method is able to capture the drying process in regions with complex terrain in a robust and efficient way. This approach ensures that the time step, which is connected to the velocities by the Courant-Friedrichs-Lewy (CFL) condition, is not overly restricted when simulating large time spans. Our proposed source term approximation coincides with the simple and economical source term of Hou et al. (2013b) in fully wet regions. For shallow flow over abrupt topography, the novel scheme outperforms previous schemes based on the HR method. We remark again that discontinuous bottom steps are not within the theoretical assumptions on the derivation of the shallow water (SW) model, however, they appear by construction in first-order and second-order FVMs. Thus, the correct handling of bottom steps in the numerical approximation, including the reconstruction procedure as well as the source term discretization, is important. We propose to reconstruct the HR water depth from bottom slopes instead of water level slopes in regions where the topography changes abruptly. This novel "adaptive " second-order reconstruction allows us to significantly increase accuracy of shallow flow down a bottom step, when compared to other second-order HR-based schemes.
The paper is organized as follows. In Section 2 , we discuss the model equations, the HR schemes and present the second-order scheme based on a new HR ( Chen and Noelle, 2017 ) in detail. In Section 3 , we present numerical experiments highlighting the advantages and disadvantages of our proposed scheme. We extensively verify the presented scheme on multiple benchmark tests, including a dam break over bottom steps, seven Riemann problems, the parabolic bump, and the parabolic basin. We validate the scheme on the Malpasset dam break and a river flood event. Finally, in Section 4 , we conclude the findings of this work and give a brief outlook into future works. This paper contributes with the following key points: • a two-dimensional well-balanced scheme based on an improved hydrostatic reconstruction • a novel second-order reconstruction which yields superior results for shallow downhill flows over a step, • an economical approximation of the source term to speed up computation, • reduction of unphysically high velocities at wet-dry zones.

The shallow water equations
In this section, we describe the shallow water model and the finite volume method (FVM) for SW schemes. The hyperbolic conservation law described by the two-dimensional shallow water equation (SWE), also referred to as the Saint-Venant system, with geometric source term can be written as where h represents the water height, hu is the discharge along the x -axis, hv is the discharge along the y -axis ( Fig. 1 a), u and v are the average flow velocities in x and y -direction respectively, g is the gravitational constant, and b is the bathymetry (assumed to be time-independent). Subscripts represent partial derivatives, e. g., U t stands for . In vector form the system writes where = [ ℎ, ℎ , ℎ ] is the vector of conserved variables, F and G are flux functions. The bed slope term S models the fluid's acceleration due to the gravitational forces. An additional friction term S f ( U ) can be included on the right hand side of (2) , which is introduced in Section 2.10 . In two dimensions, the SWEs allow for complicated steady state solutions, however we restrict ourselves to two important steady-state equi- libria. Following Chen and Noelle (2017) , there is the still-water equilibrium, i. e., where w denotes the water level = ℎ + , and the lake at rest equilibrium, which includes dry shores, i. e., ℎ , ℎ = 0 and ℎ ∇ = 0 .
If a numerical scheme is capable of balancing source and numerical flux terms for these two stationary solutions it is called well-balanced and thus preserves the lake at rest and also the still-water equilibrium.

Discretization
We choose a uniform grid x ≔ Δx and y ≔ Δy , where Δx and Δy are the cell sizes. We denote by C j,k the cell , ∶= [ −1∕2 , +1∕2 ] × [ −1∕2 , +1∕2 ] . The SWEs are discretized by the method of lines. The FVM is chosen for the spatial discretization on top of the uniform grid. An FVM discretizes the conserved variables U as cell averages, e. g., U j,k for the finite volume C j,k . This yields a system of ordinary differential equations for the cell averages where ∓ 1 2 , and , ∓ 1 2 are the discretized interface fluxes and S j,k is an appropriate source term discretization.

Hydrostatic reconstruction
To achieve well-balancedness, it is necessary to introduce a special reconstruction for the Riemann states that are fed into the approximate Riemann solver. Assuming time-independent bathymetry values b j,k ≔ b ( x j , y k ) at the cell centers, the essential idea of the hydrostatic reconstruction (HR) technique is to redefine the interface bottom values used for deriving the Riemann states in order to ensure wellbalancedness and positivity. The name HR method originates from the fact that the associated HR scheme balances the hydrostatic pressure and the topographic source terms at each interface in the still water steadystate. We illustrate the technique only in x -direction, the application to the y -direction can be done analogously.
In the following, we briefly summarize the original first-order HR of Audusse (Aud) ( Audusse et al., 2004 ) and the modification of Chen and Noelle (CN) ( Chen and Noelle, 2017 ). In both first-order schemes, there is only one hydrostatically reconstructed bathymetry value at each interface. The original first-order HR evaluates the interface bottom values b * in an upwind fashion ( Audusse et al., 2004 ), * , Aud ) .
The left-and right-sided interface heights with respect to the interface +1∕2 , between cell C j,k and +1 , are then defined as ℎ * , Aud , This definition ensures the nonnegativity of the water depths h . The CN scheme ( Chen and Noelle, 2017 ) improves the original reconstruction on partially-wet cases where an adjacent water level is lower than the bottom topography, i. e. +1∕2 , ( ) ∶= min The CN scheme defines the interface bottom values b * as * , CN The interface heights h * are given by For a constant water level W , the interface depths are then continuous across interfaces, that is max ( − * , Aud . The continuous depths then make it easy to show well-balancedness for consistent fluxes, see Section 2.9 . We remark that the hydrostatic interface heights h * of the two first-order schemes do not differ, only the HR interface bathymetry values differ in partially wet cells ( Chen and Noelle, 2017 ), compare also Fig. 2 .

Second-order reconstruction
For second-order accuracy, left-and right-sided point values have to be computed at the cell interface midpoints through slopes taking into account the neighbouring values. The left-sided point values of a cell C j,k are then denoted by subscripts − 1∕2+ and − 1∕2+ in x -and ydimension, respectively. We reconstruct = [ ℎ, , , ] instead of the conserved variables ( Audusse et al., 2004;Bouchut, 2007 ). From now on, we omit time dependence in the equations, since all reconstructed variables are time-dependent. To suppress unphysical oscillations, the generalised minmod-limiter is applied to the slopes ( Nessyahu and Tadmor, 1990;Sweby, 1984;Van Leer, 1979 ) where between 1 and 2. The minmod limiter is given by Nessyahu and Tadmor (1990) The value of controls the amount of dispersion added to the system and is chosen to be 1.3 in our simulations. This choice of slope limiting conserves the maximum principle. The left-and right-sided water depth point values are given by and the left-and right-sided water level point values by

Velocity reconstruction
According to Bouchut (2007) , we reconstruct velocities [ , ] = [ ℎ ∕ ℎ, ℎ ∕ ℎ ] instead of discharges hu, hv to avoid high velocities near dry cells. To satisfy the conservativity requirement dimension-wise on the discharges, the point values have to satisfy and and analogously for hv . Thus, the velocity point values are modified accordingly to these equations giving and , − 1 The same modification is applied for the velocity v . We remark that the depth point values in (14) are not yet hydrostatically reconstructed.
If the interface water depth, e. g., ℎ +1∕2− , , is smaller than some dry threshold dry , we set the respective velocity to zero in the reconstruction step, i. e. hu j,k / h j,k is set to zero if h j,k < for all j, k .

Second-order hydrostatic reconstruction
Here, we describe the procedure used in Audusse et al. (2004) to get second-order HR water heights h * at the interfaces. The left-and right-sided second-order bottom point values are given by subtracting the water depth from the level, i. e.
From the second-order bottom point values at the interfaces, the HR interface bottom values are set to * , Aud Then the hydrostatic interface heights h * are reconstructed by .
The procedure is visualized in Fig. 3 a and b.

"Adaptive " second-order hydrostatic reconstruction
We introduce a new second-order reconstruction in combination with the recent hydrostatic reconstruction introduced by Chen and Noelle (2017) . We propose a reconstruction that is additionally based on bottom values in case of "large " discontinuities in bathymetry and water levels, as explained for the x -dimension in the following paragraph. If we might land in a "partially wet " situation, i. e.
then, at this cell, we check if the bottom slope is greater than the water level slope D x w j, k . If the conditions and (23) hold, we reset the water level slope We proceed by recalculating the second-order water level point values (15) . Condition (25) ensures a correct treatment at partially wet cells and at wet-dry regions of an advancing wave front on a sloped we additionally reconstruct the bottom slope if interfaces are partially wet. d) In this case, the water level slope is recomputed from the water depth and bottom slope and used to derive the HR bottom values (red asterisks). We see that the left HR interface depth ℎ * −1∕2+ differs in the two second-order HR reconstructions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) bottom. We remark that the bottom slopes, e. g., D x b j,k in (24) , can be precomputed for time-independent bathymetry values.
Then, we derive a second-order HR bottom value per interface, * , CN from the water level point values (15) and Audusse's HR bottom values b * ,Aud (21) , which are also depending on the water level point values.
We set the second-order HR left and right interface depth values h * to using the water depth (14) and level (15) 27) and (28) are the second-order analogs to Eqs. (9) and (10) . The HR reconstructed depth values (28) and (22)  This adaptive reconstruction strategy is necessary, since a naive second-order reconstruction can not be applied to all cells as in the vicinity of strong bottom jumps a back wave might emerge at the top of the step. This unphysical behavior is caused by an unphysical bottom reconstruction in the upper cell of the bottom jump.

Source terms
We revisit the bathymetry source term S in a cell C j,k for the secondorder HR scheme of Audusse et al. (2004) . Clearly, the depth source term S h is zero, i. e. ℎ , = 0 .
The momentum source terms S hu and S hv are split into interface parts and a centered part The centered source terms read .
ensure second-order consistency in regions where the solution is smooth. If the solution is varying a lot, the source term is distributed towards the interfaces ( Bouchut, 2007 . The smoother the solution, the smaller the differences of the variables across the interfaces. Thus the source term is mostly influenced by the cell centered source term. However, since the source terms (35) only depend on the smoothness of the water level w and depth h , discontinuities in the bottoms might be interpreted only as a centered source term reflecting the water level slope ( Fig. 4 ). We approximate the source term without a centered source term ( Hou et al., 2013a;2013b ) where the left and right interface source term is given by and analogously for ℎ , . The simple approximation slightly reduces the computational burden and is easier to implement. The difference between the two source term approximations is visualized in Fig. 5 for fully wet cells. We remark that our source term treatment consisting of (36) -(38) also leads to differences in wet regions, when compared to the original HR method ( Audusse et al., 2004 ). A similar approximation of the source term without the centered part can also be found in the works of Hou et al. (2013a,b) with the following source terms ℎ (36) and (36) Hou ,ℎ In the robust and simple scheme of Hou et al. (2013a,b) only the bed elevations at the lower side are modified resulting in different bottom values at the left and the right interface to maintain well-balancedness, i. e. * , Hou Thus, there is no acceleration in the upper cell's source terms coming from the bed slope. Thus, as in the scheme of Audusse et al. (2004) , it does not fully account for bathymetry steps at shallow flow conditions. If in two adjacent cells there is shallow flow and we are in a partially wet situation (23) , then the simple source term approximation of Hou et al. neglects the difference +1∕2 , − * , Aud +1∕2 , . Thus, the contribution of the bottom jump to the source term at the upper cell is neglected, in contrast to our approach. This is further highlighted in the numerical experiments, see Section 3.1 and 3.3 .
Analogously, speed values in y -direction can be derived from the Jacobian of the flux in y -dimension G . We remark that the HLL flux is consistent, i. e.  HLL ( , ) = ( ) and with the given choices for the speeds it is also able to handle dry states. The HLL flux is remarkable robust, however, it is known that the HLL flux does not resolve shear waves accurately as it ignores the contact discontinuity of the transverse velocity ( Toro, 2001 ). One way to fix this issue is to include the middle wave, which leads to the Harten-Lax-van Leer-contact (HLLC) flux. To preserve the nonnegativity of the water depths and to satisfy a discrete entropy inequality, additional subcharacteristic conditions have to be satisfied ( Bouchut, 2007 ). It is possible to approximate the HLLC flux with a simplified version, i. e., Other choices for the approximate Riemann solver (ARS) include the Roe solver. Although the Roe solver is more accurate than the HLL flux, it has difficulties with dry beds. The handling of dry beds is often incorporated by imposing internal boundary conditions, which adds complexity to schemes using the Roe solver ( Castro et al., 2005;Murillo and García-Navarro, 2012;Murillo and Navas-Montilla, 2016;Parés and Pimentel, 2019 ). Moreover, at transcritical rarefactions, i. e., if the left or right eigenvalue is close to zero, an entropy fix is needed ( Harten and Hyman, 1983;LeVeque, 1992 ). Thus, the Roe solver is a priori less robust and requires an additional parameter, which decides if speeds are close to zero and thus considered for a transcritical rarefaction fix ( Toro, 2001 ).

Time integration
For first-order time integration, an explicit Euler is used, i. e. where Quantities denoted by a superscript n depend on the state U n ≈ U ( t n ). The Courant-Friedrichs-Lewy (CFL) condition restricts the time step Δ = +1 − and is given by where and represent the maximum wave speeds in x -and ydirection at time t n . They are computed by a reduction over all interface wave speeds given by the absolute values of (42) and (43) . The CFL constant has to be positive and is not allowed to be greater 0.25 to ensure the positivity of the two-dimensional second-order accurate finite volume (FV) scheme, as we show in Section 2.8 .
Heun's method is used for second-order time integration. By denoting intermediate states with an asterisk, the state +1 at time +1 is given by where the residual R ( U ) is defined according to (46) . Clearly, by (50) , the intermediate state * , +2 does not need to be explicitly calculated. The solution U is updated dimension-wise. First, we compute reconstructed values, fluxes and sources in x -dimension. Second, we perform the computations in y -dimension. Afterwards we update the solution with the combined residual.

Positivity preserving
Our scheme preserves the nonnegativity of the water depths, i. e.
The new second-order scheme guarantees nonnegative water depth for the SWEs (1) under condition (47) with the CFL number being halved when compared to the CFL condition needed for a positivitypreserving first-order scheme associated with the homogeneous problem, e. g., the HLL scheme. This statement can be proved similarly as in Audusse et al. (2004) , Audusse and Bristeau (2005) and Bouchut (2004) . The positivity follows from the facts that the reconstructed depths are nonnegative (22) and that the chosen second-order time integration is a convex combination of two first-order time steps.
The two-dimensional scheme is positivity-preserving under half the CFL condition needed for the positivity-preserving one-dimensional scheme. Thus, by choosing a positive CFL constant not greater than 0.25 we obtain a positivity-preserving two-dimensional second-order scheme. We remark that the different source treatment does not have any influence on the preservation of nonnegative states. Furthermore, the numerical speed of the HR schemes is not higher than the one of the associated homogenous scheme, since the numerical speed is a monotone function of the water depth h , see (42) and (43) .

Well-Balancedness
We show well-balancedness in two steps. First, for the still-water steady state and, second, for the lake at rest steady state involving wetdry boundaries. We remark that since our proposed scheme does not couple dimensions for the flux and source terms, it is enough to show well-balancedness dimension-wise. In fact, each interface can be associated with a subcell for which we will show well-balancedness. We will use the following convex decomposition of the residuum , = − 1 2 + , + + 1 2 − , + , − 1 2 + + , where − 1 with analogous definitions for the other subcells.
In the still water situation, we have ∇ = 0 and ℎ , ℎ = 0 and we analyse the nontrivial case h > 0. In this case, all cells are fully wet and the reconstructed interface bottom levels and water depths agree with the ones defined by Audusse et al. (2004) . For a cell C j,k , we will show that the left residuum − 1 2 + , (58) vanishes. For a constant water level W , which is preserved by our second-order reconstruction, i. e.
by consistency of the HLL flux and since the water is at rest. Furthermore, we notice that as the water level W is greater than the bathymetry, and use it to compute the left source term Together, this shows The residuum vanishes also for all other subcells by the same reasons, yielding a vanishing cell residuum , = . Thus, the new twodimensional second-order scheme is well-balanced for the still-water steady state (3) .
In the lake at rest situation, we have ℎ ∇ = 0 and ℎ , ℎ = 0 . It is enough to show well-balancedness for wet-dry boundaries in one dimension, e. g., a dry-wet front ( Fig. 6 ), i. e. ℎ , = 0 , ℎ +1 , > 0 and , > +1 , . For a dry-wet front , = , > +1 , = +1 , + ℎ +1 , > +1 , holds. We have to show that the residuum in these cells vanishes, that is , = 0 and +1 , = 0 . At a wet-dry front, the bathymetry slope D x b can only be greater or equal the water level slope D x w , thus condition (25) is not true, and we only need to consider water levels and depths as second-order reconstructed variables. The Therefore, the residual R j,k equals zero, since the cell depth h j,k and interface depths ℎ * −1∕2+ , and ℎ * +1∕2− , are zero. For cell +1 , , we will use the convex decomposition of the residuum (57) since the water is at rest. As ℎ * +1∕2+ , = 0 , we have With ℎ +1∕2 , = 0 we conclude that ℎ +1∕2+ , = 0 . For a completely wet right interface +3∕2 , ) , holds, compare the proof of the previous theorem. However, if +1 , is degenerate, i. e., cell +2 , is dry, the second-order scheme falls back to first order, in which case well-balancedness follows from the first-order CN scheme. In this case, Eq. (72) also holds ( Chen and Noelle, 2017 ). Thus, the residua ℎ , in x -dimension vanish for all j, k . An adaptation of the previous arguments shows that the residua ℎ , in y -dimension vanish. By inspection of the flux terms, ℎ − 1 2 + , = ℎ + 1 2 − , = 0 for all interfaces in the lake-at-rest state Together, this shows that our novel scheme is well-balanced, also for the lake at rest steady state (4) . We remark that our scheme is well-balanced on a per-interface basis, thus this property holds also on unstructured grids if the second-order reconstruction keeps the water levels balanced.

Friction source terms
To provide realistic water flow, a friction term is introduced in the laboratory and real-world scenarios. The friction term S f is included via an additional source term where n is the Manning roughness coefficient. It is evaluated in a semiimplicit manner by splitting the friction source term S f into a coefficientwise product of an implicitly evaluated state and an explicitly evaluated friction term ̃ ( Brodtkorb et al., 2012 ) where Then, the integration from time t n to +1 including friction is achieved by using the following explicit update of the states * , +1 = + Δ ( ) instead of (48) and (50) .

Validation
We validate the scheme on various test cases, a dam break over a step, the parabolic bump, seven Riemann problems, the parabolic basin ( Thacker, 1981 ), the Malpasset dam break event, and a historical flooding. Additionally, we verify the order of the scheme at the parabolic basin and the parabolic bump. In the following sections, we denote the scheme of Audusse et al. (2004) by Aud and the proposed second-order scheme by BHNW. The implementation of the scheme of Audusse et al. only differs in the hydrostatic reconstruction (HR), the adaptive secondorder reconstruction and in the source term approximation. In particular, the generalised minmod slope limiter and the velocity reconstruction with a dry threshold was used for all schemes. The gravity constant g equals 9.81 in all our simulations, except for the parabolic basin where it is set to 2. The dry threshold dry is set to 10 −6 in the dam break and in the parabolic basin and to 10 −4 in the Malpasset and Lobau. The experimental order of convergence EOC is defined as where U is the exact solution and U N is the numerical approximation on a mesh with cell size Δx · Δy , while U 2 N is the numerical approximation on a mesh with half of the cell size, i. e. Δx /2 · Δy /2. We use either the discrete L 1 -norm of the water depth, or the maximum water depth difference, that is, the L ∞ -norm.

Dam break over a dry step
We describe the setup for a dam break over a dry step, as specified in Chen and Noelle (2017) based on numerical experiments from Bollermann et al. (2013) and Castro et al. (2008) . The difficulty lies in the correct approximation of the wet/dry front and the bottom step.
As noted in the introduction, discontinuities in the bottom are outside the validity range of the shallow water (SW) model. However, bottom steps necessarily occur at the discrete level in finite volume methods (FVMs) and thus motivate this test. The quasi one-dimensional test is performed on a domain with range [0, 1] × [0, 0.01].
The bottom topography b and the initial water depth h 0 is given by and respectively. We use a uniform cell size of 0.0025 m for the simulated values. The reference solution is computed on a grid with a cell size of 10 −5 m and a piecewise linear step at = 0 . 1 with 100 cells in the transition layer. The transition layer used for approximating the bottom step is thus 0.001 m wide. We display the results at a final time = 0 . 18 s in Fig. 7 . In this case, the scheme of Audusse et al. and the robust scheme of Hou et al. (2013a) produce nearly identical results. As already noted in Section 2.5 , these two schemes neglect the jump in the water levels at the interface, which leads to incorrect predictions of the velocities after the step. The improved HR method of the CN scheme together with the novel adaptive second-order reconstruction enables us to capture the water flow after the step accurately.

Parabolic bump
This section is devoted to show the performance of the scheme on a quasi one-dimensional steady-state test with a parabolic bump. The scenario is set up analogously to Audusse and Bristeau (2005) and Delestre et al. (2013) and is originally from Goutal and Maurel (1997) . The analytical solutions for the steady states can be derived using the Bernoulli relation, see Bouchut (2007) and Delestre et al. (2013) . The bathymetry is given by for a domain of length L = 20 m and a width of 4 m ( Fig. 8 ).
In the case of subcritical flow, a discharge boundary condition (BC) is specified at the inflow = 0 and a water level BC at the outflow = .
The water depth is given by where ℎ = = 2 m is the water depth at the outflow boundary. The discharge in x -direction is specified as = 4 . 42 m 2 /s at the left inflow boundary. Water levels and velocities are shown in Fig. 8 for a cell size of 0.5 m. Since in this setup all cells are always flooded, the CN HR falls back to the original HR and thus the first-order Audusse scheme agrees with the CN scheme. The BHNW scheme with simple source term produces results nearly identical to the ones of the second-order schemes of Audusse et al. (Aud RK2). This is also visible in Table 1 , showing that the simple source treatment has only a very small effect on the accuracy.    Due to the increased diffusivity in first-order schemes, the water levels at the discharge inflow are underestimated (right zoom-in Fig. 8 ). An error analysis for a range of cell sizes starting from 1 m down to 0.125 m shows that the schemes are second-order accurate in the smooth subcritical flow regime ( Fig. 9 ).

Riemann problems
We test the scheme on several Riemann problems (RPs) including resonant cases. All RPs are defined on the domain [−1 , 1] and by an initial state U 0 consisting of a left state = ( ℎ , ) with a left bottom level b L for x < 0, and a right state = ( ℎ , ) with a right bottom level b R for x > 0. The exact analytical solution is given by completing the SW system (1) with = 0 and connecting the resulting Riemann states. This extended inhomogenenous system shows a rich solution pattern. In fact, the RP may have no, a unique, or multiple solutions, depending on the given states and the bottom jump ( Han and Warnecke, 2014;LeFloch and Thanh, 2007;2011 ). We restrict ourselves to cases with a unique solution. The investigated setups are listed in Table 2 . All of them result in a flow from left to right at the bottom jump. The analytical solution is computed as outlined in Han and Warnecke (2014) . All simulations are run until 0.1 s. The cell size is set to 0.002 m. We plot the water level and Froude number to emphasize the criticality of the flow states. In the plots, the gray area represents the bottom topography and the initial water level is marked with a thin dashed line. The x-axis limits are adapted to the RPs.
Riemann problem 1 ( Fig. 10 ) is a dam break over a bottom jump. The solution consists of a left rarefaction wave, a stationary shock associated with the bottom jump and a right shock. As all interfaces are fully wet, the combination of adaptive reconstruction and different source treatment does not have any visible effect. Therefore, the second-order Audusse scheme produces almost exactly the same results as the BHNW scheme, both are in good agreement with the analytical solution.  Riemann problem 2 ( Fig. 11 ) is a two shock case over a bottom jump. In this case, all interface are fully wet after 0.002 s and therefore the second-order Audusse scheme again produces almost exactly the same results as the BHNW scheme. However both schemes fail to accurately predict the state at the top of the bottom jump. These first two Riemann problems can also be found in García-Navarro (2010, 2013) .
In Riemann problem 3 ( Fig. 12 ), we test a supercritical regime over a downward bottom step. The state at the right of the bottom jump is also supercritical, it is not accurately captured by both the Audusse and the BHNW scheme. In Riemann problem 4 ( Fig. 13 ), we test a resonant regime over a small downward bottom step. In the resonance regime, the emerging solution pattern is quite complex and involves critical intermediate states or transcritical waves. As the stationary shock associated with a bottom jump is not allowed to cross the boundaries of strict hyperbolicity, the left-most wave has to be a rarefaction wave from the subcritical left state to a critical state. Then, this critical state connects via a stationary shock to the supercritical states at the right. Again, the supercritical state at the right of the bottom jump is not accurately captured by both schemes. RP 3 and RP 4 are taken from LeFloch and Thanh (2011) .
In Riemann problem 5 ( Fig. 14 ), we test a dam break over a medium downward bottom step. The left-most wave is a rarefaction from the subcritical left state to a critical state, which then goes into a stationary hydraulic jump at the bottom discontinuity. In fact, the analytical solution shows that at the bottom discontinuity = 0 three waves are present. First, a stationary shock shifting the bottom level = 1 down to an intermediate bottom level = 0 . 5 accompanied by a supercritical intermediate state. Then, a stationary hydraulic jump causes the supercritical intermediate state to become subcritical, which is then followed by another stationary shock that shifts the bottom level down to = 0 . The left rarefaction wave is not fully captured by both schemes, instead a wrong intermediate state emerges that connects the left state with the subcritical state at the right of the bottom jump. This is an artefact of the second-order reconstruction as the first-oder CN scheme is able to capture the critical state, see Fig. 14 . We remark that simply using the classical minmod-limiter, i. e. = 1 , is not enough to recover the correct intermediate state.
In Riemann problem 6 ( Fig. 15 ), we again test a resonant regime over a downward bottom step, but this time with the right water level below the left bottom elevation. The emerging solution pattern is similar to Riemann problem 3. Again, the first left wave is a rarefaction from the subcritical left state to a critical state, which then connects via two more shocks to the right subcritical state. Since at the bottom jump, the interface is partially wet, the adaptive reconstruction enables the   BHNW scheme to capture the rarefaction wave at the left. To make this point clearer, we also compare the BHNW without adaptive reconstruction against the Audusse scheme with adaptive reconstruction in Fig. 16 . In fact, none of these two variants is able to capture the left rarefaction wave. This demonstrates the necessity for the adaptive second-order reconstruction. Moreover, the BHNW scheme provides a better estimate of the right state of the bottom discontinuity, when compared to the scheme of Audusse et al.
Riemann problem 7 ( Fig. 17 ) is a resonant regime over a large downward bottom step, that connects a subcritical left state with a supercritical right state. The emerging solution pattern is similar to Riemann problem 6, except that the right-most wave is now a rarefaction wave.   Here, also the scheme of Audusse et al. gets the left rarefaction right, however the BHNW achieves superior predictions of the right supercritical states.
Concluding this section of Riemann problems, we observe that the new BHNW scheme is able to outperform the scheme of Audusse et al. in the partially wet cases, while never performing worse than it.

Thacker's planar solution
Thacker's planar solution, sometimes also referenced as the parabolic basin, is a classical test case for validation. Thacker (1981) provides an analytical solution. It describes time-dependent oscillations of a planar water surface in a parabolic basin. It is widely used for compar- Fig. 16. Riemann problem 6. We compare two variants, the BHNW scheme without adaptive reconstruction and the Audusse scheme with adaptive reconstruction. Both fail to match the full wave pattern. Fig. 17. Riemann problem 7. Partially-wet resonant regime over large downward step connecting a subcritical left state with a supercritical right state. The BHNW and Audusse's scheme capture the full wave pattern, however the BHNW provides better estimates of the right supercritical states. Fig. 18. Parabolic basin, well-balancedness. An initial lake at rest with constant water level = 0 m for a parabolic basin is simulated for 1390 s on a grid with a cell size of 160 m. At the cross-section = 0 , which is marked with a red line, we extract the velocities of wet cells. The magnitude of the velocities is within single-precision floating-point accuracy for our proposed BHNW second-order scheme, thus it numerically preserves the lake at rest steady-state. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)  ing different numerical schemes ( de la Asunción et al., 2013;Gallardo et al., 2007;Horváth et al., 2015;Liang and Marche, 2009;Sampson et al., 2006 ). Recently, Sampson et al. (2006) extended the solution of Thacker to support bed friction. However, their solution is limited to one dimension. In this two-dimensional case, we use the same setup as Holdahl et al. (1999) , where the bathymetry is given by where = 2500 m , 0 = 1 m . First, we use a constant water level = 0 m to show wellbalancedness ( Fig. 18 ). We can see that the velocity errors are within the accuracy of single floating-point numbers, which we used in our implementation.

Malpasset dam break
The Malpasset dam in southern France collapsed in 1959, resulting in a 40 m high water wave flooding the Reyran valley. The event was exhaustively studied in recent years ( Brodtkorb et al., 2012;George, 2011;Hou et al., 2013a;2013b;Singh et al., 2011;Valiani et al., 2002 ). We investigate the dam break on a structured grid with a cell size of 20 m. Friction is included with a uniform roughness coefficient of = 0 . 033 m 1/3 /s, corresponding to weedy, stony earth channels and floodplains with pasture and farmland. We compare simulation results with laboratory experiments of a 1:400 scaled model ( Frazao et al., 1999;Hervouet and Petitjean, 1999 ). In this experiment, arrival times of the wave front ( Frazao et al., 1999 ) and maximum water levels ( Hervouet and Petitjean, 1999 ) were recorded at 14 gauge locations, labeled S1-S14 in     No data are available for the first 5 gauges, thus we use gauge locations S6-S14 in our validation ( Fig. 22 a and 22 b). Additional data is also available for the shut-down time of voltage transformers of the historical event ( Fig. 22 c). The locations of the gauges and transformers, as well as the inundated area after 2000 s is displayed in Fig. 21 . Small discrepancies between the scale model and the numerical results were also reported in other studies ( Brodtkorb et al., 2012;George, 2011;Hou et al., 2014 ), and our results are consistent with these. When compared with the scheme of Audusse et al. (2004) , minor differences only occur for gauges S8-S10 for the water levels and for gauges S11 and S14 for the wave arrival times, with the new BHNW scheme obtaining comparable or slightly better results for most of the gauges except for the maximum water level at gauge S8. Regarding performance, our scheme is slightly faster, 0.3% run time reduction for the first 2000 simulated seconds, than the second-order scheme of Audusse et al. on a parallel implementation running on a machine with a 4-core Intel i5-4960K CPU at 3.5 GHz. Our proposed scheme performs better with increased simulation time because of the different reconstruction and source term treatment in regions with small water depths and complex terrain ( Fig. 22 d).

Lobau
The Lobau is a floodplain east of Vienna, in Lower Austria, located at the left bank of the Danube. It consists mostly of floodplain forests and is regularly flooded. We simulate a flood that occured in January 2011 with a CUDA GPU implementation on a NVIDIA GeForce GTX 1070. The initial time is set to 13 January 2011, 1am, and the initial state comprises several still-water bodies and the Danube ( Fig. 23 a). The water is flowing from the Danube into the Lobau only through a small slot, the "Schönauer Schlitz ". The terrain is quite complex, featuring several small channels, which render simulations challenging.
We apply an inflow BC upstream of Fischamend and an outflow BC downstream of Fischamend ( Fig. 24 a). The inflow BC is applied as a discharge BC as in Pankratz et al. (2007) , and the otuflow BC is implemented as a flux boundary condition based on water levels ( Dutykh et al., 2011;Ghidaglia and Pascal, 2005 ). The simulation domain is roughly 10 × 7 km 2 large and the simulation cell size is set to 4 by 4 m. The bathymetry is given on a raster with 2 m resolution. The Manning roughness coefficient n varies spatially between 0.03 and 0.13 m 1/3 /s. It is estimated based on the land use. Hourly measured water levels are available at three locations, PD.LP1, PD.LP16 and PD.LP18. They are displayed alongside the simulated data in Fig. 24 b. The simulated water extent after 2.5 days is displayed in Fig. 23 b. The exact initial state is unknown and there is also an operated weir at the Gänsehaufentraverse, which might explain the small differences between the observed and the simulated flood waves at the gauge locations as it was modelled as a constant bathymetry modification. Taking into account those uncertainties, the measured water levels are predicted very well by the simulation.

Conclusion and future work
We derive and test a new formally second-order finite volume method (FVM) scheme for the shallow water equations. The scheme is well-balanced, as it preserves both the still water and lake at rest steady states, and does not exhibit any oscillations. Instead of reconstructing the discharge slopes, we reconstruct the velocity slopes to obtain robust choices of the wave speeds at wet-dry fronts, ensuring fast simulations. The scheme is particularly suited for implementations on graphics processing units (GPUs), thus enabling faster than real-time simulations for large domains.
A numerical convergence analysis demonstrates that the scheme is second-order accurate. Validation against several benchmark tests, including multiple Riemann problems, reveals that the scheme converges against the reference solutions in most cases. Still, there are some scenarios where the solver does not produce satisfactory results. The scheme is able to reproduce real-world flood events such as the Malpasset dam break and a historical river flood in Austria. On test cases with shallow flow over abrupt topography, the new scheme achieves superior results than existing schemes. These improvements are due to an improved hydrostatic reconstruction (HR) procedure and a novel adaptive second-order reconstruction strategy, which enables accurate resolution of shallow flow down a bottom step. Moreover, our proposed scheme only requires modification of a few lines of code when compared to the HR scheme of Audusse et al. To sum up, the scheme is able to capture complex flows over complex terrains accurately and efficiently as shown in the numerical test cases, all the more in the presence of thin water layers.
The scheme can be applied to unstructured grids, as the source terms are evaluated on a subcell basis, only the slope reconstruction needs to be revisited. In real world cases, the friction term plays an important role in predicting the correct evolution of the flood extent. We are planning to improve our scheme by balancing moving water in the presence of friction to gain better estimations of water levels and wave arrival times. Further work is directed to combine the scheme with an infiltration model for rainfall-runoff simulations.
Agency in the scope of COMET -Competence Centers for Excellent Technologies (854174) which is managed by FFG.