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

The FVM is chosen for the spatial discretization on top of a uniform grid with cell size $\mathrm{\Delta }x\times \mathrm{\Delta }y$. The starting point for an FVM scheme are cell averages ${\overline{\mathbf{U}}}_{j,k}=[{\overline{h}}_{j,k},{\overline{w}}_{j,k},{\overline{u}}_{j,k},{\overline{v}}_{j,k}{]}^T$ of the conserved variables ${\mathbf{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 ${\mathbf{F}}_{j+\frac{1}{2},k}$ and ${\mathbf{G}}_{j+\frac{1}{2},k}$ at the interfaces and an appropriate source term discretization ${\mathbf{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 ${\mathbf{F}}_{j\pm \frac{1}{2},k}({\mathbf{U}}_{j+\frac{1}{2},k}^{-},{\mathbf{U}}_{j+\frac{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 ${\mathbf{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.”

For the first-order CN scheme, an explicit Euler time integrator is used:where $\mathbf{R}$ = flux residual. For the second-order HW and BH schemes, we use the Heun method, which is a second-order strong-stability-preserving Runge-Kutta time integrator (Bouchut 2007; Gottlieb et al. 2001):where ${\overline{\mathbf{U}}}_{j,k}^{*}={\overline{\mathbf{U}}}_{j,k}^{E,n+1}$ The Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1967) restricts the time step $\mathrm{\Delta }t=t_{n+1}-t_n$ and is given bywhere $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 $\mathbf{C}\mathbf{F}\mathbf{L}$ constant is not allowed to be greater than 0.25 in the two-dimensional case (Bouchut 2007; Buttinger-Kreuzhuber et al. 2018). For first-order schemes $\mathbf{C}\mathbf{F}\mathbf{L}$ 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

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 piecewise-constant 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). Recently, Buttinger et al. (BH) (Buttinger-Kreuzhuber et al. 2018) successfully extended the HR technique of Chen and Noelle (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 (Horváth et al. 2015) scheme is a second-order 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 (Horváth et al. 2015). 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 second-order 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.

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 $\overline{\mathbf{u}}=(\overline{u},\overline{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 ${\mathbf{F}}_I^{(1)}=q_I$ in the case of the $x$-dimension and ${\mathbf{G}}_I^{(1)}=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 $\mathbf{F}({\mathbf{U}}_I)$ whose interface state ${\mathbf{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: $h{u}_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.

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\times 5{\mathrm{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.

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), and Kienstock (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\times 2.6{\mathrm{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 first-order 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.

The Marchfeld is a $900{\mathrm{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\times 19.5{\mathrm{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.