Are 2D shallow-water solvers fast enough for early flood warning? A comparative assessment on the 2021 Ahr valley flood event

Flash floods pose a distinct challenge compared to traditional fluvial flooding, with infrastructure-based solutions proving less effective. Effective responses hinge on advanced early warning systems providing actionable information, emphasising the necessity for computational flood forecasting models. However, hydrodynamic models, renowned for accuracy and completeness, face limitations due to computational intensity.

This study explores two 2D flood forecasting models, RIM2D and SERGHEI, both with GPU implementations which allow us to maximise the forecast lead time. While RIM2D is less computationally intensive, suitable for operational use, SERGHEI, with higher computational costs, targets large-scale high-performance computing (HPC) systems.

The assessment of applicability and trade-offs is carried out on the 2021 Eifel flood event, particularly in the lower Ahr valley. A set of simulations were performed at various resolutions from 1 to 10 m, which reveal similar accuracy among both models at coarser resolutions, yet discrepancies arise at finer resolutions due to the distinct formulations. Both models exhibit a rapid computational cost escalation, but at resolutions equal to or coarser than 5 m, forecasts are remarkably faster than the real-time ideal for operational use, paving the way for their use in early warning systems. However, higher resolutions necessitate multi-GPU and HPC capabilities, underlining the importance of embracing such technology in addressing broader flood domains.

1 Introduction

The accurate prediction and timely communication of future natural disasters, particularly floods, have become crucial components for disaster management strategies. Early warning systems play a key role in reducing the loss of life and property during such events, allowing appropriate preventive measures to be taken beforehand (Šakić Trogrlić et al., 2022). One of the key tools in these systems is computational hydrodynamic models enabling the simulation and forecasting of flooding in response to varying conditions.

Two-dimensional shallow-water equation (2D SWE) models have been around for quite some time and have been implemented in multiple use cases (e.g. Pasculli et al., 2021). The 2D SWE solvers have a long history in flood modelling (De Almeida and Bates, 2013; Hill et al., 2023) and are a promising approach for enhancing the accuracy and efficiency of early warning systems (Apel et al., 2022; Cea and Costabile, 2022; Costabile et al., 2023). However, until recently, the practical application of 2D SWE models in early warning systems has been very limited due to various challenges related to computational capabilities, data assimilation, and real-time decision-making.

The 2021 flooding event in the Ahr valley (Germany) stands as a stark reminder of the destructive power that extreme weather events can unleash. In July 2021, the region experienced a catastrophic flood event, resulting in loss of life; displacement of residents; and extensive damage to infrastructure, homes, and landscapes (Mohr et al., 2022). Out of the 184 fatalities in Germany, 133 occurred along the river Ahr – a Rhine tributary. The relatively small size of the Ahr river basin (∼900 km²) and its morphological features including narrow streams in gorges result in a stream network with limited capacity for handling sudden influxes of water which consequently makes many areas in the Ahr prone to flash floods. Flash floods are characterised by their sudden onset and fast escalation (Kelsch, 2001).

Catastrophic events such as the Ahr floods are rare and have mostly a local effect, which partially explains why they have received historically less attention than large river floods and likely remain under-represented (Paprotny et al., 2018). However, climate change is likely to make such events more frequent and more intense (Donat et al., 2016; Myhre et al., 2019), thus arguably making them more prominent even in regions in which they have been atypical. From a prevention point of view, regions potentially strongly affected by flash flood events can have very little room for structural improvement. This is the case of the Ahr valley, with urbanised areas occupying the very narrow floodplain and surrounded by steep valleys. With limited potential for structural defences, early warning systems are the key tool to allow the continued safe inhabitation of these areas so that both loss of life and economic damage may be minimised.

Early warning systems pose many challenges. The spatial and temporal scales of flash floods and the consequent short lead times make it challenging to run timely and accurate flash flood simulations producing actionable information (Merz et al., 2020). The Eifel flash floods were a severe stress test for the existing early warning system, which resulted in short lead times, untimely warnings, incomplete/outdated/inaccurate information, and inconsistent recommendations (Thieken et al., 2023b). The nature and timing of the issued flood warnings played a role in the scale of the casualties (Thieken et al., 2023a). Thieken et al. (2023a) argue that warnings communicating rainfall amounts are far less interpretable (by the general population but possibly also by managers and emergency responders) than water levels and inundated areas. However, forecasting water levels, inundated areas, flow velocities, and time of arrival of a flash flood requires, firstly, a hydrodynamic extension of the existing flood forecasts, which are based on hydrological model output at selected river gauge locations, and, secondly, a high level of sophistication in the hydrodynamic flood model employed.

This means that an appropriately high-resolution model is mandatory to capture the complex geometries of valleys, streams, and urban areas in order to reliably predict inundation areas and water levels. Second, the nature and complexity of the physical phenomena do not allow for 1D simplifications, which are far more commonly implemented (Hill et al., 2023) than 2D models. Finally, the simulation needs to be computed fast enough to allow for sufficient lead time. Until recently, this was not achievable and remains the main impediment to the wide-spread adoption of 2D models in flood modelling practice (Hill et al., 2023). However, as 2D SWE solvers are enhanced to more effectively leverage high-performance computing (HPC), new possibilities for early warning with 2D SWE models arise. In general terms HPC has enabled physics-based geoscientific modelling to achieve unprecedented detail (Alexander et al., 2020), and in particular, shallow-water solvers are now fully exploiting this with the use of GPU computing (Morales-Hernández et al., 2020), as well as leveraging massively parallel supercomputing (Caviedes-Voullième et al., 2023a; Morales-Hernández et al., 2021).

The questions that naturally follows are as follows: can HPC-enabled shallow-water solvers achieve sufficient accuracy and lead time to improve early flood warning systems in order to better manage events such as the Ahr valley floods? Does this technology translate into better and more actionable information? We explore these questions using two surface flow solvers, namely the RIM2D and SERGHEI solvers, using different mathematical models and HPC implementations to assess not only the feasibility but the trade-offs.

2 Methods

2.1 Numerical models

We use two 2D surface flow solvers in this work, namely SERGHEI (Caviedes-Voullième et al., 2023a), which solves the fully dynamic shallow-water equations, and RIM2D (Apel et al., 2022), which solves a local inertia approximation. The key advantage of SERGHEI is that it can be deployed on very large-scale HPC systems, leveraging massively parallel scientific hardware. This allows us to offset the comparatively larger computational cost of solving the full shallow-water equations. In contrast, RIM2D allows us to solve the comparatively cheaper local inertia equations, arguably requiring fewer computational resources, albeit in the current version 0.2 limited to a single GPU.

2.1.1 Full shallow-water solver: SERGHEI

SERGHEI (Caviedes-Voullième et al., 2023a) solves the fully dynamic shallow-water equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+{\displaystyle \frac{\partial \mathbf{F}}{\partial x}}+{\displaystyle \frac{\partial \mathbf{G}}{\partial y}}={\mathbf{S}}_{\mathbf{b}}+{\mathbf{S}}_{\mathbf{f}},\\ {\displaystyle }& {\displaystyle }\mathbf{U}=\left[\begin{array}{c}h\\ q_x\\ q_y\end{array}\right],\mathbf{F}=\left[\begin{array}{c}{q}_x\\ \frac{q_x^{\mathrm{2}}}{h}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\\ \frac{q_x{q}_y}{h}\end{array}\right],\mathbf{G}=\left[\begin{array}{c}{q}_y\\ \frac{q_x{q}_y}{h}\\ \frac{q_y^{\mathrm{2}}}{h}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\end{array}\right],\\ {\displaystyle }& {\displaystyle }{\mathbf{S}}_{\mathbf{b}}=\left[\begin{array}{c}\mathrm{0}\\ -gh\frac{\partial z}{\partial x}\\ -gh\frac{\partial z}{\partial y}\end{array}\right],{\mathbf{S}}_{\mathbf{f}}=\left[\begin{array}{c}\mathrm{0}\\ -{\mathit{\sigma }}_x\\ -{\mathit{\sigma }}_y\end{array}\right],\end{array}\end{array}$$

where the conserved variables U are water depth h [L] and momentum components $[L^{\mathrm{2}}/T]$ in the Cartesian directions q_x and q_y. F and G represent the fluxes. S_b is the bed source term, where z is bed elevation [L]. S_f is the friction source term, where σ_x and σ_y are the friction slopes, here computed using Manning's equation. Finally, g is gravitational acceleration $[L/T^{\mathrm{2}}]$.

SERGHEI is written in C++ with hybrid parallelisation, i.e. Message Passing Interface (MPI) for distributed computations and Kokkos for shared memory computations. Kokkos (Trott et al., 2021) is a performance portability layer enabling it to reach both CPU and GPU back ends. Consequently, SERGHEI can run on multiple GPUs and is enabled for large-scale use in large HPC systems.

2.1.2 Local inertia solver: RIM2D

RIM2D is a 2D raster-based hydrodynamic model developed by the “Hydrology” section of the German Research Centre for Geosciences (GFZ) in Potsdam, Germany. RIM2D solves the local inertia approximation to the shallow-water equations (Bates et al., 2010), which has been widely shown to perform well for fluvial floodplain inundation applications (e.g. Falter et al., 2014; Neal et al., 2011; Apel et al., 2022). The local inertia approximation neglects the convective acceleration terms and as a consequence decouples the fluxes in x and y directions. Thus, the fluxes F and G in Eq. (1) reduce to

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }\mathbf{F}=\left[\begin{array}{c}{q}_x\\ \frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\\ \mathrm{0}\end{array}\right],\mathbf{G}=\left[\begin{array}{c}{q}_y\\ \mathrm{0}\\ \frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\end{array}\right].\end{array}$$

Conceptually, the local inertia formulation offers a more precise portrayal of the issue compared to the other simplified version of the SWE equations such as the zero-inertia (diffusive wave) model (De Almeida and Bates, 2013; Caviedes-Voullième et al., 2020). This is because, in contrast to the zero-inertia form, it keeps the local acceleration terms. In the discrete context, this implies that the fluid's momentum in a specific time step informs the subsequent step, thus imprinting a local acceleration in time. Thus, in describing shallow-water flows physically, the local inertia formulation stands as the intermediary between the diffusion wave approximation and the comprehensive full dynamic equations. While the original numerical solution offered by Bates et al. (2010) is susceptible to instabilities under near-critical to super-critical flow conditions and for small grid cell sizes (De Almeida and Bates, 2013), the numerical diffusion proposed by de Almeida et al. (2012) has been additionally implemented in RIM2D.

RIM2D is written in Fortran and ported to GPUs via CUDA Fortran libraries. It is worth noting that presently, RIM2D solely supports computations on a single GPU. However, efforts are underway to incorporate multi-GPU computing capabilities into RIM2D in the near future.

2.2 Study case

The Ahr river is an 86 km long tributary of the Rhine river, located in the states of Rhineland-Palatinate and North Rhine-Westphalia (Germany) in the Eifel region. Our study domain focuses on the downstream reach of the Ahr river, spanning approximately 30 km between the towns of Altenahr and Sinzig. In the first third of the reach the river valley is still very enclosed but opens upstream to the town of Bad Neuenahr-Ahrweiler into a wider valley floor. The area consists of mostly rural areas, with a handful of small settlements and the comparatively larger urban area of Bad Neuenahr-Ahrweiler (population of approximately 26 500) (Truedinger et al., 2023). The average annual precipitation level of the region is below the German mean at around 675 mm (Truedinger et al., 2023).

The nearly stationary low-pressure system “Bernd” resulted in heavy rainfall events in western and central Europe in mid-July 2021 which triggered severe and sudden flooding especially in Belgium, the Netherlands, and Germany (Schäfer et al., 2021). The Ahr valley was one of the locations in Germany which was severely affected, accounting for overall 70 % of all fatalities in Germany (Truedinger et al., 2023), 189 in the area around the Eifel, making it the second largest water-related disaster in recent history in Germany (Thieken et al., 2023b) in terms of casualties. Numerous factors contributed to this extreme impact. Firstly, the Eifel embodies a low-mountain terrain characterised by steep slopes and narrow valleys, extensively settled and cultivated by communities over an extended period. Consequently, the limited space results in a concentration of both population and structures in vulnerable zones. Furthermore, such areas are inherently susceptible to significant issues like mass movement, rapid erosive discharge, and substantial debris accumulation. These conditions notably caused extensive blockages, resulting in the destruction of numerous bridges along the Ahr river in July 2021, exacerbating the flood surge (Truedinger et al., 2023). During the 14 July 2021 event, water levels in the Ahr reached their highest values at the available gauging stations since the beginning of their measurements. Although the exact water levels are unknown, as most gauging stations along the Ahr river were damaged or destroyed during the event, there are estimates of water levels of around 9 m at the Altenahr gauge (Mohr et al., 2022), where the normal water depths of the Ahr are less then 1 m.

2.3 Data and model set-up

2.3.1 Spatial data

Three digital elevation model (DEMs) products provided by the German Federal Agency for Cartography and Geodesy (BKG) were used for model set-up. The datasets DGM1, DGM5, and DGM10 with grid resolutions of 1, 5, and 10 m were available. DGM5 and DGM10 are available finished products published by BKG. The DEMs were directly employed as the foundation for the simulations without undergoing any additional alterations or crafting solutions for potential artefacts or lack of features. This is intentional so that the simulations only rely on readily available datasets. Consequently, the simulations fail to realistically depict the riverbed, instead portraying the average water surface in the Ahr river, which usually measures less than 1 m (Apel et al., 2022). This approach is justified because both models used in this study operate based on water levels as boundary conditions rather than water depths and discharge. As a result, even with the presumed bed elevation, the water levels at the model boundary will consistently remain accurate, ensuring overbank flow and floodplain inundation happen in the correct locations and at the appropriate times. The buildings in the simulation domain were cut out from all three DEMs on the basis of building shape files provided by OpenStreetMap. An example of this can be seen in Fig. 1 (white colouring in the lower panel). Consequently, building surfaces acted as closed reflective boundaries in the simulations.

Figure 1 The red line delineates the boundary of the simulation domain, while the lower panel depicts the topography. In the upper figure, purple points indicate the positions of the Altenahr and Bad Bodendorf gauge stations. The blue line represents the maximum observed flood extent during the flooding event in 2021. Satellite imagery: © Google Earth 2024.

Manning roughness values were assigned to the domain based on the 2020 Germany land cover classification derived from Sentinel-2 data (Riembauer et al., 2021). The databases for the classification are atmospherically corrected Sentinel-2 satellite data (with the MAJA algorithm; data provided by EOC Geoservice of the German Aerospace Centre – DLR) and training data from reference data (e.g. OpenStreetMap) and the Sentinel-2 scenes themselves. This land cover was chosen for this study due to its relatively high grid resolution (10 m). In addition to the mapped land use classes, the main Ahr river channel was added as an additional land category. Based on the literature review, an appropriate Manning roughness value was chosen and assigned to each land cover class in the simulation domain. Table 1 shows the assigned Manning roughness values and the percentage coverage of each land cover type in the simulation domain.

Table 1 Land cover categories, their respective area fraction in the domain, and their corresponding Manning roughness values.

The simulation exercise is performed intentionally in a blind fashion without calibrating parameters such as roughness coefficients. The rationale for this choice is that the objective is to evaluate how feasible the use of these solvers is for early warning, and it cannot be assumed that a comprehensive calibration exercise would be available for every valley that the early warning system oversees. Consequently, a blind approach based on available spatial data and standard parameterisations would be the only choice. Of course calibration would be desirable, but with the typical absence of calibration data (flood mapping) and the occasional need for a quick model set-up in an operational case, an uncalibrated model is rather the standard use case in reality. Therefore we present the uncalibrated simulation results and do not dive into an in-depth model calibration in this study.

2.3.2 Flood event data for the inflow boundary

Inflow to the models is provided by the (official) reconstructed water levels (in metres above sea level) at the Altenahr gauge provided by the flood warning centre of Rhineland-Palatinate (Mohr et al., 2022). The reconstruction is needed because the gauge was destroyed during the 2021 event. For model set-up, observed water levels are assigned to the inflow cells in the domain. These cells are chosen on the river channel on the west boundary of the domain. In order to consider overbank flow, cells neighbouring the river channel and with elevations below the maximum water level of the flood hydrograph were additionally selected. Water depths are assigned to the selected cells only when the river water levels exceed the cell elevation.

It is relevant to point out that a stage hydrograph was selected as an upstream boundary because this is what works natively best with RIM2D. Consequently, for comparability, the same boundary was used in SERGHEI, although SERGHEI can handle an inflow hydrograph.

2.3.3 Observation data for validation

For validation purposes in this study we rely on (i) the documented maximum flood extent provided by the State Agency for the Environment (LfU – Landesamt für Umwelt) of Rhineland-Palatinate, against which we evaluate the model skill in terms of flood extent; (ii) the reconstructed stage hydrograph at Bad Bodendorf (Mohr et al., 2022), against which we compare the arrival time of the flood wave; (iii) and water depths derived from 65 high-water marks reported by residents (Apel et al., 2022), from which water depths were derived to compare against simulated maximum water depths.

2.4 Inundation performance metrics

To quantitatively evaluate flood inundation in a domain, a diverse set of metrics are used to identify over- and under-predictions and their proportions. To compute these metrics, the maximum inundation maps of the simulations are evaluated against each other and the observed flood extent. At first, cells are classified with respect to Table 2. This is done by comparing the simulation results of RIM2D to SERGHEI. In addition, the results of each model are also compared to the observed inundation extent. From each comparison a confusion map is generated. From this map, the total counts of the indices shown in Table 2 are computed and used to calculate the domain-wide inundation metrics shown in Table 3. These metrics are adapted from Wing et al. (2017) and Bernini and Franchini (2013). It is important to note that when contrasting RIM2D with SERGHEI, the outcomes generated by RIM2D are considered observed results, as indicated in Table 2

Table 2 Inundation confusion matrix. Each cell in the domain for a given simulation is compared to the corresponding cell in the observed grid and classified according to this table.

Table 3 Flood inundation performance metrics.

3 Results and discussion

3.1 Computational performance and runtime

One of the core questions of this study is whether these solvers are fast enough for their use in early warning systems. Consequently, we first examine the runtime and computational resources required to perform these simulations.

All simulations reported here were computed on NVIDIA A100 GPUs on the JUWELS Booster supercomputer at the Jülich Supercomputing Centre, as well as in the GFZ Linux Cluster.

Figure 2 shows the absolute (Fig. 2a) and relative simulation runtimes for RIM2D and SERGHEI across the four resolutions (relative to each other in Fig. 2b and relative to the event duration in Fig. 2c). Notably, at coarser resolutions (dx=5 and 10 m), both models result in very short runtimes, clocking in at least 99 times faster than the duration of the 2021 flood event. This level of efficiency renders both models highly suitable for enhancing existing operational flood forecast systems while maintaining exceptional forecast lead times. Consequently, this capability facilitates detailed flood impact forecasting and swift responses.

Figure 2 The absolute simulation runtimes (a), the ratio between the simulation runtimes of RIM2D to SERGHEI (b), and the ratio of the 2021 event duration to the simulation runtimes (c) for the dx=10, 5, 2, and 1 m simulations.

As resolutions become finer, the differences in runtime between the two models become more apparent. SERGHEI, employing multiple GPUs, results in runtimes up to 6 times faster than RIM2D, which in the current version 0.2 relies on a single GPU. At finer resolutions, i.e. large number of grid cells to be computed, the computational requirements surpass the parallel computing capabilities of a single scientific-grade GPU, necessitating multi-GPU implementations and some HPC capabilities for operational deployment. It is also notable that at the dx=10 m resolution RIM2D does exhibit a slightly faster runtime compared to SERGHEI. This can be attributed to its less computationally intensive formulation, additionally indicating one GPU to be adequate for simulations at that resolution.

In terms of the usability of these models for flood early warning, Fig. 2c shows that all simulations were faster than the duration of the event. However, this ratio of event duration to runtime varies between 1 and 400 (for RIM2D) and 10 and 300 (for SERGHEI), depending on the resolution. It is also worth noting that the dx=10, 5, 2, and 1 m resolution models each consist of 1.3, 5.5, 34.7, and 139×10⁶ cells, respectively.

It is relevant to highlight that no specific performance optimisation of the models was carried out for this particular case. Such optimisations could include compiler flags and hardware-based optimisations, domain decomposition strategies, and so on. These can potentially reduce runtimes even further, but they are not necessarily generalisable across cases, software stacks, and hardware. Consequently they are not particularly relevant for the objectives of this study. Nevertheless, such optimisation would be required for operational purposes, which would potentially boost performance even further.

Moreover, continued development in the implementation of the solvers will increase computational efficiency (e.g. by implementing a multi-GPU solver for RIM2D or by dynamically balancing the load across GPUs), so this performance is expected to improve.

3.2 Flood model skill

The flood indicators illustrating the accuracy of both RIM2D and SERGHEI in replicating flooded areas across various simulations are depicted in Fig. 3. Overall, both models demonstrate commendable performance, achieving high scores across all indicators. Notably, they exhibit relatively similar performance at coarser resolutions, but differences become more pronounced at finer resolutions. For instance, when considering the critical success index (CSI), both SERGHEI and RIM2D yield comparable results with CSI values above 0.94 at dx=5 and 10 m resolutions, whereas at finer resolutions (dx=1 and 2 m), the CSI values drop into the eighties, highlighting more discernible disparities.

Figure 3 Comparison of flooded areas with the indices' critical success index (CSI), hit rate (HR), false alarm (FA), error bias (EB), and bias percentage indicator (BPI).

These variations at finer resolutions are evident in the error bias (EB) indicator as well. Specifically, in the 1 and 2 m simulations, the scores for the two models diverge significantly, registering low scores of 54.28 and 16.11, respectively. Notably, the hit rate (HR) indicator stands out as an exception, with scores improving with better resolutions. This disparity is primarily attributable to SERGHEI depicting larger flooded areas in the finer-resolution simulations compared to RIM2D, resulting in a lower false negative (FN) value (as indicated in Table 2) and consequently leading to a higher HR score.

3.3 Maximum flood depth

Figure 4 shows the difference in maximum depth between both models for all four resolutions. Areas in which only one of the models predicts wet areas are categorised. It is important to recall that this is not the difference in water depths at any particular time but the difference in the maximum depths reached during the entire event (which may be predicted at a different time by each solver; see Sect. 3.4). The comparisons behave differently along the valley and are strongly affected by resolution. The narrower valley upstream of Mayschoss has reaches with very large differences in water depth, with SERGHEI predicting water depths up to 2.4 m higher than RIM2D at dx=10 m and up to 4 m with dx=1 m. Near Rech there is a trend of SERGHEI predicting much lower water depths than RIM2D, with larger discrepancies at coarser resolutions. Conversely, upstream of Dernau SERGHEI again predicts higher water depths than RIM2D, but the differences are much smaller, on the order of ∼0.4 to ∼1.4 m depending on the resolution. The differences in this narrow river valley with high water depths and flow velocities in the simulated flood events are likely caused by the different mathematical foundation of the models. Under these flow conditions the neglected convective acceleration in RIM2D might play a substantial role in the flow dynamics.

Figure 4 Difference in maximum water depth between the SERGHEI and RIM2D flood envelopes for all four resolutions. Positive values imply SERGHEI predicts higher maximum depths, and negative values imply RIM2D predicts higher maximum depths. The figure only compares true positive cells (flooded in both models). Gray colours show false positives and false negatives. Note the different ranges and colour scales for each spatial resolution.

Consequently, in Bad Neuenahr-Ahrweiler, where the valley widens and the water depths and flow velocities reduce, the differences are significantly smaller, with a mix of positive and negative differences. In the region around and downstream of Bad Bodendorf SERGHEI tends to predict shallower depths than RIM2D. Additionally, at higher resolution there are more areas which are flooded by SERGHEI than RIM2D than at coarser resolutions. Of particular interest is that going from 5 to 2 m generates additional flooded areas by SERGHEI in Ahrweiler.

3.4 Time to maximum depth (lag)

Figure 5 shows the difference in time to maximum water depth (henceforth lag for brevity) between SERGHEI and RIM2D for the different spatial resolutions used, and Fig. 6 shows the probability density functions of the lag. The lag is computed as follows: for both solvers, the time at which a particular cell reaches the maximum depth during the simulation is registered, and afterwards the difference (lag) between the time obtained by SERGHEI and RIM2D is computed.

Figure 5 Difference in time to maximum water depth (lag) between the SERGHEI and RIM2D flood envelopes for all four resolutions. Negative values imply SERGHEI predicts earlier maximum depths, and positive values imply RIM2D predicts earlier maximum depths. The figure only compares true positive cells (flooded in both models).

There are both positive (RIM2D predicts earlier maximum depths) and negative (SERGHEI predicts earlier maximum depths) lags. Overall, negative lags only occur upstream of Mayschoss in the narrowest part of the river valley. Clearly, the lag mostly increases from upstream to downstream (i.e. delays accumulate downstream). There are some local regions in which this does not hold (e.g. with dx=5 m, between Mayschoss and Rech). The second point is that the lag range reduces with increasing resolution. At 10 m resolution the lags are significant, up to ∼4 h, roughly 8 % of the duration of the event. At 1 m resolution the lag drops to maximums of ∼2 h, roughly 2 % of the event duration.

In the reconstructed water level graph derived from the Bad Bodendorf gauge (Mohr et al., 2022), the highest water level occurs at 27.75 h after the start of the simulation period (14 July 2021), which is 2.5 h after the peak in the inflow hydrograph at Altenahr. Herein we refer to the time difference between the peak at these two stations as hydrograph lag, and we use this 2.5 h value as a reference. We computed the same hydrograph lag between both points for the simulations and report it in Table 4. We also compute the difference between the simulated hydrograph lag and the 2.5 h hydrograph lag estimated by the reconstructed hydrographs. Finally, this difference is expressed as an error relative to the reference hydrograph lag.

Table 4 Simulated hydrograph lag between the Altenahr and Bad-Bodendorf gauges and the difference relative to the 2.5 h lag estimated from the reconstructed hydrographs.

Table 4 shows that the hydrograph lag in SERGHEI reduces significantly with increased resolution, whereas the RIM2D hydrograph lag is far less sensitive. For SERGHEI, the lag difference is always positive; i.e. the peak at Bad Bodendorf is simulated later than the reference in SERGHEI. For RIM2D it is the opposite, it is always negative, meaning that RIM2D simulates a faster peak at Bad Bodendorf than the reference. The relative error is rather constant across resolutions for RIM2D, around −40 %, whereas for SERGHEI, as it is very sensitive to resolution, there are very good results at high resolution but rather poor results at 10 m resolution.

These results suggest that the higher-resolution SERGHEI simulations capture better the flood wave advancement, and decreasing resolution increasingly results in underestimates of the flood wave movement. In contrast, RIM2D seems to overestimate the flood propagation speed but is quite insensitive to resolution. It is worth mentioning that optimising each case individually through individual calibration would very likely lead to improved results because simulated flow velocities and arrival times with different resolutions are sensitive to the roughness parameterisation (Bomers et al., 2019; Caviedes-Voullième et al., 2012; Ozdemir et al., 2013). Our results (together with contextual knowledge from the literature) also suggest that RIM2D may be more sensitive to roughness calibration (which is reasonable since the local inertia simplifications give a somewhat higher weight to the friction model) and that it may be calibrated at a given resolution and results across resolution should improve. In contrast, whereas SERGHEI seems less affected by the lack of calibration, roughness parameters may need to be calibrated for each resolution.

To further explore the difference in the predicted lag beyond a single gauge point, Fig. 6 shows the probability density function of the lag between the SERGHEI and RIM2D flood envelopes across all four resolutions. Negative values in the graphs indicate that SERGHEI forecasts an earlier peak in maximum depths, while positive values mean that RIM2D predicts an earlier peak. The comparison in the figure is limited to true positive cells (i.e. areas flooded in both models).

Figure 6 Probability density function of the difference in time to maximum water depth (lag) between the SERGHEI and RIM2D flood envelopes for all four resolutions for positive and negative lag values. Negative lag values imply SERGHEI predicts earlier maximum depths, and positive lag values imply RIM2D predicts earlier maximum depths. The figure only compares true positive cells (flooded in both models). The lag range is limited to $[-\mathrm{5},\mathrm{5}]$ for readability.

Broadly, the trend indicates that RIM2D consistently forecasts earlier maximum depths compared to SERGHEI across all four resolutions (positive lag values), as already hinted by the lags at the Bad Bodendorf gauge point. The lag between RIM2D and SERGHEI is more pronounced at coarser resolutions than at finer ones. As resolutions become finer, these disparities diminish, and the differences tend to converge toward zero.

The comparisons shown in Fig. 5, in Fig. 6, and with the reconstructed flood hydrograph at the Bad Bodendorf gauge imply that RIM2D simulates faster flood wave propagation speed, which is insensitive to the model resolution. The insensitivity to model resolution can be seen positively. The overestimation of the flood propagation speed, however, needs to be considered when interpreting the results particularly in operational flood response if this roughness parameterisation is used. SERGHEI simulates a flood propagation in line with the reconstructed hydrograph at 1 m resolution and tends to underestimate it with increasingly coarser resolutions. This again is also worth considering when using the model at a particular resolution with the presented roughness parameterisation for a particular purpose.

While studies like Martins et al. (2017) and De Almeida and Bates (2013) suggest that the local inertia approximation results in slower flood propagation speeds compared to the full dynamic equations, it is important to note that Fig. 6 solely depicts the variance in time to reach maximum water depth, which integrates additional processes and not only wave propagation phenomena. Therefore, we argue that this lag disparity should not be construed as a metric for wave propagation. In the evaluation of the flood propagation simulation it is also worth noticing that the reconstruction of the flood hydrograph at Bad Bodendorf is also a hydrodynamic modelling result and thus also prone to errors in terms of water depths and timing and is thus not an absolute quantitative reference for the evaluation of the model results.

3.5 Comparison to maximum flood marks

Figure 7 shows a scatter plot comparison of field observations versus simulated maximum water depths at recorded post-flood observations of maximum water marks on buildings. Figure 8 shows these same points explicitly in space and colour-codes the relative difference between observations and simulations.

Figure 7 Scatterplot comparison of observed maximum water depths against RIM2D and SERGHEI simulations for all four resolutions.

Figure 8 Error (as percentage) in the simulated water depth compared to that observed for RIM2D and SERGHEI for all four resolutions. Positive values indicate an overestimation in the water depth, and negative values show an underestimation. The × marks represent points which fall onto the building footprint rasterised onto the Cartesian grid. Satellite imagery: © Google Earth 2024.

Figure 7 shows a slight general trend to under-predict rather than over-predict water depths across the domain. Most importantly, it allows us to see how the comparison shifts with resolution. R² values clearly overall deteriorate for the coarser resolutions. SERGHEI trends towards underestimating more with coarser resolutions, but the opposite is the case for RIM2D.

RIM2D shows closer agreement with observed water marks with coarser-resolution models. As resolution increases, the discrepancy between simulated and observed depths becomes more pronounced. This is particularly clear from Fig. 7, where it is possible to see how for the 1 and 2 m resolutions RIM2D tends to under-predict maximum depths. At 10 m resolution, there are a few RIM2D points which greatly overestimate observations.

This discrepancy primarily stems from RIM2D's tendency to generate smaller flooded areas in higher-resolution set-ups compared to coarser ones, resulting in under-predicted depths or missing inundation in areas further from the main river channel (see Sect. 3.6 for details), which corresponds to the broad behaviour of the points in space (Fig. 8).

In contrast, for SERGHEI, a distinct trend among the four resolutions is not apparent, and all model configurations tend to produce deviations within a similar range, arguably with better estimations at higher resolution (Fig. 7a). Some additional insights can be drawn by also accounting for the spatial distribution of the comparison points, as shown in Fig. 8. There is a somewhat improving trend towards higher resolution, in which the points located farther from the main river channel which predominantly exhibit under-predicted water depths somewhat improve. In certain cases this is because the predicted inundated area falls short of the location of the points.

Close inspection of the location of the recorded water marks shows that many of the predicted points with the lowest scores, especially at coarser resolution, are the result of poor representation of the buildings in the computational grid. This is illustrated in Fig. 9, where it can be seen that for a coarse resolution (e.g. 10 m) many of the observation points fall in cells which are identified as buildings, although the point itself is not in the building. As resolution increases more of these points fall into valid areas of the computational domain. This is likely to happen since these observation points are often water marks on walls or urban furniture close to buildings. It is important to highlight that these building representation challenges are present in both the RIM2D and SERGHEI simulation scenarios. To offset this issue, we also allow a search for valid (non-building) cells adjacent to the cell containing the observed point. This allows some leeway to capture more points in the analysis. Moreover, aside from the issues relating to observed points, Fig. 9 highlights the effect that resolution can have on properly capturing the complex urban environments, even in a fully inundated area as shown in this image. It highlights that although overall metrics may suggest that the 10 m resolution is sufficient to broadly capture the flood, the inundation dynamics in complex urban built-up areas are prone to errors with raster resolutions too coarse to match the urban complexity. In the presented case study in small towns and villages this appears to be the case at 10 m resolution, but this might be different in other urban fabrics.

Figure 9 Detailed view at Ahrweiler of the maximum water depth at four different resolutions (10, 5, 2, 1 m from top left to bottom right) predicted by SERGHEI, together with the location of some observation points (gray dots) and the building footprints (black lines). White areas are grid cells excluded from computations as they are flagged as buildings. The figure shows how the observation points may fall in cells which are flagged as buildings at coarser resolution.

3.6 Flood evolution

Figure 10 shows a comparison of the evolution of flooded areas with increasing water depths for both solvers and all resolutions. Complementarily, Fig. 11 shows the fraction of flooded area larger than a certain depth threshold relative to the full extent of the flood. In physical terms, Fig. 10f corresponds to the very deeply flooded areas and of course includes the main channel. This is of course a rather small fraction of the flooded area (less than 10 % for most of the simulations, as shown in Fig. 11). In contrast, Fig. 10a reflects most of the flooded area (only excluding areas flooded with less than 5 cm of water). Arguably, water depths below 10 cm only reflect an inconvenience in terms of flood impact. However, depths of around 50 cm already include flooded underground and ground floors in buildings, have transport potential to move unsecured objects, and represent a danger to human life. A very large fraction (between ∼80 % and ∼90 % at the peak) of the flooded areas is indeed flooded with more than 50 cm of water, and up to around 40 % to 50 % of the flooded area exceeds 2 m of depth. This strongly underlines the considerable impact of this flood event.

Figure 10 Inundation areas with water depths above 5, 10, 50, 100, 200, and 500 cm during the dx=1, 2, 5, and 10 m simulations. The values have been measured with a 900 s temporal resolution.

Figure 11 Flood area ratio of water depths above 10, 50, 100, 200, and 500 cm during the dx=1, 2, 5, and 10 m simulations compared to the flooded area of the corresponding simulations with water depth above 5 cm. The values have been measured with a 900 s temporal resolution.

In comparative terms, the flooded-area evolution for all water depths shows similar behaviours, especially in terms of interpreting the results for flood impact and warning. Nonetheless, some deeper reading of the differences proves insightful.

The SERGHEI simulations show a clear trend of decreasing peak flooded areas and delayed peaks with coarser resolution. This is expected and consistent with well-known hydrograph attenuation and delay due to numerical viscosity (diffusion) (Caviedes-Voullième et al., 2012).

Interestingly, for RIM2D the effect of resolution is the opposite as in SERGHEI. The local inertia solution results in higher peak areas for coarser resolutions. Additionally, no significant delay of the peaks is observed in the RIM2D flooded-area curves. The insensitivity in the timing to resolution is consistent with the behaviour of diffusive wave (zero-inertia) formulations as discussed in Sect. 3.4, and these results suggest that the local inertia approach keeps this property. It is possible that roughness calibrations could alleviate this issue.

Comparing across solvers for the same resolution shows that (i) for the coarser grids (5 and 10 m) SERGHEI results in smaller flood extents than RIM2D across all depth thresholds, (ii) for the finer grids (1 and 2 m) SERGHEI results in larger flood extents than RIM2D across all depth thresholds, and (iii) the peak of the flooded-area curves is somewhat earlier for RIM2D than for SERGHEI for all resolutions. Observations (i) and (ii) are explained by the previous discussion on the effects of resolution on the different solvers.

Observation (iii) is consistent with the discussion in Sect. 3.4. Although the lack of convective terms in the local inertia equation typically leads to slower wave propagation in comparison to the full shallow-water equations (De Almeida and Bates, 2013; Martins et al., 2017), this is not reflected in Fig. 10. It is likely that the complex dynamics of wave propagation and flood buffering in the channel and floodplains may play a more significant role than the attenuated wave propagation speeds. Moreover, as noted by De Almeida and Bates (2013), the relevance of this wave slowdown is greater for higher Froude numbers. In this event, the simulations show that most of the flow field experiences sub-critical conditions (in fact, mostly with Froude<0.6). This suggests that the wave slowdown in the local inertia solver may not be very significant except for very local areas with higher Froude numbers.

Another important aspect in the evaluation and discussion of the simulation results above is that all simulations used the same set of roughness parameters. Many studies (e.g. Costabile et al., 2023; De Almeida and Bates, 2013; Pappenberger et al., 2005) emphasised that roughness used in surface flow solvers is not absolute but has to be regarded as effective roughness. This means that roughness is the main calibration parameter for hydraulic models, which can compensate for effects of model formulations (solvers) and model set-ups (resolution) on simulation results (Caviedes-Voullième et al., 2012, 2020; Costabile et al., 2017). With dedicated calibrations of both RIM2D and SERGHEI models for different resolutions, it can be expected to reduce the differences in the model results. However, this is out of the scope of this study, which aims at exploring the differences in simulation results caused by solvers and spatial resolution in 2D hydrodynamic models using standard roughness values, just as a modeller might do when exploring potential floods in a new setting and context (precisely what our simulation exercise was intended to represent). A comprehensive study of the sensitivity of roughness coefficients on the RIM2D local inertia solver is expected as future work.

4 Conclusions

In this study we demonstrate that the state of the art in 2D surface flow modelling currently allows for simulations of flash flood events significantly faster than real time, such as the July 2021 Ahr valley flood event. Evidently, runtime remains a function of the model used (in our case the local inertia solver RIM2D and the full shallow-water solver SERGHEI), the target resolution for the forecast (here between 1 and 10 m) and the computational hardware (here we used between one and eight scientific-grade NVIDIA A100 GPUs).

We show that for this particular event, it is currently possible to generate flash flood forecasts 304 times faster than real time at 10 m resolution and 99 times faster than real time at 5 m resolution. Using HPC resources with SERGHEI it is possible to achieve simulations 8.2 times faster than real time even at 1 m resolution. This holds particular significance, especially regarding the Ahr valley floods, where the type and timing of flood warnings were pivotal in determining the extent of the casualties, together with the shortcomings of existing warning systems (Thieken et al., 2023a), including short lead times, untimely alerts, outdated or inaccurate data, and inconsistent guidance (Thieken et al., 2023b). Traditionally, many areas rely on early warning systems that communicate information primarily based on rainfall amounts and water levels or discharges at a limited number of river gauges. However, the Ahr valley floods highlight the potential limitations of such systems, particularly in scenarios where detailed information on the inundation extent, expected water levels in the inundated areas, and water arrival times is essential for effective response and decision-making. The models employed in this study demonstrate the ability to simulate water levels, inundated areas, and flood propagation with a high level of detail and accuracy.

The detailed analysis of the two solvers applied to a range of different spatial model resolutions using the same set of hydraulic roughness showed large similarities in simulation results in terms of inundation extent and depths. Some differences were also observed in terms of timing of flood peaks and wave propagation. These differences can be explained by the different mathematical foundations of the models (i.e. local inertia formulation versus full shallow-water equations) and the resulting differences in simulated wave propagation and dependencies of simulation results from spatial resolution. Knowing about these differences as laid out in the “Results and discussion” section helps in selecting the appropriate model and spatial resolution for the problem to be studied, as well as for interpreting the results. This mainly accounts for flood propagation speed and flow velocities and less for simulated water depths and flood extent, which are traditionally the main concern in flood forecasts. From a practical point of view of deciding on model complexity for these types of events, our results place the comparative behaviour of the local inertia approximation relative to the full shallow-water equations in the expected ranges, with flood extension being not very sensitive to the selected model, whereas hydrodynamic fields are more sensitive (Caviedes-Voullième et al., 2020; Costabile et al., 2019).

Another relevant practical insight of this study is the value of resolution in flood simulation for early warning. Although the skill metrics are in general terms acceptable for the different resolutions, higher resolution still generally improves results. Broadly, there seems to be some significant change in the behaviours below 5 m resolution, which may be attributable to better-resolved topography or buildings. Arguably, 10 m resolution (although providing good model skill) may be too coarse to provide accurate details in urban areas simply because relevant features are not resolved, and with the computational efficiency shown here, it is absolutely feasible to move to higher resolutions to avoid this risk.

Considering that the presented models are not calibrated and that the inherent uncertainties in flood forecast chains originate from uncertainties in rainfall forecasts and hydrological modelling, the uncertainties introduced by the choice of the hydraulic model and spatial resolution are comparatively low (Apel et al., 2008; Sampson et al., 2014). Thus, the choice of the hydraulic model can be rather based on the required simulation runtimes, spatial resolution, and available computational resources rather than on the specific hydraulic properties of a particular solver. For the presented test case in the Ahr valley, spatial resolutions of 5 m and even 10 m would yield forecasts sufficient for actionable flood response, with both solvers providing valid simulation results with simulation runtimes short enough for use in operational flood forecasts. Calibrated models are of course expected to perform even better.

In summary, the key outcome of this work is a proof of concept that this technology is mature enough to be adopted in early warning systems, ensuring sufficient lead time and providing far more informative and actionable results than traditional flood early warning systems. High-resolution and time-resolved depth and velocity fields provide a far better picture of flood severity and allow for additional analytics to derive impact metrics that are much more easily interpretable by the general public, managers, and emergency responders compared to warnings based on communicating, for example, rainfall amounts (Thieken et al., 2023b). The key next steps involve implementing these solvers into workflows which more broadly cover the data and modelling chains for early warning systems.