Numerical simulation of 2D real large scale ﬂoods on GPU: the Ebro River

. Modern ﬂood risk management and mitigation plans incorporate the presence of numerical models that are able to assess the response of the system and to help in the decision-making processes. The shallow water system of equations (SWE) is widely used to model free surface ﬂow evolution in river ﬂooding. Although 1D models are usually adopted when simulating long rivers due to their computational e ﬃ ciency, 2D models approximate better the behaviour in ﬂoodplains of meandering rivers using a ﬁne mesh which implies una ﬀ ordable computations in real-world applications. However, the advances on parallelization methods accelerate computation mak-ing 2D models competitive. In particular, GPU technology o ﬀ ers important speed-ups which allow fast simulations of large scale scenarios. In this work, an example of the scope of this technology is presented. Several past ﬂood events have been modelled using GPU. The physical domain (mid-dle part of the Ebro River in Spain) has a extent of 477 km 2 , which gives rise to a large computational grid. The steps followed to carry out the numerical simulation are detailed, as well as the comparison between numerical results and observed ﬂooded areas reaching coincidences up to 87.25 % and speed en-hancements of 1-h of simulation time for 1-day ﬂood event. These results lead to the feasible application of this numerical model in real-time simulation tools with accurate and fast predictions useful for ﬂood management.


Introduction
A UN Survey [24] reveals the severity of flooding episodes regarding their consequences, not only in terms of material damages but also concerning human losses. Prediction and prevention procedures are planned in order to mitigate these effects and are essential to restraint the limits of those losses. With this aim, simulation can be a competitive prediction tool if numerical models are able to provide accurate results at affordable computational time.
The simulation of flooding events is a common practice for companies and river basin administrations that nowadays demand more accurate and faster models involving larger temporal and spatial scales. One-dimensional (1D) hydraulic models are the most used for these purposes [1,7,15,16,18,21], mainly due to their efficiency. Nevertheless, when combined with real-time gauging stations to build a Decision-support Service, there may be a loss of accuracy in their predictions due to their limitation in simulating the floodplain-river interactions. On the other hand, two-dimensional (2D) hydraulic models have been accepted as a way to overcome the 1D constraints. These models are able to simulate river flood-plain and complex geometries involving wet/dry boundaries. However, as they require a great level of detail in the form of number of cells to ensure a reliable accuracy, their disadvantage resides in the necessity to refine the computational grid and the large amount of computations, which turns into unaffordable simulations for real-time systems or large space and time domains.
However, computational development is enhancing some methods and some alternatives have emerged recently. In particular, the GPU technology, which uses the Graphic Processing Unit for computations, is based on the parallelization of the calculation and thus its acceleration. The CUDA-language implementation (needed for the use of GPU) turns 2D computation into an affordable technique for large space and time domains due to the computational time reduction [4,11,20].
The physical domain (477 km 2 ) includes a part of the Ebro river (125 km long), with a drop from 295 m to 119 m. This low slope, together with the nature of the materials going through, is responsible for one of the main peculiarities of the river, as is its pronounced meandering. It is supplied with an Automatic Hydrological Information System (SAIH) providing meteorological and water level data at gauging points all over the Ebro basin. The Hydrographic Confederation of the Ebro (CHE) uses this system to predict on real time the main flood variables and alert agencies of Civil Protection. Last flooding events in the Ebro basin have raised several questions about the quality of the information provided by CHE. Hence, CHE needs more efficient tools that help minimizing damage and financial losses occurring in flooding events.
In the present study, the reach is represented by a 2D mesh in which Shallow Water Equations for free surface flow are solved and numerical results corresponding to the simulation of the 2015 flood event of more than 2000 m 3 /s maximum discharge are shown and compared with field data. The computational time required demonstrates that this type of models can be used in real time prediction systems.

Governing equations and numerical scheme
In a compact form, the 2D shallow water system of equations that governs free surface flow can be written as in [14]: where U is the vector of conserved variables U = (h, hu, hv), h is the water depth, (u, v) are the components of the depth-averaged velocity vector in the (x, y) plane, t is the time variable, E = (F, G) is the flow vector of conserved variables in the form being g the gravity acceleration. On the right side of (1) the source terms include the bed slope and the friction term: The bed and friction slope terms in x and y directions are formulated as follows where z b is the bottom depth and n is the Manning's coefficient, which depends on the terrain and comes from experimental analysis [22,23]. In every cell in which the physical domain is discretized, the system of equations (1) is solved using a finite volume numerical method where Ω refers to the finite volume (cell), dΩ is the cell boundary and n denotes the unit normal vector to each cell edge, k. The Jacobian J n of the normal flux (E · n) is diagonalized in terms of a matrix formed by its eigenvaluesλ m k . Then, Roe's linearization is used to decouple the original hyperbolic system and to define an approximate matrixJ nk at each cell interface k [11,13].
The change on the vector of conserved variables and the source term vector across each cell edge, k, can be expressed in terms of the eigenvectors,ẽ m [12], of J n as: Therefore, the conserved variable vector can be updated at cell i as follows: where ∆t is the time step size, A i is the cell surface area, l k is the length of each computational cell edge andγ represents the flux and source terms contributions from the adjacent cells [12]. Finally, NE is the number of edges in the cell (NE = 3 in a triangular unstructured mesh). The explicit nature of the numerical scheme requires the time step ∆t to be restricted by the CFL condition [13] in order to ensure stability where The CFL number remains constant during the whole simulation in order to guarantee the numerical stability. CFL=0.9 has been used for all the cases.
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. The idea of accelerating the calculations in unsteady hydraulic simulation using multiple CPU was reported in [10,19] or [3] as well as using GPU in [2,9,17] or [4]. Although a very good compromise between number of CPUs used and performance is offered by the former option, the cost of using multiple CPU is significant due to the hardware investment and associated use. Alternatively, the GPU technology offers the performance of smaller clusters with less disbursement [8]. The main difficulty, and apparent drawback, when porting codes from CPU to GPU, is the cell order required by the GPU to process data efficiently. This drawback is not present when dealing with structured meshes due to the inherent order and a simple and efficient implementation is relatively easy to be obtained. On the other hand, CUDA is a programming model designed by NVIDIA to run parallel solutions on their GPUs. In this work, the GPU implementation of the numerical method has been made using CUDA. In this work, GPU implementations, already developed and detailed in previous papers, are tested in the simulation of real size flooding events comparing their results with field data, showing the possibilities that this technique can offer to the users of similar models.

Case study
The selected case of study is a reach in the Ebro River basin and encompasses a total area of more than 470 km 2 . Starting at Castejón de Ebro (Navarra, Spain) and ending at Zaragoza (Aragón, Spain), the whole domain covers 125 km of the river, as seen in Figure 1, where a few locations have been highlighted, some of them with gauging stations. For a proper computational representation of the terrain a Digital Terrain Model (DTM) of 5 m resolution and a bed roughness distribution have been used to characterize the river channel and all the hydraulic structures present. Since terrain level is measured by means of laser technology, the water surface interferes on the measures and the river bed level is not properly represented in DTM's, as seen in Figure 2 (a), where elevation, z b , is illustrated. Thus, the river channel has been represented by the use of cross sections, as seen in Figure  2 (b). And finally, a land use map has been used to provide information of roughness spatial distribution, as seen in Figure 2 (c), where n is displayed. Different soils origin have been taken into account by means of different Manning's roughness coefficients.
Recent works [5] have shown the benefit of using unstructured meshes in unsteady hydraulic simulations over irregular topography. That is the reason an unstructured 678000 cells (10 m long in the river and 150 m in the boundaries) has been used in this study. The quality of the numerical results is sensitive to the grid resolution. In that sense, an adaptive mesh refinement has been applied. In this work, several flood events which took place in the Ebro river have been simulated to calibrate and compare the numerical model. Nevertheless, here only the 2015 event is shown since it is the most important in recent years. Figure 3 (a) shows the inlet hydrograph used to characterize the discharge evolution during the flood episode. A gauging curve was imposed as outlet boundary condition at the end of the domain (Figure 3 (b)).

Results
With the aim of comparing field and simulated data in the 2015 flood event, some observation points and cross sections have been distributed at strategic locations where gauging stations register water surface elevation (WSE) and discharge (see Figure 1). On the other hand, the convex hull of the observed flooded area during the episode has been provided by CHE and the coincidence between maximum observed and calculated flooded areas is also computed by [6].
where A O and A s stand for the maximum observed and simulated flooded area, reaching a 87.25 % of coincidence.   Figure 5 shows the comparison between measured and calculated data in terms of temporal evolution of WSE and discharge at gauging stations. The continuous line represents the simulated data, while the dotted line illustrates the real measurements. The differences between the measured and computed hydrographs in Figure 5 (d) are justified by the difficulties experienced during that flooding event at the Castejon gauging station (inlet section to our model). The results indicate that the assumed inlet hydrograph does not correspond to the actual inlet flow, hence leading to discrepancies in the discharge at the downstream control sections. This fact was corroborated by the river basin administration (CHE). Table 1 shows the total simulation time and the speed-up reached by the GPU (GeForce GTX TITAN BLACK). The real flood event lasted 21 days. If GPU technology had not been used, the simulation time would have been unaffordable, as seen in Table, where the simulation time using a CPU (Intel CORE I7-4770) parallelized in 8 cores is also written.

Conclusions
Numerical simulation has become an important tool in recent years to help in the prediction and management of flood events. Acceleration technology developments have turned two dimensional models into affordable improving the accuracy of the results. In this work a large flood event has been reproduced by means of a 2D SW model and, after comparing the results, several conclusions can be drawn. Thanks to the model, the results provided by the simulation present a high level of accuracy in terms of WSE and discharge. The variables evolution has been compared with field data and results present a good agreement with observations. The 2D models have demonstrated to provide more information than 1D models, as flooded areas. However, these models have been traditionally rejected due to its high computational cost. In this work, a powerful tool to overcome this problem has been presented. By using GPU technology, important speed-ups have been reached and fast and accurate results have been provided by the model. Not only the mathematical model is important to provide accurate results, but also the available initial data used to represent the terrain. Due to the lack of recent or accurate information the final flooded area reached a 87.25 % of coincidence with the observed measurements, which could be even improved if field data of higher quality were available.
The application of these very useful programming techniques makes possible the simulation of explicit schemes on large spatial resolution cases over long time scales of interest in Hydrology. Moreover, the reduction of the computational time required to perform a simulation makes easier the calibration of parameters that participate in the model.