An Extended Flood Characteristic Simulation Considering Natural Dependency Structures

Abstract

The design of a river-basin-specific flood hydrograph generator based on gauge records enables the generation of synthetic flood hydrographs for the dimensioning of hydraulic structures. Based on selected flow time series, flood waves can be described using four parameters based on flood characteristic simulations, as described by Leichtfuss and Lohr (1999). After successfully adapting suitable distribution functions, dependencies in the load structure are quantified in the next step using copula functions. This newly developed approach builds on the procedure proposed by Bender and Jensen (2012), which assumes hydrological independence. Using copula functions results in increased accuracy in the extended flood characteristic simulation. Moreover, considerable enhancements are achieved through the utilization of genetic algorithms, wherein the descending branch of the flood hydrograph is adjusted by employing an additional variable factor. Subsequently, any number of synthetic flood hydrographs can be generated by combining these parameters. In keeping with the principle of Monte Carlo simulation, a sufficiently high number of synthetic events results in extreme conditions with a low probability of occurrence being reliably represented. Hence, this endeavor has the potential to enhance result reproducibility and prediction quality. As a result, this expanded approach can be employed to provide dependable assessments regarding inflows aimed at optimizing reservoir capacity, for instance.

1. Introduction

Extreme flood hydrographs are frequently used in water management and hydraulic engineering due to their importance. In the design of hydraulic structures, it is often not sufficient to rely only on peak discharges; naturally, additional factors need to be considered. It is therefore crucial to consider the volume associated with the particular design of a hydrograph. In Germany, for example, the relevant regulations require that both flood peaks and volumes are taken into account for the design of dams [1]. This is mandatory because they provide more accurate information about the inflows that maximize the capacity of the reservoir. It should be noted here that these extreme flood hydrographs are usually derived from synthetic design events, which are rarely observed in practice. Still, synthetic hydrographs offer advantages beyond reservoir design, as they can be used for sizing drainage structures and pumping stations or as input to unsteady flow models, particularly for levee systems.

Currently, there are several methods available for describing synthetic hydrograph shapes. Starting with a description using a triangular shape, the mapping of synthetic hydrographs has been continuously developed and refined [2]. However, these methods are subject to certain limitations or require a large database of information. For instance, the approach proposed by Fouroud and Broughton [3] and Bhuyan et al. [4] is often of minor use due to the catchment area size.

Additionally, several other issues such as precipitation data, knowledge about the catchment area, geomorphological aspects, numerical and statistical analyses of peak flows based on observed values, and the utilization of digital terrain models need to be kept in mind and should be considered additional sources of information [3,4,5,6,7,8,9,10,11].

Table 1. Initial data determining the flood hydrograph.

Author	Initial Parameters
Fouroud and Broughton (1981) [3]	Precipitation data
Bender and Jensen (2012) [5]	Catchment area conditions
Wyncoll and Couldby (2013) [6]	Geomorphological characteristics
Candela et al. (2014) [7]	Land cover data, digital elevation model
Chatzichristaki et al. (2015) [8]	Numerical simulations
Bhuyan et al. (2015) [4]	Statistical analysis based on observation data
Aranda and Garcia-Barthual (2018) [9] Shatnawi and Ibrahim (2022) [10] Fischer and Schumann (2023) [11]

provides an exemplary overview characterizing the aforementioned input parameters for determining a flood hydrograph in each respective method. Additionally, there are further statistical methodologies, such as those used to simulate coastal storm surge events, that are not discussed in the present study [12,13].

Hence, building upon the approach established by Leichtfuss and Lohr, which is similarly employed, for instance, in the “Determination of design discharges according to DIN 19700 in North Rhine-Westphalia” guideline issued by the MUNLV (the Ministry of Environment, Agriculture, and Consumer Protection of the state of North-Rhine Westphalia, Germany), an alternative method has been introduced as a substitute for conventionally intricate precipitation-runoff modeling [5,14,15]. This method aims to generate synthetic design hydrographs based only on observed flow time series. With sufficient data availability in the form of flow time series, this approach can be applied to any application area, regardless of the catchment size. Given that the application to a case study revealed shortcomings in the method’s ability to accurately represent prolonged peak discharges in a natural manner, Bender and Jensen formulated an enhanced approach [5]. This revised approach effectively addresses this weakness of the methodology.

The main objective of this research is to highlight the prevailing shortcomings of both approaches and improve simulation results accordingly. In order to achieve this objective, copula functions and genetic algorithms are employed. The fundamental questions are as follows: How can synthetic scenarios be generated effectively considering the hydrological boundary conditions? And what is the impact of applying genetic algorithms to the mathematical representation of the descending branch?

2. Extended Flood Characteristic Simulation according to Bender and Jensen

The primary flood characteristic simulation methodology proposed by Leichtfuss and Lohr is based on a four-parameter approach, which involves describing the temporal evolution of a flood event using parameters, including the rising duration, t_r; the shape parameter of the rising branch, s_r; the peak discharge, Q_P; and the shape parameter of the descending branch, s_d. Synthetic hydrographs are then generated based on these parameters [14].

However, the enhanced approach introduced by Bender and Jensen further improves flood characteristic simulations by incorporating an additional parameter, namely, the peak duration, t_P. The shape parameters s_r and s_d are dimensionless input parameters used in the hydrograph function and describe the characteristics of the rising and falling branches of the flood wave, respectively. The interaction between all these parameters (rising duration, t_r; the shape parameter of the rising branch, s_r; peak discharge, Q_P; peak duration, t_p; and the shape parameter of the descending branch, s_d) and the resulting shape of the hydrograph is illustrated in Figure 1 [5].

2.1. Hydrograph Functions

Within the flood characteristic simulation, the mathematical portrayal of the hydrograph entails a differentiation between the rising and descending branches of the hydrograph. These components are parameterized using distinct and independent functions and are represented in Equations (1)–(4), as outlined below. The rising branch is modeled using the Kozeny function, while the descending branch follows a hyperbolic function [16,17].

In these equations, t represents time, ${\mathrm{t}}_{\mathrm{r}}$ indicates the duration of the rising branch, $\mathrm{Q}\left(\mathrm{t}\right)$ describes the total discharge at a given time, and ${\mathrm{Q}}_{\mathrm{B}}\left(\mathrm{t}\right)$ represents the base discharge at the same time. ${\mathrm{S}}_{\mathrm{r}}$ and ${\mathrm{s}}_{\mathrm{d}}$ are both shape parameters of the rising and descending branches and are valid between 2 and 10, as well as 0.05 and 4. ${\mathrm{T}}_{\mathrm{S}}$ determines the starting point of the descending branch, b is the axis section, and c is the compression factor with a constant of one (1). ${\mathrm{Q}}_0$ represents the maximum discharge, and ${\mathrm{Q}}_{\mathrm{P}}$ represents the peak discharge of the event.

$$\mathrm{Q}\left(\mathrm{t}\right)=\left[{\mathrm{Q}}_{\mathrm{P}}-{\mathrm{Q}}_{\mathrm{B}}\left(\mathrm{t}\right)\right]\times [({\frac{\mathrm{t}}{{\mathrm{t}}_{\mathrm{r}}})}^{{\mathrm{s}}_{\mathrm{r}}}\times {\mathrm{e}}^{{\mathrm{s}}_{\mathrm{r}}\times \left(1-\frac{\mathrm{t}}{{\mathrm{t}}_{\mathrm{r}}}\right)}]+{\mathrm{Q}}_{\mathrm{B}}(\mathrm{t})$$

(1)

$${\mathrm{Q}}_{\mathrm{D}}\left(\mathrm{t}\right)={\mathrm{Q}}_0\times \left\{1-\left[\raisebox{1ex}{$({{\mathrm{e}}^{\mathrm{k}})}^{\mathrm{c}}-{{(\mathrm{e}}^{-\mathrm{k}})}^{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{$({{\mathrm{e}}^{\mathrm{k}})}^{\mathrm{c}}+{{(\mathrm{e}}^{-\mathrm{k}})}^{\mathrm{c}}$}\right.\right]\right\}$$

(2)

$$\mathrm{k}=\frac{\left({\mathrm{t}}_{\mathrm{s}}+\mathrm{t}\right)}{[\mathrm{b}+{\mathrm{s}}_{\mathrm{d}}\times {(\mathrm{t}}_{\mathrm{s}}+\mathrm{t})]}$$

(3)

$${\mathrm{t}}_{\mathrm{s}}=\frac{1}{\frac{2}{\ln \left(2\times \raisebox{1ex}{${\mathrm{Q}}_0$}\!\left/ \!\raisebox{-1ex}{${\mathrm{Q}}_{\mathrm{P}}$}\right.-1\right)}-{\mathrm{s}}_{\mathrm{d}}}$$

(4)

The equation incorporates the peak discharge (Q_P) and shape parameters (s_r and s_d) as the primary parameters. The other input variables remain nearly constant and are fixed values. The maximum discharge, which needs to be determined at the beginning of the parameter estimation process, is purely a calculated quantity. A standard value of Q₀ = $2\times {\mathrm{Q}}_{\mathrm{P}}\left(\mathrm{H}\mathrm{H}\mathrm{Q}\right)$ is often used as an example [15]. The term “HHQ” refers to the highest known value, which signifies the maximum recorded value for the particular flow time series used.

The peak duration (t_P) generally forms a slightly decreasing or increasing plateau. Consequently, the applied method records all values lower than 0.5% below the flood peak (Q_P). In cases where a stagnant peak discharge is not observed, a peak duration of t_P = 0 is assigned [5].

Based on the existing equations for calculating the respective discharge value, the total volume (V_tot) can be calculated using integrals. In this context, it is necessary to distinguish between the rising branch (V_r, as modeled using the Kozeny function) and the descending branch (V_d, as represented using a hyperbolic function). The volume calculation at the peak can be calculated relatively easily with ${\mathrm{V}}_{\mathrm{P}}={\mathrm{Q}}_{\mathrm{P}}\times {\mathrm{t}}_{\mathrm{P}}$.

Equation (5) describes the calculation of the total volume (from time t = 0 to the last point of the hydrograph, t_end), subdivided into its components, for the respective hydrograph.

$${\mathrm{V}}_{\mathrm{t}\mathrm{o}\mathrm{t}}={\mathrm{Q}}_{\mathrm{P}}\times {\int }_0^{{\mathrm{t}}_{\mathrm{r}}}{\left(\frac{\mathrm{t}}{{\mathrm{t}}_{\mathrm{r}}}\right)}^{{\mathrm{s}}_{\mathrm{r}}}\times {\mathrm{e}}^{(1-\frac{\mathrm{t}}{{\mathrm{t}}_{\mathrm{r}}})\times {\mathrm{s}}_{\mathrm{r}}}\mathrm{d}\mathrm{t}+{\mathrm{Q}}_{\mathrm{P}}\times {\mathrm{t}}_{\mathrm{P}}+{\mathrm{Q}}_0\times {\int }_{{\mathrm{t}}_{\mathrm{S}}}^{{\mathrm{t}}_{\mathrm{e}\mathrm{n}\mathrm{d}}}\left\{1-\left[\frac{{\left({\mathrm{e}}^{\mathrm{k}}\right)}^{\mathrm{c}}-{\left({\mathrm{e}}^{-\mathrm{k}}\right)}^{\mathrm{c}}}{{\left({\mathrm{e}}^{\mathrm{k}}\right)}^{\mathrm{c}}+{\left({\mathrm{e}}^{-\mathrm{k}}\right)}^{\mathrm{c}}}\right]\right\}\mathrm{d}\mathrm{t}$$

(5)

2.2. Parameter Identification and Creation of a Synthetic Hydrograph

In parameterization, flow time series are analyzed for flood events with a peak flow (Q_P) exceeding a certain threshold. According to Leichtfuss and Lohr, a reasonable threshold for parameterization is ${\mathrm{Q}}_{\mathrm{P}}\ge 2\times {\mathrm{Q}}_{\mathrm{M}}$ [14]. As the current study aims to generate extremely rare synthetic hydrographs with a low probability of occurrence, the principle of the annual maximum series (AMS or AMAX) is used to parameterize the hydrographs [18]. Additionally, it is possible to use critical water levels for river levee dimensioning as a threshold for parameterization.

It is important for accurate parameter determination to ensure that the analyzed flood events are hydrologically independent. This ensures that any flood event influenced by a previous event will not affect the representativeness of the watershed in terms of its shape and peak characteristics [19]. The process for generating synthetic design hydrographs can be summarized as follows:

Initially, a suitable gauge dataset with sufficient information is selected. Then, all parameters are determined based on the principle of the annual maximum series, following the guidelines outlined in Section 2.1. Following this approach guarantees that the sample size remains consistent for each parameter. Appropriate distribution functions are then assigned to these parameter samples (t_r, s_r, Q_P, t_p, and s_d). Upon successful fitting, a large number of synthetic hydrographs can be generated by using the principle of Monte Carlo simulation [20]. Including a sufficient number of synthetic events ensures a desired design event (e.g., HQ₁₀₀₀₀) among the generated hydrographs. The probability of occurrence of the synthetic hydrographs is directly related to the probability of occurrence of the flood peak (Q_P).

2.3. Constraints and Enhancement Strategies for the Methodology

The presented procedure was applied to the observed gauge data in Bavaria, Southern Germany (e.g., River Middle Aurach, Emskirchen), spanning the years 1967–2022. This region exhibits diverse hydroclimatic characteristics, encompassing the alpine region in the south, uplands in the north and east, and flatlands in the central and western parts of the area. Furthermore, the methodology was employed in the analysis of the subsequent gauge datasets: the Deggendorf (Kollbach) and Rappoldshofen (Aisch) gauges.

During the implementation, certain deficiencies in the extended flood coefficient simulation became evident. These limitations arose from the aspects described below.

2.3.1. Limitations in Modifying the Descending Branch of the Hydrograph

The previous approach utilized two independent functions—functions (1) and (2)—to assess the fit of the rising and descending branches. The shape parameters of both the rising and descending branches were employed to achieve a high degree of conformity with the observed hydrographs. The goodness of fit between the observed and simulated hydrographs could be assessed, for instance, by utilizing the root-mean-square error (RMSE) [21,22]. Nevertheless, it was found that certain hydrographs exhibited deficiencies in modifying the descending branch. This tendency resulted in an underestimation of the observed hydrograph. Therefore, the current approach introduces an additional parameter, the compression factor, c, which influences the shape of the descending branch and was previously held constant. These parameters (s_r, s_d, and c) are now determined using an iterative, generation-based optimization method.

• Genetic algorithm

The genetic algorithm runs through multiple iterations to approach an optimal solution, with each iteration representing a generation. Each generation consists of a population of solution approaches. To enhance the solution routes, the best solutions are combined to form a new generation. The selection of the best routes is based on a fitness score, which, in this case, reflects the deviation between observed and simulated runoff.

The process (Figure 2) begins with an initial generation where solution possibilities are randomly generated for the rising (s_r) and descending branches (s_d and c). Each member is assigned a route with scores in any order. Elite selection identifies members with the highest fitness and forms the elite pool. From this pool, parents are chosen for the new generation, and each child inherits part of its solution from one parent while filling in the missing part from the other. The population is replenished, and the cycle of fitness calculation and reproduction is repeated over multiple generations. With each generation, the deviation between observed and simulated values diminishes.

Progress eventually reaches a stagnation point when the best possible solution is achieved [23,24].

2.3.2. Seasonal Statistics

In widely used flood statistics, annual maximum flows are commonly used to calculate design variables and quantiles. However, these flows comprise winter and summer events, each with their distinct genesis. In the northern hemisphere, the available data are further divided into the hydrological winter (November–April) and summer periods (May–October).

Various options exist to classify observed floods into different event types. These classifications can be based on factors such as snowmelt from snow cover, catchment characteristics, flood-triggering rainfall amounts, or rainfall durations, among others [11,25]. Deriving such indicators requires data beyond just the runoff series, such as snow depths or precipitation data with a sufficient temporal resolution. However, obtaining such additional data is often not feasible, leading to the subdivision of events into summer and winter categories in the present studies.

For more details on different type-based approaches, one can refer to the given studies in [26,27]. These sources provide additional possibilities for considering different event classifications.

2.3.3. Assumption of Statistical Independence for Individual Parameters

So far, the parameters (t_r, s_r, Q_P, t_p, and s_d) in the collected data have been assumed to be statistically independent of each other [14]. Further investigations should determine the extent to which adjusting the parameters is beneficial for the Monte Carlo simulation. Consequently, the correlation among all existing parameters, including the newly adapted compression factor, c, was subsequently assessed.

• Correlation

To determine the presence of dependencies between descriptive parameters, the correlation between random variables can be assessed using two statistical measures: linear correlation (Pearson’s r) and rank correlation (Kendall’s τ). Both measures range from minus one to one (−1 to 1), with a value of one (1) indicating complete positive dependence and a value of minus one (−1) indicating complete negative dependence [28,29].

Consequently, the dependencies for the Emskirchen gauge (Middle Aurach) are described using Kendall’s τ (Figure 3). The left side shows the correlation matrix for the summer period, while the right side shows the winter period correlation matrix. Although the resulting values are comparable, they exhibit variations in some parameters. This discrepancy arises from the distinct genesis of the flood hydrographs, as explained in Section 2.3.2, and it describes the respective dependency measure between the measured values.

Nonetheless, both matrices reveal correlations of varying magnitudes, which were considered in the previously valid procedure. This limitation can be addressed by incorporating copula functions.

• Copula functions

Copula functions are mathematical tools used to model the dependence structure between correlated random variables, independent of their marginal distributions. In the case of n-dimensional variables, copulas are defined as follows:

$$\mathrm{C}:{\left[\mathrm{0,1}\right]}^{\mathrm{n}}\to \left[\mathrm{0,1}\right]$$

(6)

Sklar’s theorem establishes the relationship between a copula function and an n-dimensional multivariate distribution as follows [30]:

$${\mathrm{F}}_{{\mathrm{X}}_1,{\mathrm{X}}_2,{\dots ,\mathrm{X}}_{\mathrm{n}}}\left({\mathrm{x}}_1,{\mathrm{x}}_2,\dots ,{\mathrm{x}}_{\mathrm{n}}\right)=\mathrm{C}\left[{\mathrm{F}}_{{\mathrm{X}}_1}\left({\mathrm{x}}_1\right),{\mathrm{F}}_{{\mathrm{X}}_2}\left({\mathrm{x}}_2\right),\dots ,{\mathrm{F}}_{{\mathrm{X}}_{\mathrm{n}}}\left({\mathrm{x}}_{\mathrm{n}}\right)\right]$$

(7)

To apply copulas, the marginal distributions of the variables (t_r, s_r, Q_S, t_p, s_d, and c) are first fitted independently to each sample using univariate statistical methods [31]. The observed pairs of variables (X and Y) are subjected to a transformation process that maps them into a range between zero and one. This transformation is achieved by directly applying the obtained distribution functions to the observed values. The corresponding copula model is fitted to these transformed values. Parameter estimation in common copula models often employs Kendall’s τ [32].

For hydraulic and hydrological applications, Archimedean copulas are commonly recommended [13,33]. These single-parameter copula functions are characterized by the parameter θ and can represent a wide range of dependency structures. Commonly employed one-parameter copulas include Clayton, Frank, Joe, and Gumbel copulas [34,35,36,37,38]. Further information on the copulas listed and their generator functions, the value ranges of the copula parameter θ, and the functional relationship between θ and Kendall’s τ can be found, for example, in the following studies [13,31,39].

The selection of appropriate marginal distributions can be determined using the methods described by, e.g., Naghettini [31]. Additionally, the goodness of fit of the copula model should be assessed. Various methods similar to those used for assessing univariate distributions can be employed. Simple measures, including RMSE, between the empirical and parametric copulas, as well as information criteria, such as the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) considering the model parameters, can be used. Further possibilities are also provided by bootstrapping procedures comparing the parametric and non-parametric estimations of copulas based on common tests, such as the Kolmogorov–Smirnov and Cramér von Mises tests [40,41]. It is important to note that the model selected, based on formal fitting tests, may not necessarily be the most hydrologically meaningful. Therefore, hydrological considerations should complement fitting tests. A visual assessment of the overall goodness of fit is also recommended [38].

In addition to analyzing bivariate correlations, it is possible to consider multivariate aspects (e.g., Vines) in which an unlimited number of characteristics can be examined theoretically. Vines are graphical structures that generalize dependency trees and represent joint distributions [42,43,44]. The fundamental concept behind vine copulas involves the creation of high-dimensional copulas through a systematic blending of (conditional) bivariate copulas.

There exist diverse types of Vines [32]. In our scenario, we sequentially investigate dependencies, examining the relationship between the first and second elements, the second and third elements, and so forth. Figure 4 provides a visual representation of building a 3D vine copula with parameters X₁, X₂, and X₃, where their pairwise dependencies are captured through the bivariate copulas ${\mathrm{C}}_{{\mathrm{X}}_{\mathrm{1,2}}}$ and ${\mathrm{C}}_{{\mathrm{X}}_{\mathrm{2,3}}}$. Mathematical mapping is illustrated in the figure [45].

2.3.4. Multi-Peak Flood Events and Uncertainties

In keeping with Bender and Jensen, the methodology employed in this study is applicable to the generation of single-peaked flood hydrographs only [5]. It does not account for the modeling of hydrographs with multiple peaks. The challenge of modeling such complex hydrographs has already been addressed by Klein. By using multivariate statistics, Klein developed a method that enables the superposition of multiple Kozeny functions [45]. In addition, Fischer and Schumann also considered multiple-peaked flood events. In this work, the probability density function of a mixture distribution was used to account for multiple-peaked events for the generation of synthetic design floods [11].

Furthermore, the potential impacts of climate change and its uncertainties are not directly integrated into the methodology. This is because the procedure relies on observed flow time series and, consequently, cannot categorically eliminate uncertainties such as measurement errors. For a more in-depth exploration, we refer the reader to additional literature [46,47].

3. Results

The results focus on two main aspects. Initially, the outcomes of the genetic algorithm, specifically highlighting the enhancements in the goodness of fit of the descending branch achieved through variable combinations of s_d and c. Furthermore, the findings derived from the examination of hydrological dependencies within the data and their impact on the results of the extended flood characteristic simulation are discussed in detail. Both adjustments are based on a direct comparison of the simulation accuracy using the detected flood events from the Emskirchen gauge (Middle Aurach).

3.1. Consideration of Combined Genetic Algorithms in Adapting the Descending Branch

Through the utilization of combined (s_d and c) genetic algorithms, a greater level of fitting accuracy can be attained in the descending branch, as compared to the previously suggested method by Bender and Jensen, wherein solely the shape parameter is fitted (Figure 5). This improvement is demonstrated through a correlation analysis using Kendall’s (τ) method and the Kling–Gupta efficiency (KGE) metric [29,48]. Here, the observed and simulated flood volumes are considered. In Figure 5, both the observed and simulated flood hydrographs are displayed, depicting the peak discharge (Q_P) and volume (V).

In the upper section, the genetically adjusted shape parameter for the descending branch remains constant with a compression factor of c_fix = 1. In contrast, the lower section shows a variable adjustment of the shape parameter and compression factor, c_var, for the hyperbolic function of the descending branch. For both winter and summer events, a clear improvement in the simulated flood volumes becomes evident due to the variable compression factor. For example, KGE for summer events improves from 0.666 with compression factor c_fix = 1 to 0.814 with a variable compression factor. The correlation also improves from 0.89 to 0.951.

Figure 6 illustrates the distinctions in the trajectory of the descending branch, depending on the chosen methodology for the compression factor, using an exemplary flood hydrograph. This highlights the observed flood hydrograph (gray) and the simulated hydrographs with s_d and c_fix (green), as well as s_d and c_var (orange). By utilizing a variable compression factor c_var $(0\le \mathrm{c}\le 1)$ and s_d, there is an improvement in the adaptation of the descending branch to the observed flood hydrograph. This can be justified by acknowledging that both variables have a direct impact on the outcome of the hyperbolic function. With c_fix, it is not always possible to achieve a complete correction solely through adjustments of the shape parameter (s_d) for certain hydrographs. Thus, the fixed adjustment (c_fix) in the example below tends to slightly overestimate the hydrograph in the descending branch.

Figure 7 illustrates a simplified depiction of the procedure for refining the coefficient of variation (c_var) within the range of 0.1 to 1 with two constant shape parameters for the descending branch (s_d = 0.05 and 0.40).

It is worth emphasizing that increasing the bandwidth of s_d can also encompass a significantly broader spectrum of hydrographs. In contrast, this improvement can not only be corrected in the same way via the variable adjustment of the compression factor but can also be amplified. Consequently, it is possible to not only ascertain the extent of the descending branch within a specified period but also facilitate significant stretching of the originally compressed shape, exerting a substantial influence on the rate of the descent of the descending branch. The simultaneous application of both variables leads to an improvement, as exemplified in Figure 5.

3.2. Consideration of the Dependency Structure between Correlated Random Variables in Terms of Copulas

3.2.1. Utilizing D Vines for the Present Application

In the current application, a six-dimensional copula can be constructed as an example according to the concept of Vines based on the correlational parameters (t_r, s_r, Q_P, t_p, s_d, and c) for the winter period [49]. The detailed structure, including the appropriate copulas, is further described by using the subsequent visualization (Figure 8).

To this end, the parameters were tailored using univariate distribution functions. In the current scenario, the GEV, normal, and Weibull distributions were chosen due to their goodness of fit based on the Kolmogorov–Smirnov test [40].

Considering the dependency structure between the correlated random variables in terms of copulas results in a new sample. Adhering to the principles of Monte Carlo simulation, a sufficiently large number of synthetic events can be generated as per Leichtfuss and Lohr [14,20].

3.2.2. Plausibility Check

The enhancement in the proposed approach can be assessed by comparing the replication accuracy based on detected flood events recorded via tide gauges. This involves a comparison between the observed volume and the synthetically simulated values (n = 1000) obtained from both the existing flood characteristic simulation (EFCS) approach without considering natural dependencies and the newly developed enhanced flood characteristic simulation with the copula (EFCSC) approach. A preliminary assessment of the improved results can be gleaned from Figure 9.

The figure displays the volume and peak discharges derived from the observed data (gray), as well as the volume and peak discharges generated synthetically using the EFCS approach (blue) and the newly developed EFCSC approach (orange). Both simulations are based on the adjustment of six variables using genetic algorithms (Section 3.1). The data are further divided into the summer and winter half-year periods.

It is evident that there is a further improvement in the simulation results for both the summer and winter periods. From a computational standpoint, the improvement in the results in the summer period is reflected in an increase in the correlation coefficient from 0.342 (EFCS) to 0.403 (EFCSC). This represents an improvement of 17.8% in the summer and 6.3% in the winter period.

4. Discussion

The improvement in the approach presented is firstly determined through a direct comparison of the accuracy of reproducing detected flood events from the Emskirchen gauge on the Middle Aurach River. In this process, all observed flood hydrographs are parameterized and synthetically reproduced using the extended flood characteristic simulation (EFCS) by Bender and Jensen, as well as the newly developed extended flood characteristic simulation considering natural dependency structures (EFCSC). Notably, to enhance precision in reproducing dependency structures, copula functions are employed as part of the improvement, along with the integration of genetic algorithms.

An initial impression of the enhanced results can be gained by directly overlaying the relationship between the observed and simulated volume and peak discharges, as depicted in Figure 9. This is expected due to the existence of interdependencies between individual parameters that were previously disregarded. Through the application of copula functions, weak (<0.35), moderate (0.36 to 0.67), and strong correlations (>0.68) can be accounted for [50]. This expanded information promptly leads to improved outcomes. Naturally, the degree of improvement strongly relies on the available data.

An additional improvement by incorporating another parameter (compression factor, c) through the adaption of genetic algorithms is shown in Figure 4. This leads to an enhancement in the goodness of fit for the descending branch (Figure 5). Thus, the volume of the flood hydrograph, for instance, in the dimension of hydraulic structures (e.g., weir or river levees) can be associated with higher reliability by integrating this additional parameter.

In conclusion, the presented modification of the extended flood characteristic simulation considering natural dependency structures represents a valid upgrade of the existing method. Although the approach cannot capture multi-peaked hydrographs, the incorporation of an additional parameter (compression factor, c) allows for a significantly more accurate simulation of single-peaked flood events. Moreover, by employing copula functions, the goodness of fit can be further enhanced based on the data foundation for simulated flood hydrographs.

Thus, with the developed methodology, the mapping accuracy can be further improved through the incorporation of copula functions and an additional parameter. Further, it should be noted that, given sufficient data availability in the form of low time series, this method can be applied to any catchment area independently of additional characteristics, e.g., precipitation data, knowledge about the catchment area, geomorphological aspects, and digital terrain models (refer to Table 1). This constitutes a noteworthy advantage; furthermore, a user’s workload experiences only a marginal increase, even as substantial improvements in results are realized.