Exploring heavy-tailed behavior of flood frequency curves: an analysis of dike overtopping effects

1 Introduction

“Heavy tail” is a statistical term used to describe the property of a distribution with a tail (the part of the distribution most distant from the mode) heavier than the tail of an exponential distribution (e.g., Foss et al., 2009; Nair et al., 2020; Papalexiou and Koutsoyiannis, 2013; Wietzke et al., 2020) (Figure 1a). A heavy-tailed distribution implies that extreme events are more likely than predicted under an exponential distribution. This has significant implications for the design of structural flood defenses and risk management. If the underlying (unknown) population distribution is heavy-tailed, but sample distributions are estimated to be light-tailed, this can lead to a severe underestimation of the likelihood of extreme flood events and insufficiently robust infrastructure. Hence, understanding the mechanisms controlling the heavy-tail behavior of flood distributions is crucial for the proper estimation of design values.

Merz et al. (2022) reviewed the factors that cause or contribute to heavy tails of flood peak distributions. Heavy-tail behavior can be related to various atmospheric, catchment, and river system processes, such as rainfall mechanisms, runoff generation, and reservoir construction (Tarasova et al., 2019). Recent studies have further examined these controls: Macdonald et al. (2022) identified event and catchment characteristics governing tail behavior, Wang et al. (2023) showed that heavy tails can be inferred from the nonlinearity of hydrograph recession dynamics, and Macdonald et al. (2025) demonstrated that the spatial variability of rainfall and runoff generation can strongly influence flood peak tail heaviness. Building on these analyses, Macdonald et al. (2024) showed that runoff generation, particularly threshold processes, exerts a strong control on tail heaviness (Viglione et al., 2009; Macdonald et al., 2024), while Wang et al. (2025) mapped the spatial distribution of heavy-tailed flood behavior across Europe and the United States from common streamflow dynamics. While the influence of atmospheric and catchment processes on the upper tail of flood peak distributions has been the focus of many studies, the relation between river system processes and heavy-tail behavior of flood peak distributions is much less well investigated (Merz et al., 2022).

Merz et al. (2022) formulated three hypotheses how river system processes affect the upper tail of flood peak distributions. First, construction of reservoirs tends to increase the tail heaviness. Reservoirs can substantially reduce the peaks of small and medium floods, but have little effect on extreme floods, thus increasing the upper tail of the flood frequency curve (Assani et al., 2006; Mei et al., 2017; Volpi et al., 2018). The difference between the flood frequency curves of natural and regulated rivers depends on factors such as the reservoir location on the river, its storage capacity, and release capacity (Ayalew et al., 2013). Second, confluences in the river network are locations where a significant change in the flow occurs, potentially affecting the upper tail. For example, heavy-tail behavior would occur when the flood waves from the upstream rivers that meet at a confluence would coincide for large floods, but not for small floods. However, the hypothesis that flood wave superposition at confluences leads to the downstream heavy-tail behavior has limited evidence and low plausibility (Geertsema et al., 2018; Guse et al., 2020; Merz et al., 2022). Third, when rivers are diked to prevent flooding, this can notably impact the flood hydrographs. The reduction in the river cross-sectional area can cause water levels to rise, resulting in a faster flood wave propagation and less attenuation. Consequently, the downstream flood peaks may rise resulting in a heavier tail compared to the pre-dike situation (Merz et al., 2022). During extreme floods when dikes fail and a considerable volume of water is retained in the floodplain, the pre-existing natural conditions may be reinstated, potentially leading to a reduction of flood peak (Farrag et al., 2022) and consequently transforming the upper tail of the flood frequency curve back to the pre-dike situation with a lighter tail (Figure 1b).

Only a few studies have investigated the processes of dike overtopping when river water levels exceed the crest height. Depending on various factors such as the design, condition, and maintenance of the dike, the hydraulic conditions of the flood, and the duration of the overtopping, the dike may either withstand the flood or breach, leading to the uncontrolled flow of water (Vorogushyn et al., 2009). Overtopping and dike breach both lead to a displacement of water from the river to the floodplain that can significantly influence the flood frequency distribution, particularly at the upper tail. However, the heavier tail behavior occurs only till a certain point, the dike breach point or the moment when overtopping occurs (Merz et al., 2022), after which the hydraulic load will be relieved, and the heavy-tail effect is not significant anymore. Only a few studies have analysed the effects of dike breaching (Apel et al., 2009; Vorogushyn et al., 2010) and overtopping (Farrag et al., 2022) at upstream locations on downstream flood peaks. Apel et al. (2009) compared the flood frequency curves at different locations along the Lower Rhine based on observed streamflow and simulated peaks. They found that the retention effect caused by dike breaching only appeared in the flood frequency curves derived from the simulation model, as this model explicitly considered the retention effect of dike breaching. Modeled data resulted in significantly lower peaks downstream of the breach location.

The inclusion of dike overtopping and breaching, and their impact on downstream areas, is usually not considered when extrapolating flood frequency curves to large return periods, whether by statistical or model-based methods. This can lead to an overestimation of tail heaviness (Figure 1b) resulting in unnecessarily high costs of flood protection infrastructure. In this study, we use a simulation-based approach to investigate the changes in the upper tail of flood peak distributions as a result of dike overtopping. We hypothesize that the tail heaviness is reduced by dike overtopping and floodplain inundation, and that the strength of this effect depends on the overtopping volume upstream of the point of interest. We go beyond the investigation of Apel et al. (2009) by considering a large number of gauges in a large river basin, thus accounting for the effects of hydraulic interactions and interplay of tributaries in a complex river network. We use the coupled hydrologic-hydraulic modeling system RFM (Regional Flood Model), which is capable of simulating dike overtopping and inundation. To be able to robustly estimate extreme events and tails of flood peak distributions, we use a multi-site weather generator to produce very long time series of precipitation and temperature for our study area. These synthetic meteorological fields are used to drive the hydrologic-hydraulic modeling system. At many locations in the river system, we compare the flood frequency curves, and specifically their tail behavior, resulting from two scenarios: the 1D scenario, where only the 1D hydrodynamic model is considered, and the 1D-2D scenario, where the 2D hinterland model is activated once the river water level exceeds the dike level. The first scenario thus represents the situation where a model is used to estimate extreme floods without considering the retention effects of dike overtopping. This situation also mimics the statistical extrapolation of a frequency curve to high return periods. The proposed model setup allows to investigate how dike overtopping affects the upper tail of flood peak distributions and how this effect depends on the overtopping volume. Unlike earlier studies that examined individual reaches or relied on statistical extrapolation alone, the novelty of this work lies in quantifying the effect of dike overtopping on tail heaviness across a large, hydraulically complex river network, explicitly resolving floodplain retention and the interplay of tributaries, and using synthetic series long enough to constrain the far tail.

2 Methods

To investigate the effect of dike overtopping on tail heaviness at a regional scale, we use the derived flood frequency analysis approach (Blazkova and Beven, 2004; Falter et al., 2015; Lamb et al., 2010). We use the Regional Flood Model (RFM) (Falter et al., 2016; Farrag et al., 2021, 2022) to simulate long-term synthetic flood time series with and without dike overtopping. RFM is a comprehensive model chain (Figure 2) for simulating flood processes and includes the Regional Weather Generator (RWG) (Hundecha et al., 2009; Nguyen et al., 2021), the mHM rainfall-runoff model (Kumar et al., 2013; Samaniego et al., 2010; Thober et al., 2019) and the coupled 1D-2D hydrodynamic model RIM (Farrag et al., 2022), which also considers dike overtopping. Finally, RFM contains the multivariate flood damage estimation model FLEMO (Kreibich et al., 2010; Thieken et al., 2008), which is, however, not used in this study.

2.1 Regional weather generator (RWG)

The RWG is a multisite, multivariate stochastic weather model that generates daily precipitation time series at multiple locations using a multivariate auto-regressive model (MAR-1). It also generates other important variables, including daily minimum, average, and maximum temperature conditioned on precipitation state (wet/dry). Temperature variables are used to derive potential evapotranspiration. Developed by Hundecha et al. (2009) and further advanced and evaluated by Nguyen et al. (2021), RWG uses an extended Generalized Pareto distribution (Naveau et al., 2016) for characterizing at-site precipitation. The model also accounts for spatial and temporal dependence of precipitation. Marginal precipitation distributions are parameterized on a monthly basis to account for seasonal variability. More detailed information on RWG can be found from Nguyen et al. (2021).

2.2 Distributed hydrological model (mHM)

In this study, we use the mesoscale hydrological model mHM (Samaniego et al., 2010), a grid-based hydrologic model that has been extensively tested in various river basins. The model incorporates all major hydrological processes, namely canopy interception, snow accumulation and melting, soil moisture dynamics, evapotranspiration, surface and subsurface runoff generation, deep percolation, and flood routing in river reaches. The runoff generation process is conceptualized as two reservoirs linked in series, with the upper one accounting for fast and slow interflow components, while the lower one, fed by the upper reservoir through a percolation flux, is responsible for base-flow generation. The Muskingum routing algorithm is used to route the total runoff generated at every grid cell to its downstream neighboring cell. In order to capture the spatiotemporal dynamics of hydrological processes, over 20 parameters need to be estimated through calibration. mHM uses the multiscale parameter regionalization (MPR) approach (Samaniego et al., 2010) to estimate model parameters based on input data through non-linear transfer functions. In this way, an optimal, spatially consistent parameterization is achieved by optimizing the parameters of the transfer functions. For more information on the model, the reader may refer to Kumar et al. (2013), Samaniego et al. (2010), and Thober et al. (2019).

2.3 Regional inundation model (RIM)

The Regional Inundation Model (RIM) integrates a 1D hydrodynamic routing model (RIM1D) and a 2D hinterland inundation model (RIM2D). RIM1D simulates the flood wave propagation in river channels and on floodplains between river dikes or elevated banks. It solves the diffusive wave approximation of the shallow water equations (SWEs; Equations 1 and 2) with the first two terms in Equation 1, local and convective acceleration, respectively, being neglected.Where h is the water depth (m), A is the cross-section area (m²), g is the gravitational acceleration (ms⁻²), S o is the bed slope (m⁻¹), S f is the friction slope (m⁻¹), x is the distance (m), and t is time (s), and q is the net lateral flux at the specific cross-section (m³ s⁻¹), defined as difference between lateral inflow and overtopping flow. An explicit finite difference scheme is used to solve the equations with adaptive time stepping conditioned on the Courant–Friedrichs–Lewy stability criterion. The 1D river domain is represented by simplified 8-point cross-sections (Farrag et al., 2022).

∂ Q ∂ t + ∂ ∂ x ( Q 2 A ) + gA ∂ h ∂ x − gA ( S o − S f ) = 0 (1)

∂ Q ∂ x + ∂ A ∂ t = q (2)

We couple the distributed hydrological model mHM with the hydrodynamic routing model RIM1D to route the flow continuously through the entire river network. First, we run mHM separately, where the runoff is calculated for each cell in the catchment and then routed using the Muskingum method to the next neighboring downstream cell following the topography of the DEM (Figure 3a). To avoid routing the water twice using Muskingum in mHM and using RIM1D, the previous step is carried out for all cells in the domain except those overlapping the 1D domain used in RIM1D (grey cells in Figure 3b). The discharge calculated by mHM is then fed into the RIM routing scheme in two ways: (i) as upstream boundary conditions at the headwater nodes of the river network and (ii) as lateral flow added directly to the nearest cross-section (Figure 3c) and represented as q in Equation 2.

In case overtopping flow over dikes or elevated banks occurs in RIM1D, the overtopping volume is estimated by the broad-crested weir formula and used as boundary condition for the RIM2D model to simulate floodplain inundation. RIM2D explicitly solves the inertia formulation of the SWEs on a regular grid with separate treatment of fluxes in x and y direction. The momentum equation of the inertia formulations disregards the convective acceleration term (second term) in Equation 2. RIM2D is efficiently parallelized in Graphical Processor Units (GPUs) and can dynamically simulate large-scale inundation events (Falter et al., 2016). For more information about the 1D computational engine and the 2D hinterland inundation model used in RIM, readers are referred to Farrag et al. (2022) and Apel et al. (2022).

2.4 Analysis of tail heaviness

The RFM model chain produces long streamflow time series at many locations in the river system, from which flood time series are obtained. Several methods exist to quantify the tail behavior of flood time series. (e.g., El Adlouni et al., 2008). These include estimating the shape parameter of the Generalized Extreme Value distribution (GEV) of the annual maximum series (AMS) as a tail indicator (Bernardara et al., 2008; Farquharson et al., 1992; Villarini and Smith, 2010; Macdonald et al., 2022), the shape parameter of the Generalized Pareto Distribution for the Peak-over-Threshold series (Coles, 2001; Fischer, 2018), the L-coefficient of skewness (McCuen and Smith, 2008; Salinas et al., 2014), the upper tail ratio (ratio of the flood-of-record and the magnitude of the 10-year flood) (Smith et al., 2018) and the Gini index (Wietzke et al., 2020). In addition, Embrechts et al. (1997) have used the generalized Hill ratio plot and the mean excess plot as a graphical method to estimate the upper tail behavior from data. Among these options, the GEV shape parameter has several favorable properties (Wietzke et al., 2020) and is probably the estimator most widely used in hydro-meteorological applications.

In this study, we apply the Generalized Extreme Value (GEV) distribution (Equation 3) to analyze the flood frequency of the annual maximum discharge series.Where δ, ζ, and ξ denote the scale, location, and shape parameters, respectively, with ξ controlling the tail behavior. Positive values of ξ correspond to the heavy-tailed Fréchet distribution, which is a special case of GEV, whereas non-positive values indicate light-tailed distributions. The GEV parameters are estimated using the L-moments method in R.

F ( x ) = { e − ( 1 + ξ ( x − ζ δ ) ) − 1 ξ , ξ ≠ 01 + ξ ( x − ζ δ ) ≥ 0 e − e − ( x − ζ δ ) , ξ = 0 (3)

Estimating the upper tail behavior is challenging for various reasons. Particularly in the context of extreme events, the limited sample size makes it difficult to accurately estimate the upper tail, especially if the distribution is heavy-tailed. As a result, the return period of rare events can be overestimated or underestimated (Davidson, 2012; Merz et al., 2022; Papalexiou and Koutsoyiannis, 2013; Wietzke et al., 2020). The non-stationarity of the underlying hydrological or environmental conditions can also affect the upper tail (Sivapalan et al., 2005). For example, climate change can alter the frequency and intensity of extreme rainfall, transferring into changes in flood frequencies (Ansa Thasneem et al., 2021). The fundamental assumption of stationarity may in fact be violated but is often difficult to detect and account for. Further, the underlying model assumptions (the choice of the probability distribution) used to fit the annual maximum discharge data can affect the estimation of the upper tail. Rogger et al. (2012) suggested that flood peak distributions may have complex shapes; for example, most flood events in a catchment might be generated by subsurface stormflow, while for events with extremely high rainfall volume, the catchment storage capacity might be exceeded leading to saturation excess overland flow. This change in runoff generation processes might lead to a step change in the slope of the flood frequency curve. Such effects can only be observed when the sample of events is sufficiently large to obtain such process changes. Our approach of generating very long time series, including events with and without dike overtopping, thus helps in understanding the upper tail behavior.

3 Setup of the regional flood model for the Rhine basin

3.1 Case study area

The Rhine River basin, with a total area of about 185,000 km², is one of the largest river basins in Europe, extending across Austria, Switzerland, Germany, France, Netherlands, and Belgium. The extensive network of tributaries drains the Rhine basin from its headwaters in the Alps down to the lowlands in north-western Germany and the Netherlands (Figure 4). The Rhine has experienced several catastrophic flood events, for example in 1993, 1995 and 2021, which caused significant socio-economic damage and fatalities and prompted extensive flood management efforts (Bubeck et al., 2012; Birkmann et al., 2023). The countries within the Rhine River basin have implemented various measures to mitigate the impacts of flooding. International cooperation plays a crucial role in managing flood risks in the Rhine River basin. The International Commission for the Protection of the Rhine (ICPR) was established to coordinate efforts among the basin countries. A wide range of flood protection infrastructure, such as dikes, flood walls, and controlled and uncontrolled flood retention basins were implemented to regulate the river flow during high rainfall or snowmelt periods (Vorogushyn and Merz, 2013; Kreibich et al., 2017).

3.2 Setup and evaluation of the regional weather generator RWG

The RWG is set up for a very large area of approximately 600,000 km², that includes all of Germany and parts of neighboring countries, and spans across five major river basins: Ems, Weser, Upper Danube, Elbe, and the Rhine (Nguyen et al., 2021). Here, we use the daily gridded E-OBS dataset v. 25.0e (Cornes et al., 2018) with 0.25° spatial resolution for the period from 19,502,021 to set up and calibrate the stochastic model. The marginal distributions for daily precipitation and temperature (minimum, average, and maximum) are fitted on a monthly basis at 528 grid nodes of the E-OBS dataset. This RWG setup was extensively evaluated with regard to various performance metrics such as wet/dry day probability, simulation of average and extreme precipitation and temperature overall showing very good performance (Nguyen et al., 2021). The spatial correlation was however found to slightly overestimate the observed correlation, which can lead to an overestimation of the areal precipitation with increasing spatial scale (Nguyen et al., 2021). We use RWG to generate 3,000 years of synthetic daily time series of precipitation and temperature to force the hydrological model mHM.

3.3 Setup and calibration of the hydrological model mHM

mHM is set up and calibrated for the same areas as RWG, i.e., entire Germany and riparian countries covering the five major river basins of the Rhine, Elbe, Upper Danube, Weser, and Ems (Guse et al., 2024). It contains three spatial layers. Level-0 represents the physiographic characteristics, such as elevation, soil properties and land cover, with a spatial resolution of 500 m. For this, the European Digital Elevation Model (EU-DEM; www.eea.europa.eu/data-and-maps/data/eu-dem) is aggregated as described in Samaniego et al. (2019). The meteorological data of level-2 are provided at 0.25° (daily E-OBS product, Haylock et al., 2008). Level-1 represents the simulation of the hydrological processes with a resolution of 5 km. This spatial resolution is a trade-off between a high spatial resolution and a tractable simulation time for the entire simulation domain.

One spatially consistent parameter set resulting from the MPR-based calibration is used for the entire simulation domain. 13 gauges distributed across Germany and covering dry and wet catchments are used for model optimization using a modified version of NSE (mNSE) as suggested by Hundecha and Merz (2012). 10,000 parameter sets are sampled based on Latin-Hypercube Sampling to identify the best performing parameter set among all calibration gauges. The model performance is evaluated at 141 gauges, with a focus on medium to large catchments (> 500 km²) due to the model resolution of 5 km. The model performance is evaluated with average values for mNSE of 0.57 and KGE of 0.61 for these 141 gauges.

3.4 Setup and calibration of the river network model RIM1D

RIM1D is implemented for the river network shown in Figure 4, following the setup of Farrag et al. (2022). The overbank cross-section geometry, derived from the 10 m digital elevation model (DEM) using a semi-automatic supervised procedure (Farrag et al., 2022), is considered in this study. Cross-sections cover the river floodplain between dikes or elevated banks and stretch from the dike crest (or elevated bank) to the opposite dike crest (or elevated bank) normal to the flow direction. Hence, dike height is explicitly considered in the cross-section geometry, and cross-section edge points represent potential dike overtopping locations. The bankfull depth (dBF), which is not represented in the DEM, is initially estimated using a decision tree model trained on more than 1,200 surveyed cross-sections provided by the Federal Institute of Hydrology (BfG) for the Rhine, Mosel, and Necker River. Given the uncertainty related to these estimates and the uncertainties of the DEM, we then calibrate the estimated bankfull depth by varying dBF in the range ± 0.5 m to maximize mNSE (Farrag et al., 2022) for the simulated discharge and water level within the calibration period. The performance ranges between mNSE of 0.51 and 0.89 for peak discharge and 0.4 and 0.8 for water level. For more details on the calibration procedure of the bankfull depth and the dike height of the cross-sections, readers are referred to Farrag et al. (2022).

3.5 Setup of the inundation model RIM2D

RIM2D is set up for the entire area behind dikes and elevated banks based on the DEM with 100 m resolution. A uniform floodplain roughness value of 0.03 is used as described in Farrag et al. (2022). For each cross-section of the RIM1D model, interface cells in RIM2D are defined which receive overtopping flow volume as boundary condition. The overtopping flow is computed using the broad-crested weir equation. This volume is subtracted from the 1D modeling domain in case of the coupled 1D-2D simulation. The overtopping flow stops as soon as the water level in the interface cells in the RIM2D domain is equal to the water level in the corresponding RIM1D cross-section or the water level in the RIM1D drops below the dike crest height. In the course of the continuous simulation, the RIM2D model is triggered only in case of overtopping events. It is initialized with the dry state and the water depths are set to zero 10 days after the last dike overtopping (Farrag et al., 2022).

3.6 Design of computational experiments

To investigate the effects of dike overtopping along the river network on the tail heaviness, we compare two simulation scenarios, the 1D and the 1D-2D scenarios. In the 1D scenario, only the RIM1D hydrodynamic model is used to route the flow along the river network without allowing any water to overtop the dikes and spill into the floodplain. In this case, the water depth in the cross-section can build up without a limit. Cross-sections are treated as having vertical walls at their edges. In the 1D-2D scenario, the RIM1D module compares the water level at each cross-section with the dike level and activates the overtopping flow calculation in case of exceedance.

Daily meteorological fields of 3,000 years, generated by RWG, are used to force the mHM model. mHM calculates flow hydrographs at various points along the Rhine River network, which are used, in turn, as boundary conditions and lateral flow for RIM1D. We then use RIM to simulate 3,000 years of daily discharge across the Rhine basin for the 1D and 1D-2D scenarios. The 3,000 years are split into 50 realizations of 60 years to obtain time series that are representative of a typical gauge record length. For each of the 50 timeseries, we derive the annual maxima and fit a GEV distribution to estimate the GEV shape parameter. These 50 values thus provide a measure of sampling uncertainty from the ensemble. We finally compare the statistical properties of the shape parameter for both scenarios along the main river course considering the confluences with tributaries. In this way, we gain insight into the effects of dike overtopping on the tail behavior of downstream flood peaks and the role of the overtopping volume.

4 Results

4.1 Influence of dike overtopping on the shape parameter

To test our hypothesis that dike overtopping and floodplain inundation reduce tail heaviness, we compare the GEV shape parameters for the 1D scenario and the 1D-2D scenario. In the 1D scenario, dike overtopping is ignored, while in the 1D-2D scenario, dikes can be overtopped, thus reducing the downstream flood peaks. Taking this hydrodynamic interaction into account substantially reduces the GEV shape parameter (Figure 5). For the 1D scenario, most of the gauges and realizations show positive shape parameters, i.e., heavy tails. This is particularly true for the gauges along the tributaries. The gauges on the main stem of the Rhine show smaller shape parameters. The two most upstream gauges, Maxau and Speyer, show negative shape parameters for all realizations, indicating a bounded tail. For the other gauges along the main Rhine (from Worms to Rees), the sign of the shape parameter varies with a mixture of positive and negative values. For the 1D-2D scenario, the shape parameters are notably reduced and often drop to negative values for most of the gauges. Some gauges, such as Kemmern, Rockenau or Marburg, are less affected by the overtopping effect and still show a considerable number of realizations with large positive shape parameters.

Overall, every gauge shows a positive mean difference between the 1D and 1D-2D scenarios, with an average reduction of 0.11 in the shape parameter. The largest decreases occur at tributary gauges such as Boos (Nahe), Kemmern (Main), Marburg (Lahn) (mean Δ ≈ 0.15–0.20, see Supplementary Table A-1), whereas far-upstream main-stem gauges (Maxau, Speyer) show only modest changes (Δ ≈ 0.05).

Figure 6 presents the mean GEV shape parameter for all gauges under both scenarios. The 1D scenario is consistently above the 1D-2D scenario, confirming the decrease in tail heaviness due to overtopping seen in Figure 5. The decrease in the mean shape parameter is 0.12, when averaged across all gauges. The decrease is smallest (≈ 0.05–0.08) along the heavily regulated Main River, and largest (≈ 0.12–0.18) at unregulated, flood-prone tributaries such as the upper Moselle, Nahe and Lahn (for details see Supplementary Table A-1). Far-upstream main-stem gauges (Maxau, Speyer) show only modest changes (≈ 0.05).

In the Moselle and Main tributaries, a clear upstream-downstream gradient is visible: the mean shape parameter is highest at the most upstream gauges and steadily declines downstream. This pattern aligns with the findings of Macdonald et al. (2022), who noted that smaller, upstream catchments tend to exhibit heavier tails, although catchment size alone provides only a weak explanatory control. In contrast, the Neckar, Nahe, and Lahn tributaries do not show such a gradient, suggesting more complex or heterogeneous influences on tail behavior. Along the main Rhine channel, the two most upstream gauges (Maxau and Speyer) exhibit low or negative shape parameters. A sharp positive jump in the shape parameter is observed at Worms—just downstream of the Neckar confluence—indicating enhanced flood extremes, likely due to tributary contributions. This increase is seen in both scenarios, but is much smaller in the 1D-2D scenario, which considers overtopping.

Further downstream, the shape parameters vary moderately but remain below those seen in upstream tributaries. The Main River, in particular, shows remarkably stable shape parameters along its course, which can be attributed to the presence of 34 river barrages that regulate flow for navigation and hydropower. Despite these differences, the consistent finding across all gauges is that the 1D-2D scenario yields lower mean shape parameters than the 1D scenario, reaffirming the dampening effect of overtopping and floodplain retention.

To assess the variability and distribution of the GEV shape parameter across realizations, we compare results from the 1D and 1D-2D scenarios for key gauges along the main Rhine River (Figure 7), combining kernel density estimates, boxplots, and jittered individual realizations. This combined visualization provides a comprehensive picture of the sampling uncertainty and tail behavior in both scenarios, which is especially valuable when comparing ensembles of model simulations to point estimates.

Most gauges exhibit approximately symmetrical distributions with substantial spread, reflecting the sampling uncertainty inherent in fitting GEV parameters to 60-year time series (see Figure 4 for gauge locations and river network). While the two distributions overlap, especially at mid-reach gauges like Worms and Mainz, the shape parameters from the 1D-2D scenario are consistently lower compared to the 1D scenario (median shape parameters fall by ≈ 0.12 on average across the main-stem gauges). The drop in the median is largest on the unregulated Ruhr and Nahe tributaries (0.14–0.4), and smallest on the barrage-controlled main river (0.06–0.08) (Supplementary Table A-1). The interquartile range (IQR) changes follow the same spatial logic: Ruhr and Nahe stations exhibit the biggest IQR expansions (−0.10 to −0.21), whereas Main-River gauges cluster near zero (IQR Δ ≤ ±0.03). Figure 7 confirm that median values drop everywhere on the Rhine main stem, from +0.035 (Maxau) to +0.139 (Worms). The IQR widens modestly at every main-stem gauge (IQR Δ ≈ −0.01 to −0.04), except at Frankfurt on the Main (where IQR Δ ≈ +0.03, indicating a slight contraction).

In contrast to the universal reduction in median shape parameters, the 1D-2D scenario shows increased variability across realizations at most gauges (25 out of 30), with the IQR widening particularly strongly at unregulated tributaries such as the Ruhr and Nahe. The exception is the heavily regulated Main River, where barrage control suppresses this additional variability (Supplementary Table A-1).

Notably, historical shape parameter estimates (yellow diamonds) often lie closer to the 1D-2D ensemble, suggesting that simulations including overtopping may better reproduce observed tail behavior at several gauges.

In order to understand the effect of overtopping between consecutive gauges, Figure 8 shows gauge-to-gauge scatter diagrams of the GEV shape parameter, using all 50 realizations for the two scenarios. For the tributary pair Lauffen-Rockenau (Neckar) and the main-stem pair Maxau-Worms in the upstream part of the Rhine basin, points are widely scattered around the 1:1 line. In both panels the red cloud lies systematically below the blue cloud, meaning overtopping reduces the downstream gauge’s shape parameter more strongly than the upstream one. The slope (of the best-fit line) therefore flattens, signaling a downstream dampening of tail heaviness once flood-plain storage is activated. For the Middle Rhine and Lower Rhine pairs (Kleinheubach-Frankfurt, Mainz-Kaub, and Andernach-Cologne); points for both scenarios align tightly along the 1:1 line, showing that neighboring main-stem gauges inherit nearly identical shape parameters. This suggests that overtopping is insignificant along these reaches and has a limited effect on shape parameters of the next downstream gauges.

For the Lower Rhine reach (Cologne-Rees), the relationship remains close to linear but shows greater vertical spread: realizations with highly negative shape at Cologne do not always map to equally negative values at Rees. The red points again sit slightly below the blue cluster, indicating an additional—but modest—tail dampening over this reach.

4.2 Effects of overtopping on the flood frequency curve

Figure 9 shows the flood frequency curves for the 1D and 1D-2D scenarios including the overtopping volume that occurs for each flood event in the area between the specific gauge and the upstream gauge. At gauge Maxau, close to the upstream boundary of the modeling domain, the flood frequency curve of both scenarios is identical, as no overtopping occurs within the 3,000 years of simulation. At gauge Worms, downstream from Maxau, the two curves diverge slightly around the 100-year return period. Overtopping starts with very small volumes and becomes higher with higher return period events, which leads to more divergence between the two curves. The area between the gauges Maxau and Worms is one of the locations where large inundation occurs in the 50 simulated realizations. The lower left panel in Figure 10 shows the inundation resulting from one of the events in the first realization where 400 million m³ of water inundated the floodplain between Maxau and Worms. Further downstream at gauge Kaub, the difference between the two curves is still rather pronounced, although fewer overtopping events with smaller overtopping volumes occur between Kaub and its upstream station. This effect is a consequence of the propagation of the reduction in flood peaks from the overtopping at upstream locations. A similar situation can be seen at the gauges Andernach and Cologne. Between Cologne and Rees, there are many overtopping events with overtopping volume reaching a value of 4,320 million m³ in individual scenarios, increasing the differences in the flood frequency curves for the two scenarios. This reach is not only prone to significant overtopping but also has very high floodplain storage capacity resulting in extensive inundation (Figure 10, upper left and right panels). Figure 9 shows that overtopping reduces extreme discharges only where the associated river reach has appreciable floodplain storage. At Maxau, Kaub, and Cologne the local overtopping volume is negligible and the 1D and 1D-2D curves are very close to each other. By contrast, the reach immediately upstream of Worms stores up to ~10⁹ m³, lowering the 1,000-year discharge by about 34%. A smaller but still clear reduction by about 12 and 15% occurs between Cologne and Rees, respectively, where Lower Rhine floodplains store more than 4 × 10⁹ m³. Intermediate behavior is seen at Andernach, where moderate storage in the Kaub–Andernach reach reduces the 1,000-year discharge by roughly 10%. The higher difference between frequency curves for 1D and 1D-2D models is related to the overtopping volume and the floodplain storage capacity in the immediately upstream reaches of the target gauge.

Figure 10 maps the simulated flood depths for selected events, illustrating how overtopping manifests spatially and results in inundation areas. Upstream of the gauge Rees (top-left) extended floodplains on both banks inundate up to > 5 m, occupying tens of square kilometers. Between the gauges Düsseldorf and Cologne (right) inundation is confined to isolated meander cut-offs and harbor basins. In the Maxau–Worms segment (lower-left) the floodplain near Maxau remains dry, but large storage cells develop immediately upstream of Worms. These spatial patterns reveal strongly variable floodplain storage capacities during inundation events and correlate well with the results in Figures 8 and 9.

4.3 Overtopping volume vs. shape parameter

Figure 11 shows the relationship between the cumulative overtopping volume, i.e., the total overtopping volume from the entire upstream river network, and the reduction in the shape parameters resulting from the overtopping in 1D-2D scenario compared to 1D. Since not all realizations result in overtopping or inundation at all gauges, for each gauge, we consider only those realizations where overtopping occurs and calculate the correlation with the shape parameter reduction. The number of points thus differs from one gauge to another.

Figure 11 shows a generally positive relationship between cumulative overtopping volume and the reduction in tail heaviness which is consistent with the analysis in Figure 9 for the overtopping volume immediately upstream of the target gauge. However, the spread of points in Figure 11 for the cumulative volume is extremely high. This highlights that specific processes such as timing of overtopping during the flood wave propagation and flood peak attenuation downstream may play a significant role. Considerable reductions in the shape parameter are obtained with relatively small overtopping volumes in a number of realizations. Points for the gauges Speyer, Kaub, Andernach and Cologne cluster tightly near the origin, indicating little storage with shape parameter reduction mostly below 0.15. Worms and Rees display a broad spread: typical overtopping events store less than 0.4 × 10⁶ m³ and cut the shape parameter somewhat stronger by up to 0.2, while occasional realizations reach overtopping volumes of 3–4 × 10⁶ m³ and shape reduction up to 0.3. We observe two types of outliers in Figure 11. First, realizations with large overtopping volume, e.g., 1–3 × 10⁶ m³ at Cologne and Rees, lead to relatively small shape parameter reduction. Here floodplain storage fills late, mostly during the falling hydrograph limb. This results in total large storage volumes but small peak reduction. This situation may occur when only few locations in the dike lines are overtopped and the inundation process is comparably slow. Second, realizations with small total overtopping volume and large shape parameter reduction, e.g., 0.28 at gauge Rees. In these cases, upstream retention, not local overtopping, has already trimmed the peak.

5 Discussion

Modeling experiments presented in this work, comparing scenarios with a 1D river routing scheme and a coupled 1D-2D scheme allowing for dike overtopping, suggest general but spatially variable decreases in shape parameters of flood peak distributions due to floodplain retention. The substantial reduction in the magnitude of GEV shape parameters along the main Rhine channel, accompanied by increased variability across realizations, suggests that overtopping has a system-wide dampening effect on flood extremes. This supports the idea that dike overtopping and floodplain storage act as a buffer, attenuating individual flood peaks. However, the increased variability in shape parameter estimates across realizations indicates that overtopping also introduces additional event-to-event heterogeneity, widening the uncertainty associated with estimating rare events. The heterogeneous response among tributaries—where some gauges become bounded-tailed while others remain heavy-tailed—points to the controlling role of local factors such as catchment size and floodplain storage capacity. Moreover, the consistent tail dampening at some downstream locations, even without local overtopping, reveals the significance of upstream-to-downstream hydrodynamic transfers. These findings underline that overtopping effects are not purely local but can propagate in a complex river network further downstream, with implications for the spatial design of flood risk management strategies.

The spatial patterns in GEV shape parameters emphasize that both natural and anthropogenic factors modulate the tail behavior of flood distributions across the Rhine basin. The upstream-downstream gradients observed in some tributaries (e.g., Moselle, Main) align with the notion that smaller headwater catchments tend to produce heavier tails, but this relationship is inconsistent across the network. The lack of clear gradients in other tributaries (e.g., Neckar, Lahn) suggests that local controls such as topography, land use, and hydraulic structures play a role. Notably, the remarkably stable tail behavior along the Main River highlights the significant influence of river regulation: the extensive system of barrages appears to buffer flood variability and suppress tail heaviness. The sharp step change in the shape parameter at gauge Worms (Figure 6) also points to the effect of tributary contributions (Neckar-Rhine) in shaping flood extremes. A previous study based on observational data, however, found that confluences are generally unlikely to cause heavy-tail behavior of flood peaks downstream (Guse et al., 2020). Our study is based on simulations which include a larger sample of extreme events than observations which might explain this discrepancy.

Along the main Rhine reach, overtopping uniformly lowers the median shape parameter, yet increases variability. This is due to the fact that the retained volume and the decrease in flood peak vary from event to event. The central tendency of the shape parameter drops at every gauge (Figures 5–7), yielding lighter-tailed flood distributions. Second, the ensemble spread typically widens: the inter-quartile range increases at roughly 83% of gauges (25 gauge out of 30), indicating that the range of shape parameter estimates—and thus the width of any 95% confidence band around extreme-flow quantiles—can grow even as the median declines. This trade-off is not uniform across the network. At heavily regulated reaches such as Main, barrage control strongly constrains peak variability, so overtopping slightly narrows the IQR. Elsewhere, especially on unregulated tributaries, the additional storage pathway introduces event-to-event heterogeneity that amplifies uncertainty. These results indicate that considering overtopping will generally reduce design discharges but may also widen the uncertainty― a dual effect that risk managers must reconcile.

The paired scatter plots (Figure 8) highlight how local floodplain processes affect tail heaviness at gauges along the network. Where floodplain volume is large and activated between neighboring gauges (upstream of Worms and downstream of Cologne), there is a consistent reduction in shape parameters due to overtopping along with a considerable spread, i.e., increased variability. In reaches with small storage capacity, a small decrease in shape and no changes in spread is found. Typical inundation patterns along the reaches with constricted and wide floodplains (Figure 10) support the statistical results presented in Figures 5–8. Where storage is patchy or heavily constrained by urban levees (Cologne reach) or virtually absent (Maxau), both the tail behavior and the extreme discharge quantiles remain largely unchanged. This behavior is consistent with Tull and Passalacqua (2025), who showed that floodplain attenuation becomes storage-limited as discharge rises, so that retention is most effective below the point at which the available floodplain storage is exhausted.

Figure 11 suggests that cumulative overtopping volume alone cannot explain the reduction in the shape parameter at specific gauges, despite the overall tendency for shape parameter reduction with the presence of overtopping. Timing effects, such as when or how close to the peak flow an overtopping occurs, control how quickly the outflow volume is generated and affect the peak reduction.

Our results extend prior findings that river system processes, particularly dike overtopping, have a significant influence on the upper tail of flood peak distributions. While earlier studies (Apel et al., 2009; Farrag et al., 2022) highlighted the retention effect of dike failures at selected locations, our study systematically quantifies this effect across a large, complex river network using high-resolution modeling. Compared to Apel et al. (2022), we expand the spatial scope to many gauges along the Rhine and its tributaries, capturing interactions between tributaries and main channel dynamics.

The observed reduction in the GEV shape parameter in the 1D-2D scenario corroborates the hypothesis by Merz et al. (2022) that overtopping and floodplain inundation reduce tail heaviness by mimicking more natural, pre-dike flow conditions. Similar results have been shown in reservoir studies, where peak attenuation due to upstream retention leads to lighter tails (Mei et al., 2017; Ayalew et al., 2013). The spatial patterns we find—e.g., higher shape parameters in smaller, upstream catchments—are also consistent with results by Macdonald et al. (2022), who found that tail heaviness tends to decrease with catchment size.

Our findings are subject to several limitations. First, the modeling framework assumes that overtopping begins strictly when water levels exceed predefined dike thresholds. In reality, dike integrity depends on multiple factors including structural condition, maintenance history, and thus dikes can fail due to other breach mechanisms (Vorogushyn et al., 2009), which are not explicitly modeled.

Second, the digital elevation model (DEM) and estimated bankfull depths used in RIM1D and RIM2D simulations introduce geometric uncertainties. Although we calibrate dike heights and bankfull depths within ±0.5 m ranges, residual errors in floodplain topography can impact the volume and extent of simulated overtopping.

Moreover, the sensitivity of tail heaviness to these assumptions has not been explicitly quantified in this study. Future work could employ ensemble modeling approaches with perturbed dike heights or DEMs to assess how such uncertainties propagate into flood frequency estimates. Likewise, the simulation assumes that the overtopping water does not trigger cascading breaches or secondary flood dynamics, such as an overtopping event triggering a breach at the unprotected dike side, re-entering the river downstream.

Lastly, while we use 3,000-year-long synthetic meteorological series to drive the hydrologic-hydrodynamic model chain, potential biases in the weather generator and the hydrological-hydraulic coupled model may affect the realism of the derived hydrologic extremes. A potential bias of the weather generator in representing the spatial correlation of precipitation might overestimate total catchment event precipitation (Nguyen et al., 2021) and consequently flood peak flows. The relative changes and shape parameter reduction, which are the focus of this study, are likely to be not significantly affected by this overestimation.

6 Conclusion

Our results highlight that ignoring dike overtopping in flood modeling can significantly overestimate the upper tail of flood distributions, particularly at downstream locations. This has important implications for flood risk assessments and infrastructure design. Using model setups that allow for overtopping can yield more realistic estimates of design floods, especially for high return periods.

From a policy perspective, incorporating overtopping and floodplain retention into national and international risk mapping initiatives, e.g., carried out by the EU member states following the EU Flood Directive, could improve the accuracy of hazard estimates and ensure that flood defenses are neither under-designed nor overdesigned. Furthermore, identifying locations with high sensitivity of shape parameters to overtopping volumes could inform targeted investments in controlled retention measures.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frwa.2026.1873312/full#supplementary-material

Mean GEV shape parameter and interquartile range per gauge for the 1D and 1D-2D scenarios. Supplementary material at the source

Footnotes

Edited by: Upaka Rathnayake, Atlantic Technological University, Ireland

Reviewed by: Shoukat Ali Shah, Wuhan University, China

Mariacrocetta Sambito, Kore University of Enna, Italy