1 Introduction

Finite volume methods, initially developed for solving problems in gas dynamics, have been accepted as a reliable and accurate tool for the numerical solution of the shallow water equations [20]. From pioneering works to present, the solution of the Euler equations for gas dynamics is based on the physics and the mathematical properties of the hyperbolic differential equations. Hyperbolic equations provide solutions using a set of elementary waves derived from the equations themselves. Numerical techniques were first derived on that basis for the inviscid Euler equations and have been considered a starting point in the numerical modelling of shallow water flows traditionally paying little attention to the presence of source terms associated to bed friction and variations of the bed slope. Unsteady problems were tackled with explicit schemes where the stability was controlled by the Courant–Friedrichs–Lewy (CFL) condition. Nevertheless, when solving real problems all sorts of situations are likely to be encountered and, in a wide range of them, flow is dominated by bed friction and bed variations.

Numerical experience has shown that appropriate discretization of source terms is necessary. It can be argued that in some cases even naive discretizations of the source terms may work, but there are well documented situations in which only sophisticated schemes can perform adequately. When incorporating the presence of the source terms in a given specific finite volume scheme, the main focus has been traditionally put on the preservation of still water equilibrium. This requirement has led to the notion of well-balanced schemes that have resulted in efficient explicit finite volume models of shallow water flow. But, when well-balanced numerical schemes are applied to real cases with irregular geometries and transient flow, involving large variations in the bed slope, advance over dry areas and also the drying of initially wetted regions, these schemes can compromise the quality and reliability of the solution [23].

In complex cases, flow features impose heavier restrictions than the classical Courant–Friedrichs–Lewy (CFL) condition on the time step size, derived for the homogeneous case without source terms, that may lead to inefficient computations that cannot be explained when coming back to the fundamentals of the basic numerical scheme. One of the most dramatic consequences is the appearance of negative values of water depth, not only in wet/dry fronts, but also in initially wet areas [13].

Therefore, even though well-balanced schemes have proved successful in a large variety of cases, they still fail in some situations. In order to avoid instabilities and a dramatic reduction of the time step, it is usually recommended the use of different tuning parameters in the numerical simulation, as limiting values of bed resistance or minimum cut-off depths values that compromise mass conservation and have to be calibrated in each case. As a result, the numerical results can be strongly affected by the selection of these values. It is remarkable that the lack of conservation may invalidate the numerical results.

The presence of source terms requires the construction of new numerical schemes appropriate to the nature of the equations, instead of using extensions of schemes constructed for the simple, homogeneous case. It is important to remark that, in presence of source terms, the number of waves involved in the equations is different from the number of waves provided in absence of them. Considering that numerical schemes constructed departing from the homogeneous case are based on the definition of this reduced number of waves, it is clear that numerical schemes for hyperbolic equations with source terms deserve their own analysis and development.

The definition of augmented numerical schemes (augmented in the sense of considering an extra wave arising from the presence of source terms) allows explaining all problems detailed before and provides the correct steps to estimate correctly the different types of source terms involved in the problem [14]. This means that no tuning parameter neither any flux redistribution is necessary. Augmented numerical schemes allow recovering the stability region given by the basic CFL condition and no time reduction is necessary. Exact conservation is guaranteed. It can be said that the renewed descriptions made for augmented numerical schemes are the key to explain and avoid important difficulties in real applications. These advances, in combination with the geometric flexibility inherent to unstructured-grid models, provide an excellent tool in environmental engineering.

Physically based simulations of complex systems usually require large computational facilities to be completed in a reasonable time. Moreover when the simulated phenomenon is unsteady and based on a dynamical estimation of the updating time step, the computational performance is an important topic to be taken into account. One of the most widespread strategies to reduce the computational cost is the use of parallel techniques, involving a suitable number of processors. Geophysical flows usually involve large size domains and long time scales. Practical applications require a compromise between spatial accuracy and computational efficiency. In order to achieve the necessary spatial resolution, rather fine grids become necessary in many cases requiring more data storage, increasing proportionally the number of operations and reducing the allowable time step size for explicit calculations. When, at the same time, a reasonable computational time is desired, the use of GPU codes is one of the options for computing large space and temporal domain problems. The idea of accelerating the calculations in unsteady hydraulic simulation using parallel implementations and using GPU was recently reported in [3, 9,10,11, 18, 22].

2 Mathematical model

2.1 Governing equations

The mathematical model considered is based on the shallow flow equations, where the general three-dimensional conservation laws are depth averaged. The pressure distribution is considered hydrostatic and several frictional terms can be assumed so that the 2D hydrodynamic equations are written in global coordinates as follows:

$$\begin{aligned} \frac{\partial {\mathbf {U}}}{\partial t}+\frac{\partial {\mathbf {F(U)}}}{\partial x}+\frac{\partial {\mathbf {G(U)}}}{\partial y}= {\mathbf {S}}_{\tau }+ {\mathbf {S}}_{b} \end{aligned}$$
(1)

where

$$\begin{aligned} {\mathbf {U}}=\left( h,hu,hv\right) ^{T} \end{aligned}$$
(2)

are the conserved variables with h representing flow depth in the z coordinate and (uv) the depth averaged components of the velocity vector. The fluxes are given by

$$\begin{aligned} {\mathbf {F}}=\left( hu,hu^{2} + \frac{1}{2} g_{ \psi } h^{2} , huv \right) ^{T} \nonumber \\ {\mathbf {G}}=\left( hv, huv, hv^{2}+ \frac{1}{2} g_{ \psi } h^{2} \right) ^{T} \end{aligned}$$
(3)

with \(g_{ \psi } = g \cos ^{2}{ \psi }\) and \(\psi\) the direction cosine of the bed normal with respect to the vertical. The physical basis of this gravity projection is explained in [7] and it is of utmost importance for not ruining the numerical predictions when the simulation involves the presence of steep slopes.

The term \({\mathbf {S}}_{\tau }\) denotes the frictional effects in the bed, and is defined as

$$\begin{aligned} {\mathbf {S}}_{\tau }=\left( 0 , - \frac{ \tau _{b,x} }{\rho } , - \frac{ \tau _{b,y} }{\rho } \right) ^{T} \end{aligned}$$
(4)

with \(\tau _{b,x},\tau _{b,y}\) the bed shear stress in the x and y directions respectively and \({\rho }\) the density. The main rheological properties are governed by the frictional forces. These interactions are computed by means of semiempirical friction laws depending on the case. One of the most widely used friction formulations is Manning’s law where:

$$\begin{aligned} \tau _{bx}= \rho g h \frac{n^{2} u \sqrt{u^{2}+v^{2}}}{h^{4/3}}, \quad \tau _{by}= \rho g h \frac{n^{2} v \sqrt{u^{2}+v^{2}}}{h^{4/3}} \end{aligned}$$
(5)

in terms of the Manning roughness coefficient n.

On the other hand, the term \({\mathbf {S}}_{b}\) is the term including forces associated to the bed slopes

Thanks to the hyperbolic character of (1) it is possible to obtain a Jacobian matrix, \({\mathbf {J_{n}}}\), which is built by means of the flux normal to a direction given by the unit normal vector \({\mathbf {n}}\), \({\mathbf {E_n}}={\mathbf {F}}n_{x}+{\mathbf {G}}n_{y}\),

$$\begin{aligned} {\mathbf {J_{n}}}= \frac{\partial {\mathbf {E_n}}}{\partial {\mathbf {U}}} =\frac{\partial {\mathbf {F}}}{\partial {\mathbf {U}}}n_x + \frac{\partial {\mathbf {G}}}{\partial {\mathbf {U}}}n_y \end{aligned}$$
(6)

whose components are

$$\begin{aligned} {\mathbf {J}}_n=\left( \begin{array}{lll} 0 &\quad n_x &\quad n_y \\ (g_{\psi } h-u^{2})n_x-uv n_y &\quad v n_y+2u n_x &\quad u n_y \\ (g_{\psi } h -v^{2})n_y-uv n_x &\quad v n_x &\quad u n_x+2v n_y \end{array} \right) \end{aligned}$$
(7)

The eigenvalues of this Jacobian matrix constitute the basis of the upwind technique which is detailed in the next subsection.

On the other hand, the bed evolution can be modeled through the Exner equation, which is basically a movable bed continuity equation where the bed level time variations are due to the solid fluxes which cross the control volume. In this work the authors only focus on highly concentrated bed-load phenomena and, consequently, the 2D Exner equation is:

$$\begin{aligned} \frac{\partial z}{\partial t} + \xi \frac{\partial q_{s,x}}{\partial x} + \xi \frac{\partial q_{s,y}}{\partial y} = 0 \end{aligned}$$
(8)

where \(\xi = \frac{1}{1-p}\), p is the material porosity and \(q_{s,x}\), \(q_{s,y}\) are the solid discharges in the x and y directions respectively. They are computed as a function of excess bed shear stress with respect to the critical value and taking into account the bed shear stress direction. This bed load transport is often expressed through the following dimensionless parameter [4]:

$$\begin{aligned} \varPhi = \frac{|{\mathbf {q}}_{s}|}{\sqrt{g_{ \psi }(s-1)d_m^{3}}} \end{aligned}$$
(9)

where \(s=\rho _s/ \rho\) is the ratio of solid material (\(\rho _s\)) over water (\(\rho\)) densities, and \(d_m\) is the grain median diameter. According to the numerical assessment performed in [7], several empirical formulae can be chosen for computing the dimensionless bed load discharge. Among them, the present work includes results using the Meyer-Peter and Muller [12] formula:

$$\begin{aligned} \varPhi =8( \theta -\theta _c)^{3/2} \end{aligned}$$
(10)

where \(\theta\) is the dimensionless shear stress

$$\begin{aligned} \theta =\frac{n^{2}}{(s-1)d_m h^{1/3}}\left( u^{2}+v^{2}\right) \end{aligned}$$
(11)

and \(\theta _c\) is the critical shear stress. This formula is only applied when the shear stress is larger than the critical shear stress. Otherwise there is no sediment transport.

2.2 Numerical method

The numerical scheme is constructed by defining an approximate Jacobian matrix \(\widetilde{{\mathbf {J}}}\) at each k edge of length \(l_k\) between neighbouring cells defined through the normal flux \({\mathbf {E_n}}\) so that the volume integral in the cell i of area \(A_i\) at time \(t^{n+1}=t^{n}+\varDelta t\) is expressed as:

$$\begin{aligned} {\mathbf {U}}^{n+1}_{i} = {\mathbf {U}}^{n}_{i} - \sum _{k=1}^{NE} \sum _{m=1}^{3} ( \widetilde{\lambda }^- \widetilde{\alpha } - \widetilde{\beta }^- )^m_k {\mathbf {\widetilde{e}}} ^m_{k} l_k \frac{\varDelta t }{A_i} \end{aligned}$$
(12)

The superscript minus in (12) implies that only the incoming waves are considered for updating the flow variables of each cell, defining \(\widetilde{\lambda }^-=\frac{1}{2}(\widetilde{\lambda } -|\widetilde{\lambda }|)\) where \(\widetilde{\lambda }\) are the discrete representation of the Jacobian matrix eigenvalues and \({\mathbf {\widetilde{e}}}\) the corresponding eigenvectors. Coefficients \(\alpha\) and \(\beta\) are used to properly formulate the Riemann solver and can be found in [14]. Furthermore, special care is considered when calculating wet/dry fronts with the present method. The strategy proposed is based on enforcing positive values of interface discrete water depths coming from a detailed study of the Riemann problem [14, 15]. When they become negative, the numerical values of the friction and bed slope source terms is reduced instead of diminishing the time step.

Equation (8) is also integrated in a grid cell \(\varOmega _i\). Using Gauss theorem:

$$\begin{aligned} \frac{d}{dt} \int _{\varOmega _i} z {d}\varOmega + \oint _{\partial \varOmega _i }\xi q_{s{\mathbf {n}}} {dl} = 0 \end{aligned}$$
(13)

where \({q}_{s{\mathbf {n}}}=(q_{s,x}n_x+q_{s,y}n_y)\).

Assuming a piecewise representation of the variable z and that the second integral can be written as the sum of fluxes across the cell edges, the bed level is updated as:

$$\begin{aligned} z^{n+1}_{i} = z^{n}_{i} -\sum _{k=1}^{NE}\xi q_{s{\mathbf {n}},k}^* \frac{\varDelta t \; l_k}{A_i} \end{aligned}$$
(14)

where:

$$\begin{aligned} q_{s{\mathbf {n}},k}^* = \left\{ \begin{array}{ll} q_{s{\mathbf {n}},i} &\quad {{\rm if}} \;\; \widetilde{{\mathbf {\lambda }}}_{s}>0 \\ q_{s{\mathbf {n}},j} &\quad {{\rm if}}\;\; \widetilde{{\mathbf {\lambda }}}_{s} <0 \\ \end{array} \right. \end{aligned}$$
(15)

where \(q_{s{\mathbf {n}},i}\) and \(q_{s{\mathbf {n}},j}\) are the bed load discharge computed at the neighbouring cells i, j, and \({\mathbf {\widetilde{\lambda }}}_{s}\) is the numerical bed celerity estimated as:

$$\begin{aligned} \widetilde{{\mathbf {\lambda }}}_{s} = \frac{\delta q_{s{\mathbf {n}},k}}{\delta z_k} \end{aligned}$$
(16)

with \(\delta q_{s{\mathbf {n}},k}=q_{s{\mathbf {n}},j}-q_{s{\mathbf {n}},i}\) and \(\delta z_{k}=z_{j}-z_{i}\) provided that the bed is not flat. Otherwise, the numerical bed celerity is estimated as the average liquid flow celerity.

2.3 Numerical stability and time step calculation

The explicitly updated conserved variables are defined through the fluxes obtained within each cell, so, the computational time step has to be chosen small enough for ensuring a stability region. Traditionally, the numerical stability has been controlled through a dimensionless parameter, CFL,

$$\begin{aligned} \varDelta t = \text {CFL} \; \frac{ \min (\chi ) }{ \max _{k,m} |\widetilde{\lambda }^{m} |_k }\quad \text {CFL} \le 1 \end{aligned}$$
(17)

where \(\chi\) is a relevant distance between neighbouring cells [14] and \(\widetilde{\lambda }^{m}\) are the hydrodynamic celerities. The stability criterion is revisited for including a discrete estimation of the bed celerity, \(\widetilde{\lambda }_{s}\), as in [8],

$$\begin{aligned} \varDelta t = \text {CFL} \; \frac{ \min (\chi ) }{ \max _{k,m} |\widetilde{\lambda }^{m},\widetilde{\lambda }_{s} |_k }\quad \text {CFL} \le 1 \end{aligned}$$
(18)

With this numerical strategy, the stability condition takes into consideration the most restrictive numerical wave speed coming from the hydrodynamical and morphodynamical solvers. The resulting global time step is used for updating the whole set of conserved hydrodynamic and morphological variables in the system of equations.

3 Numerical results

3.1 Case 1. Unsteady flow around a conical island

In this section, the model is applied to a laboratory test case originally performed by [2]. A tsunami wave, generated by a lab-scale wavemaker, passes around a conical island, located in the center of a \(26\,{\mathrm{m}} \times 27.6\,{\mathrm{m}}\) flat and frictionless domain. The base diameter of the island is 7.2 m and the top diameter 2.2 m with a height of 0.625 m. A still water level of \(H_0=0.32\) m is assumed as initial condition and the inlet wave is set as a boundary condition for water level H in one of the 27.6 m sides as follows:

$$\begin{aligned} H=H_0+A\,sech^2\left[ \frac{C_1(t-T)}{C_2}\right] \end{aligned}$$
(19)

where

$$\begin{aligned} C_1=\sqrt{g H_0}\left( 1+\frac{A}{2 H_0}\right) ,\quad C_2=H_0\sqrt{\frac{4 H_0 C_1}{3 A\sqrt{g H_0}}} \end{aligned}$$
(20)

being A the solitary wave amplitude and T the time at which the wave peak enters the domain (\(A=0.032\) m and \(T=2.84\) s for the case considered in this work).

The purpose of this test is to prove the robustness of the 2D shallow water model when dealing with transient and long-wave runup over a 2D topography with abrupt changes in the dry–wet interfaces. Experimental data were measured at three gauging points [2] located near the island as shown in Fig. 1. Figure 2 shows a 3D plot of the water surface evolution at different times. The numerical results are compared with observed data in Fig. 3.

Fig. 1
figure 1

Case 1. Water level probes location

Fig. 2
figure 2

Case 1. 3D representation of the water level at \(t=0\), 9 s, 10 s, 15.5 s, 25 s and 30 s (z axis is exaggerated 5\(\times\))

Fig. 3
figure 3

Case 1. Water level (numerical vs. observed) for all the gauging points

Figure 2 shows a 3D representation of the water level all over the domain at several times. As seen, stability and symmetry is fully preserved at all times of the simulation. Figure 3 shows the comparison between numerical and observed data for the three water level gauges. In general, the shallow water model accurately predicts the main water level peaks and reflections, specially for gauges 2 and 3.

3.2 Case 2. Rainfall-runoff simulation in urban areas: triangular grids and local mesh refinement

This test case has been selected from [21] (Copyright: Environment Agency, 2009) that discusses the theoretical background to 2D flood inundation modelling and makes recommendations for benchmark test cases to differentiate between 2D model types in terms of performance and predictive capability. The work is devoted to definition of the benchmark test cases, the application of a range of software packages to these tests and comparative reporting of the outcome from the tests.

This example has been chosen to test the model capability to simulate shallow inundation originating from a point source and from rainfall applied directly to the model grid, at relatively high resolution. This section explores the use of two different kind of triangular grids combined with the possibility of locally refining the mesh. The test considered here, is a rectangular domain of \(960\,{\mathrm{m}}\times 400\,{\mathrm{m}}\) in the city of Glasgow (UK). A spatially uniform rainfall event of 400 mm/h of intensity between \(t= 1\) min and \(t= 4\) min is assumed all over the domain. It is combined with an isolated source inflow located in the northeastern part as displayed in Fig. 4. The local inflow discharge is assumed in a Gaussian shape centered at \(t=37\) min of maximum discharge \(Q_{max}=5\,{\mathrm{m}}^{3}{/}{\mathrm{s}}\) and an standard deviation of \(\sigma =7.5\) min. These elements are included as mass source/sink terms in the continuity equation. Free open boundaries are considered. Two Manning’s roughness coefficients are used, distinguishing between road and pavements \(n_1 = 0.02\) and the rest of the soil uses \(n_2 = 0.05\). Some gauge points are placed inside the domain, shown in Fig. 4.

Fig. 4
figure 4

Case 2. Domain of study and gauge location

Despite the wide use of structured meshes, complex geometries for internal or external boundaries are problematic to be represented if not using unstructured meshes. Moreover, when dealing with topographic representation some recent works [5] have shown the benefit of using unstructured meshes in unsteady hydraulic simulations over irregular topography. The numerical results is sensitive to the grid resolution. Hence grid refinement is clearly an option to modify the whole resolution. In that sense, adaptive grid refinement, readily available when using triangular unstructured meshes [19] and designed to follow local bed variations or irregular boundaries can be very useful. Four different mesh configurations are examined in this section: a uniform Delaunay grid, a locally refined Delaunay grid, a uniform equilateral grid and a locally refined equilateral grid. The detail of the four kind of meshes is displayed in Fig. 5.

Fig. 5
figure 5

Case 2. Detail of the meshes: Uniform Delaunay (upper left), uniform equilateral (lower left), locally refined Delaunay (upper right) and locally refined equilateral (lower right)

The results in terms of temporal evolution of water surface elevation (left) and velocity magnitude (right) at a few probes for all the grids are plotted in Fig. 6.

Fig. 6
figure 6

Case 2. Temporal evolution of water surface elevation (left) and velocity magnitude (right) at gauge points 1–5

The results are sensitive to the mesh used. This is noticeable in the water surface level but even more in the water velocity values. The computational times in this case are very similar on all the grids. To compute a 5 h event, the model took around 1 h when using a single Intel Core i7-3770K CPU 3.5GHz 30 min when using a 4 processor parallel CPU version and 1min when using a GPU implementation with GPU NVIDIA TITAN Black [16].

3.3 Case 3. Large scale inundation flow

The Ebro River basin is managed by the Ebro River Basin Authority (Confederación Hidrográfica del Ebro-CHE) (www.chebro.es), an autonomous institution depending of the Environmental Ministry which is responsible for the management, regulation and maintenance of the drainage basin of this river. Measurements registered in their meteorological and gauging stations have helped us to carry out this study. It is a river with an average discharge of \(Q = 400\,{\mathrm{m}}^3{/}{\mathrm{s}}\) reaching often \(Q = 2500\,{\mathrm{m}}^{3}{/}{\mathrm{s}}\) (Return period 5 years). In this work, we have focused our study in the middle part of this river which is the most affected zone by the floods. It is a 125 km long stretch of river that flows between the towns of Castejón de Ebro (upstream) and Zaragoza (downstream) over a total extension of \(744\,{\mathrm{km}}^{2}\). The information provided at those locations is used as boundary conditions in the numerical model.

To carry out an optimal computational representation of the domain, different features were taken into consideration: the roughness, the topography and the mesh refinement. To represent the roughness, different soil uses were considered; to introduce topographical information a \(5\,{\mathrm{m}}\times 5\,{\mathrm{m}}\) digital terrain model was used together with a river bed reconstruction [6] and a calibration process to find important zones requiring a high level of discretization. Dry bed initial conditions (\(h=0, u=v=0\)) were assumed to generate suitable steady flow conditions that could be used as realistic initial conditions for inundation flow simulations.

The footprint of the maximum flooded area and the time evolution of the depth of water in some measurement stations between Castejón de Ebro and Zaragoza (Tudela and Novillas) were also used to calibrate the mesh. A formula [1] which provides the % of adjustment between the observed maximum flood area (\(A_{Obs}\)) and the corresponding calculated area (\(A_{Sim}\)) as the ratio between its coincidence and their sum is used

$$\begin{aligned} Fit_A(\%)=100\frac{A_{Obs}\cap A_{Sim}}{A_{Obs}\cup A_{Sim}} \end{aligned}$$
(21)

This parameter penalizes the percentage in which areas do not coincide; either by excess or default. Analysing this parameter, we have also to consider that in \(A_{Obs}\) the information comes from the limits of the flooded area without considering dry land inside if it exists. For this reason, rather than considering \(A_{Sim}\), the boundary convex hull of \(A_{Sim}\) has been used.

For the 2015 event, a steady state of discharge \(Q=473\,{\mathrm{m}}^{3}{/}{\mathrm{s}}\) was used as initial condition. The upstream boundary condition corresponding to the 2015 flood was a discharge hydrograph reaching a discharge peak of \(Q_{max}=2691\,{\mathrm{m}}^{3}{/}{\mathrm{s}}\) in Castejón (see Fig. 8). The duration of the event is almost 21 days, one of the largest of the events analysed, and it has been computed in 17 h (with GPU NVIDIA TITAN Black). A snapshot plan view of the flooded surface along the river reach can be observed in Fig. 7 where a zoom view of the velocity field near Alagón is also provided to highlight the detail of information provided by the model. A summary of the data used in the two flooding events reported here is provided in Table 1.

Fig. 7
figure 7

Case 3. Ebro inundation plan view plot with a detail of the velocity field near Alagón

Table 1 Case 3. Large scale inundation events

Since Ebro River basin authorities manage interventions when a flood event takes place, their interest on discharge time evolution and hydrograph conveyance makes its representation also quite interesting. In Fig. 8 the inlet hydrographs as well as the time evolution of discharge at strategic points where gauging stations are placed is shown for two analysed cases, the 2009 (left) and the 2015 (right) event. The evolution of the hydrograph along the river is in good agreement with the observed data. The correct simulation of the arrival time of the flood wave at each point is quite important, and the model is able to reproduce it properly. The other relevant effect is the change on the wave shape due to the storage capacity of the floodplain, which also can be seen on the results due to the use of a 2D model guaranteed by the quality in the mesh construction.

Fig. 8
figure 8

Case 3. Inlet hydrograph, measured and predicted hydrographs for the 2009 Ebro flood (left) and the 2015 Ebro flood (right) at gauge points of Tudela and Zaragoza

In Fig. 9 the time evolution of the measured and calculated water surface elevation at Tudela (upper), Novillas (middle) and Alagón (lower) are represented for both the 2009 (left) and the 2015 (right) event. The prediction of the water surface levels is also satisfactory. Due to the set of initial conditions, the results differ at the beginning of the simulation. The tendency followed by the numerical results is quite similar to the measured data, and only differences of cm are seen in terms of water level. The importance of having complete, updated and accurate topographical information so that the computational mesh results are optimal cannot be ignored. River beds are constantly changing due to erosion and other phenomena and actual topography information is necessary. The simulations are both based on recent topographical data and the results corresponding to the oldest event (2009) display less agreement. An optimal computational mesh has been generated and calibrated for this reach of the Ebro River where the 2015 event has been reproduced reaching a flooded surface fit of \(84. 3\%\). All topics previously appointed to generate an accurate computational model need high performing technology to manage with large domains (\(744\,{\mathrm{km}}^{2}\)), great number of computational cells (867,672) and fast and accurate resolution. All of that is obtained using a GPU-parallelized upwind numerical scheme which simulates a hydrograph of 21 days in 17 h and makes feasible to reproduce events on a real-time basis.

Fig. 9
figure 9

Case 3. Comparison of measured and predicted water surface elevation at three points along the river for the 2009 event (left) and 2015 event (right)

3.4 Case 4. Unsteady erosive flow

This section aims at validating the erosion/sedimentation component of the model. Two dam-break test cases have been simulated and compared with the experimental data provided in [17], with the Meyer-Peter and Müller erosion model (see Table 2). The bed sediment is assumed as 100% gravel with a mean diameter \(d_m=3.0\) mm, a density \(\rho =2390\,{\mathrm{kg}}\,{\mathrm{m}}^{-3}\) and a porosity \(p=0.41\). The roughness is modelled by setting a Manning’s coefficient \(n=0.01\) for the channel structure and 0.018 for the sediment.

The case setup is a 80 m long, 1.2 m wide and 0.8 m deep glass flume, as shown in Fig. 10. The position of the probes is also shown in the same sketch. The reservoir has a fixed bed with two different longitudinal slope values (\(S_0=0\) and \(S_0=0.003\)). Downstream the erodible bed has a constant longitudinal slope of \(S_0=0.003\). A rigid wall 0.31 m high at the end of the domain allows to hold a second water level. For all the simulations, a 147, 000 cell locally refined unstructured triangular mesh is used for the spatial discretization of the domain.

Fig. 10
figure 10

(Adapted from [17])

Case 4. Sketch of the flume for the erosional test case

Table 2 Case 4. Simulations setup

Figure 11 show the water level results for all the probes located along the whole domain, together with the experimental results from [17] for Test 2 as a sample. The numerical water levels agree with the measured values in every single probe of the domain for the four test cases. Regarding the erosion numerical results, a general good behaviour is observed for the Meyer-Peter and Müller model.

Fig. 11
figure 11

Case 4, Test 2. Water level (numerical vs. observed) for all the gauging points

4 Conclusions

A complete 2D dynamic shallow water model in combination with the Exner equation, based on finite volumes over triangular grids has been used to simulate environmental flows including urban flooding. It offers the possibility to fit the calculation domain to complex geometries. The discretization method proposed is fully conservative and able to tackle highly irregular boundaries. The validation of the model as a tool has been based on water depth and water level from laboratory test cases and a field test case corresponding to documented flood events in the Ebro river. In the present paper, this technique has been shown to work well in very demanding hydraulic situations of practical interest.

The main conclusions that can be drawn from our numerical study is that water depth predictions against measured values are good indications that the mathematical formulation, based on the main hypothesis of hydrostatic pressure and Manning friction model, is a good compromise between description of physical phenomena and practical interest.

From the numerical point of view, the extended Riemann solver used is useful to ensure robust solutions controlled by the classical CFL condition and, therefore, avoids the necessity of reducing the stability region, refining the grid or using tuning parameters. It is able to cope with unsteady flow over rough, irregular and dry beds and produces stable and conservative numerical solutions. Furthermore, the accuracy provided by the first order finite volume scheme is sufficient for the kind of flow and the uncertainty of the experimental and field data available.

It must be pointed out that the accuracy of the numerical solution is related to the grid refinement but the proposed numerical model is able to provide accurate results of the water depth variable on triangular grid cells with nearly the size of the main streets. The good performance of the numerical model proposed in the situations of interest presented offers the possibility to simulate different scenarios and configurations for predictive purposes. The use of unstructured meshes has been considered, since this grid topology is the only one which avoid misleading preferential flow directions and offers an easy way to achieve local mesh refinement.

For this purpose, the shallow water equations have been discretized in Finite Volumes and the numerical schemes implemented to run on a GPU card. The resulting model can be used over realistic bathymetries in affordable computation time even when considering large domains and retaining a high level of accuracy. The saving of computational time allows to address large-number-of-cells, large-time and large-space scenarios, strengthening preventive measures and enhancing response capacities. This opens the possibility of facing the hydraulic and sediment transport analysis in a particular location for several years or the geomorphological changes in domains of a regional-size.