Determination of Hydrological Flood Hazards Thresholds and Flood Frequency Analysis: Case of Nokoue Lake Watershed

1. Introduction

Floods are currently the most frequent and damaging natural risk in West Africa. [1], [2,3,4]. They have harmful effects on the activities and populations living along the banks and involve significant security challenges for the most exposed areas. They are natural phenomena that are integral to the natural regime of water bodies (lakes) and watercourses (rivers), and protection against them requires prevention and forecasting [5,6,7]. Unlike the management approaches of the 1960s, current policies tend to better account for the significant role of floods and the means of managing the flood risk of a water body or a watercourse. Flood risk is particularly complex to understand due to its random nature associated with climate change, especially in highly developed and urbanized areas such as urban and peri-urban zones [8,9,10]. This issue is also present in the Nokoue lake watershed. The basin is indeed subject to high rainfall intensities that can be potentially devastating due to the rapid urban growth along its banks. Scientific literature has shown that, regardless of the nature of the floods, hydrological studies are often overlooked, and flood prevention structures are poorly designed due to the use of outdated empirical formulas [11,12,13]. Therefore, to reduce anxiety about the threat of flooding from Nokoue lake in Benin, the estimation of extreme water level is used in the context of public policy implementation for risk prevention or coastal management, particularly through the characterization of flood hazards. The purpose of these estimations is to provide a high level of safety in flood risk prevention. In this study, we apply probabilistic models to estimate extreme water level up to a 100-year return period for Nokoue lake. The estimates are made using an extreme value statistical analysis method. The choice of a 100-year return period is driven by political rather than scientific reasons. However, the estimation of extreme surcharges (quantiles) for return periods close to 100 years is useful for analyzing rare extreme events. The results of this study could provide valuable insights for evaluating the extreme scenarios referenced by public policies in flood risk prevention for Nokoue lake.

2. Materials and Methods

2.1. Study areas

Nokoue lake is located in southern Benin between 6°38' and 6°50' North latitude, and 2°35' and 2°55' West longitude. It extends between 150 and 170 km², respectively, during the low water period and the high-water period (Figure 1).

2.2. Statistic Description of Water Level Peaks Values

This study focuses on the statistical analysis of annual water level peaks of Nokoue lake, sourced from the Institute of Hydrology and Oceanology Research of Benin (IRHOB). An annual water level peak is defined as the maximum observed value within a year. These values were extracted from a daily data collecting covering the period from 1997 to 2022.

2.3. Standaedized Water Level Peaks Indices Calculation and Categorization of Flood Hazard

The identification of flood hazard thresholds is based on three steps: the calculation of standardized water level index, the representation of standardized index based on annual water level peaks, and the projection of normalized water level index on to the x-axis showing the annual water level peaks. The categorization of flood risk alert thresholds for Nokoue lake is based on the standardized water leval index inspired by McKee's study. By normalizing the annual water level peaks series of nokoue lake, these thresholds were determined [15,16,17]. The formula for the water level index is as follows.

$$I_h={\displaystyle \frac{X_i-\overline{X}}{\sigma }}$$

(1)

Where $I_h$ is the water level index, $X_i$ is the observed annual maximum value, $\overline{X}$ the mean of the annual peak water level values, $\sigma$ the standard deviation of the annual water level peaks values. The different risk categories such as: limited, moderate, significant, and critical are shown in Table 1.

2.4. Implementation of Frequency Models

Extreme water levels are estimated using a method of statistical fitting and extrapolation of extremes. Only the main points of the method are outlined here. For more information, please refer to [14]. The calculations were performed using the R environment because it is well known nowadays well. The methodology for establishing frequency fitting curves consists of three main steps: (i) statistical tests (stationarity, independence, and homogeneity), (ii) Selection and calculation of empirical probabilities of the extracted annual water level peak values and (iii) fitting to estimate parameters of the distribtuion and quantiles corresponding to several specified return periods associated with exceedance frequencies.

• Hypothesis testing.

The Mann-Kendall, Wald-Wolfowitz, and Wilcoxon tests were used respectively for stationarity, independence, and homogeneity. The p-value represents the risk of error if we consider that the null hypothesis $H_0$​ (the hypothesis that the sample is stationary, homogeneous, and independent) is not true. The maximum acceptable value for the risk of error is set at 5%. If the p-value is less than 5%, there is less than a one in five chance of being wrong in considering that the series of annual water level peaks is not independent, stationary, and homogeneous [19].

✓

Stationarity test

A series of random variables is said to be stationary if its statistical characteristics (mean, variance, or moments) remain invariant over time [19]. The Mann-Kendall test was used to verify stationarity, allowing us to test the following hypotheses:

-

H0: The statistical characteristics of the random variables are constant over time.

-

H1: The statistical characteristics of the random variables are not constant over time.

Given $n$ random variables $x_1,x_2,x_3,\dots ,x_n$ arranged in chronological order, the test statistic $S$ is expressed as follows:

$$S={\sum }_{i=1}^{n-1}{\sum }_{j=i+1}^{n}sign\left(X_j-X_i\right)$$

(2)

Under hypothesis $H_0$​, the statistic $S$asymptotically follows a normal distribution with a mean of zero and a variance $var\left(S\right)$. The closer the test statistic $S$ is to zero, the more the observations are considered to be stationary.

✓

Independence test

A data series is said to be independent if one data point is not influenced by the preceding data point [20]. For the implementation of this test, we used the Wald-Wolfowitz test. It allows us to verify the following hypotheses:

-

H0: the series is independent.

-

H1: the series is not independent.

When the series is sufficiently large, the Wald-Wolfowitz statistic $W$follows a normal distribution $U$with mean $\overline{w}$and variance $var\left(W\right)$. Let $x_1,x_2,x_3,\dots ,x_n$be random variables, then:

$$W={\sum }_{i=1}^{n-1}{x}_i\ast x_{i+1}+x_n\ast x_1$$

(3)

✓

Homogeneity

A sample of a random variable is said to be homogeneous when its data points come from the same distribution (collected under the same conditions). Several statistical tests (Mann-Whitney, Kruskal-Wallis, Wilcoxon tests, etc.) are used to ensure the homogeneity of a statistical series. The homogeneity test introduced by [20] was selected for this study. It is a non-parametric test that uses the ranks of observations instead of the actual values. Mathematically, the problem is formalized as follows: Given a series of observations of length $N$, from which two samples $X$and $Y$ are drawn, let $N_1$ and $N_2$ be the sizes of these samples, with $N=N_1+N_2$ and$N_1<N_2$. The values in the series are then ranked in ascending order. We then focus only on the rank of each element in the two samples within this series. If a value appears multiple times, it is assigned the corresponding average rank. Next, we calculate $W_x$ by the formula equation 4:

$$W_x={\sum }_{k=1}^{n1}rank\left(X\right)$$

(4)

The sum of the ranks of the elements from the first sample in the combined series.

• Selection and calculation of empirical probabilities of water level peaks values

The experimental probabilities associated with the observations were calculated using formula in the equation 5 with the Weibull formula, which aims to obtain unbiased exceedance probabilities for all distributions. [18].

$$p\left(x_i\right)={\displaystyle \frac{i}{N+1}}etT\left(x_i\right)={\displaystyle \frac{1}{1-p\left(x_i\right)}}$$

(5)

Where $p$ is the exceedance probability of the annual water level peak, it is the rank of the peak in the series, and $N$ is the size of the series consisting of the annual water level peaks.

• Fitting distributions to the sample of annual water level peaks values

A parametric distribution is fitted to the annual water level peaks. Adopting a distribution to study and describe annual water level peaks is undoubtedly the most critical step, introducing the greatest uncertainties [20,21]. It is wise to test other distributions within the asymptotic domain of extreme events. Various approaches have helped guide this choice, but unfortunately, no universal or foolproof method exists [22]. It is prudent to test other distribution belonging to the asymptotic domain of extreme events. Various approaches can help facilitate this choice, but unfortunately, there is no universal and infallible method. [23]. The annual water level peaks values were fitted to different probability distributions to determine the quantiles (estimated values assigned to events with a desired frequency or return period) for various return periods. We then selected the distributions that best fit the entire dataset. RStudio software was used for this part. In this study, the Generalized Extreme Value (GEV) distribution, the Gumbel distribution and the Generalized Pareto (GPA) distribution are all types of generalized extreme value that are often used to model extreme events, such as river or lake floods. A comparative study of the performance of these recommended distributions by [24] is the best approach for justifying the choice of a distribution. The linear moments method available in the (lmomco) package, based on negative logarithmic likelihood, was used for parameter estimation. The three distributions functions used in this article are as follows in these equation 6, 7 & 8:

$$\mathrm{G}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{l}:F\left(X\right)=e^{{-}^{e^{-\alpha (x-x0)}}}$$

(6)

Where $\alpha (x-x_0)$ for the Gumbel distribution, where $\alpha$ is the scale parameter and $x_0$ is the location parameter.

$$\mathrm{G}\mathrm{E}\mathrm{V}：\left\{\begin{array}{c}f\left(x\right)=e^{\left\{-{\left[1-k\left({\displaystyle \frac{x-x_0}{s}}\right)\right]}^{{\displaystyle \frac{1}{k}}}\right\}}fork\ne 0\\ f\left(x\right)=e^{-\left[-e^{\left(-{\displaystyle \frac{x-x_0}{s}}\right)}\right]}fork=0\end{array}\right\}$$

(7)

Where:

$k$: Positive, non-zero shape parameter

$x_0$​: Location parameter

$s$: Scale parameter

In the case where $k=0$, this distribution simplifies to the Gumbel distribution.

GPA : Let $X$ be a random variable with distribution function $F$ and $u$ be a threshold value. The random variable $Y=X–u$ pour $X>u$ follows the conditional distribution function:

$$\left\{\begin{array}{c}G\left(y\right)=G\left(x-\mu \right)={\displaystyle \frac{F\left(x\right)-F\left(\mu \right)}{1-F\left(x\right)}}withx>\mu \\ G\left(y\right)=1-{\left[1+k\left({\displaystyle \frac{y}{\sigma }}\right)\right]}^{-{\displaystyle \frac{1}{k}}}\end{array}\right\}$$

(8)

The linear moments used to estimate the parameters of the distribution in this study are linear combinations of weighted probability moments. The estimators are derived from solving a system of equations that equates the sample linear moments with those of the theoretical distribution to be fitted. For the linear moment method [24], which proposes numerous distributions ranging from 1 to 5 parameters, the method is formalized as follows:

Let $x$ be a random variable with distribution function $F$, and let $x_{1:n}$, $x_{2:n}$, $x_{3:n}$, …., $x_{n:n}$. ​ represent the order statistics for a sample of size $n$. The order statistics $\mathrm{1,2},3,\dots n$ for a sample of size $n$ are given in the equation 9:

$${\lambda }_r=r^{-1}{\sum }_{k=0}^{r-1}{\left(-1\right)}^k\left({\displaystyle \genfrac{}{}{0pt}{}{r-1}{k}}\right)E\left(x_{r-k:r}\right)$$

(9)

with $\left\{\begin{array}{c}r=\mathrm{1,2},3,\dots n\\ E\left(x_{j:r}\right)={\displaystyle \frac{r!}{\left(j-1\right)!\left(r-j\right)!}}\int {\left\{x\left(f\left(x\right)\right)\right\}}^{j-1}\ast \left\{{\left(1-F\left(x\right)\right)}^{r-1}\right\}dF\left(x\right)\end{array}\right.$Where $E$ representing the mean of the random variable $x$, $\lambda$ the moment and $r$ representing the order of moments.

2.5. Model Performance Metrics

The quality of the statistical extrapolation of extreme events is assessed using linear moments diagrams, Taylor diagrams, cumulative distribution functions, and the root mean square error (RMSE), as these methods are more practical and powerful compared to the χ² test, Bayesian Information Criterion (BIC), and Akaike Information Criterion (AIC).

• Root Mean Square Error criterion

The root mean square error (RMSE) is a method for objectively evaluating the performance of models. It provides a measure of the average magnitude of prediction errors, with lower values indicating better model accuracy. It is formulated as follows in the equation 10:

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N(O_n-M_n}{)}^2$$

(10)

• Linear moments diagram

par la formule de The linear moments diagram is based on the combination of skewness coefficients (${\tau }_3$) and kurtosis coefficients (${\tau }_4$) to graphically assess which distribution best fits the sample observations. Constructing the diagram requires knowledge of the function relating ${\tau }_3$ to ${\tau }_4$ ​ through the formula in the equation 11:

$$\tau 4={\sum }_{j=0}^8{a}_j\left({\tau }_3^2\right)$$

(11)

With a_j ​: polynomial approximation coefficient. The linear moments package in the R programming language was used to represent the kurtosis coefficients as a function of skewness coefficients and the experimental characteristic of the sample.

• Taylor diagram

It is a two-dimensional diagram that visualizes the relationship between observed and simulated data using three concise statistical parameters: the correlation coefficient $R^2$, the root mean square error $\left(RMSE\right)$ , and the standard deviation ($\sigma$). To create this, we used the plotrix package with its taylor.diagram function in R, which evaluates both the correlation and the RMSE between the empirical and theoretical distributions. Let $O_n$​ and $M_n$​ be the observed values and model predictions with respective means and standard deviations ${\mu }_1,{\sigma }_1$and ${\mu }_2,{\sigma }_2$​, where $n=\mathrm{1,2},3,\dots ,N$and N is the number of observations. The statistical parameters $R^2$ and $RMSE$ are given by the following equations 12 & 13:

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\sum }_{i=1}^N(O_n-M_n}{)}^2$$

(12)

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{\left(Q_{obs}^i-Q_{sim}^i\right)}^2}{{\sum }_{i=1}^n{\left(Q_{obs}^i-{\overline{Q}}_{obs}\right)}^2}}$$

(13)

3. Results

3.1. Description of the Annual Water Level Peaks Values

The statistical examination of the data series indicated that, over an 27-year period, the water level peaks of Nokoue lake ranged from 3.5 to 4.4 meters, observed in 1997 and 2022, respectively, with 3.95 meters as a mean annual water level peaks of Nokoue lake (Table 2, Figure 2).

3.2. Results of Standardized Water Level Index

The Figure 3 shows the different categories flood hazard based on the standardized water level peaks index and the interannual evolution of the annual water level peaks index in the Nokoue lake watershed. The analysis of Figure 3 reveals that the water level peaks index for Nokoue lake exhibit variability throughout the entire series. Additionally, a positive index characterizes wet years, ranging from 0.3 to 3.3 in the Nokoue lake watershed. Eight years are considered extremely wet across the watershed based on standardized water level index that exceeds +Similarly, three hydrological phases emerge from the examination of Figure The first phase corresponds to the period 1997-2001, characterized by a high frequency of positive indices; among the eight extremely wet years, three (1997, 1998, and 1999) fall within this phase. The second phase spans from 2002-2008, marked by a decline in annual peak water levels and a predominance of negative index (not covered in this study). The final phase covers the period from 2003 to 2022, characterized by a slight recovery in peak water level with positive index. This resurgence in water level peaks is accompanied by extreme rainfall events, which can cause floods and socioeconomic and environmental damages within the study area. This was the case in 2018, 2020, and 2022, during which the Nokoue lake watershed recorded positive index of +3.3.

Figure 5 shows the identified alert thresholds for water levels in the Lake Nokoué watershed.

The analysis of Figure 4 indicates that the water level thresholds classes are as follows: limited hazard classes between negative infinity and 3.94 m, moderate risk between 3.94 and 4.04 m, significant hazard classes between 4.04 and 4.14 m, critical hazard classes between 4.14 and 4.95 m. At the moderate risk threshold, flooding may occur if the watershed receives a certain amount of rainfall, with the extent of flooding being moderate. The significant and critical hazard thresholds are more prominent in the Nokoue lake watershed, helping to characterize hydrologically wet years when flood events are predominant as the watershed experiences rainfall. This is not without socioeconomic and ecological consequences in the study area.

3.3. Results of Hypothesis Tests

The results of statistical tests conducted on the extracted annual peak water level data indicate a rejection of the alternative hypothesis ($H_1$) in each test (Table 3). Given that if the p-value is greater than 0.05, the null hypothesis ($H_0$) is accepted, this outcome implies homogeneity, independence, and stationarity among the annual peak water levels.

• The hypothesis that the data series of annual peak water level is independent is accepted with a 95% confidence level. There is no correlation between the data in the series.
• The absolute value of the Mann-Kendall statistic $\left(\right|K\left|\right)$ is evaluated at 0.The hypothesis that there is no trend in data series is accepted at a 5% significance level.
• The absolute value of the Wilcoxon statistic $\left(\right|W\left|\right)$ is evaluated at 0.The mean of the two sub-samples (1997-2015 and 2016-2022) is statistically equal, meaning the series is homogeneous. Thus, the null hypothesis $H_0$ is accepted at a 5% significance level.

3.4. Results of Empirical Probability

The empirical probability density is composed of two phases: a rising phase from 0.3 to 1.2 and a declining phase from 1.25 to 0.The peak water level with the highest empirical probability densities are between 3.8 meters and 4.0 meters (Figure 5)

Figure 5. Empirical probability density of annual water level peaks.

Figure 5. Empirical probability density of annual water level peaks.

3.5. Results of the Fitting to Statistical Distributions

The L-moments estimation method applied to determine the parameters of the three best distributions such as Gumbel, GEV, and GPA fitted to the annual peak water level of Nokoue lake yielded the following results (Table 4, Figure 6, Figure 7, Figure 8, Figure 9) are presented below.

The plots of annual peak water level of Nokoue lake as a function of return period, created on a semi-logarithmic scale, do not perfectly follow the shapes of the various distribution curves. Figure 9 below shows the results of the fitting performed for each of the three best statistical distributions such as Gumbel, GEV, and GPA. The graphical analysis indicates that the scatter plot is best represented by the Gumbel distribution. The result in the root mean square error (RMSE) of the Gumbel distribution is estimated at 0.0724 (Figure 9). It is the smallest value of the observed root mean square error, indicating that the Gumbel distribution best represents the annual peak water level of Nokoue lake. It is important to note that the performance evaluation focused on the ability of different theoretical distributions to replicate the sample of peak water level values for Nokoue lake. For this purpose, cumulative distribution function curves, L-moment diagrams, and Taylor diagrams were constructed to identify the best-fitting distributions. The comparison of the cumulative distribution function of the Gumbel distribution with those of the GEV and GPA distributions fitted to the annual peak water level revealed different root mean square errors (RMSE). The Gumbel distribution performed best, with an RMSE of 0.0724, compared to RMSEs of 0.0754 for the GEV and 0.0761 for the GPA distributions, respectively (Figure 10a, Figure 10b). Recall that the L-moment diagram is a graphical representation of the coefficients of kurtosis against the coefficients of skewness. These were plotted to compare the distributions that align with the sample of annual peak water level for Nokoue lake. This representation showed that the L-moment diagrams of both the Gumbel and GEV distributions closely approximate the annual peak water levels of Nokoue lake (Figure 10d). The Taylor diagram was used to evaluate the root mean square error (RMSE), correlation coefficient, and standard deviation for the generalized extreme value (GEV) distribution applied to the sample of annual peak water level of Nokoue lake. The comparison of these coefficients revealed that the Gumbel model has the highest correlation coefficient, estimated at 97%, and the lowest values of RMSE and standard deviation ($\sigma$), estimated at 1 and 0.45, respectively (Figure 10c).

3.6. Results of the Water Level Peaks Estimates for the Gumbel, GEV, and GPA Distributions

For the preliminary determination of flood quantiles at extreme frequencies and the return period (RP) of the reference flood, the three best distributions were selected based on their superior performance in the model selection tests (Table 5, Table 6). The Gumbel distribution appears to fit the tail of the distribution better. For exceedance probabilities ranging from 0.85 to 0.01, the quantiles estimated with the Gumbel distribution closely approximate the empirical quantiles (Table 5). The quantiles associated with return periods (RP) of 2, 3, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, and 100 years are provided in Table The reference flood for the Nokoue lake station is the highest recorded quantile value, corresponding to 5.05 m/day with the Gumbel distribution, with an associated return period of 150 years. This implies that the reference flood has a 0.66% chance of occurring annually. The annual peak water level for the 100-year return period is 4.95 m/day. The box plot (Figure 11) shows that most of the annual peak water level of Nokoue lake estimated with the Gumbel distribution are below the mean.

4. Discussion

The standardized water level index for Nokoue lake, ranging from $-1.5$ to $+3.9$, allow for the categorization of hydrological thresholds, corresponding respectively to the limited, moderate, significant, and critical hazard classes [25,26,27]. This categorization of hydrological thresholds helps define the extent of flooding in the Nokoue lake watershed. Flood hazard water level thresholds vary according to the risk level, as shown in Figure The analysis of Figure 3 indicates that the periods $1997-2001$ and $2009-2022$ were characterized by a high frequency of positive index in the Nokoue lake watershed. This suggests that the study area experienced hydrologically wet conditions, marked by floods during these years. These periods of significant water level increases in the Nokoue lake watershed result from the inflows from the lake's tributaries, combined with excess rainfall over the corresponding periods. This return to a wet state aligns with previous studies conducted in West Africa, specifically in Benin by [25,28], [51] and in Burkina Faso by [26,27] and elsewhere [29,30], [52,53]. The results of the statistical fitting show that for the sample of annual water level peaks of Nokoue lake, the Gumbel, GEV, and GPA distributions are the three best fits for the data. The fitting curves for these three distributions are quite similar. The plot of annual water level peaks against the logarithm of the return period displays a concave shape with indefinite growth. This concavity is dictated by the sign of the shape coefficient for each distribution, which influences the curvature in the distribution tails. Beyond the 10-year return period, the Gumbel distribution shows higher concavity compared to GEV and GPA. From an asymptotic perspective, the Gumbel distribution provided the best fit for the sample, making it the selected model for annual water level peaks. This could be explained by its asymptotic behavior [31]. This could be partly explained by the fact that the Gumbel distribution provides acceptable results when using the robust L-moment estimation method [32,33,34], [35,36,37]. Additionally, it is supported by extreme value theory. The Fisher-Tippett-Gnedenko theorem states that for a series of independent and identically distributed observations following a common distribution $F$, the limit of their linearly normalized maxima as $N$ approaches infinity converges toward a Gumbel distribution [38,39], [40,41,42], who studied around thirty hydrometric stations in Wallonia, also demonstrated that the Gumbel distribution was suitable for determining flood flows in 95% of basins when using the annual peak flow method. The authors of this reference [43,44,45,46], [47] demonstrated that the GPA distribution was much more suitable for a series composed of values exceeding a threshold. It confirms that the tail of the GPA distribution is thicker as the value of the parameter $k$ increases. Our results corroborate with those of [48,49], [50]. These authors, in their study, demonstrated that the Gumbel distribution was more efficient than the GEV and GPA distributions. They recommended that the estimation of parameters using the maximum likelihood method for the GEV distribution and the generalized Pareto distribution (GPA) can be carried out very effectively and accurately using a global optimization tool that can bypass various local optima.

5. Conclusions

The Gumbel model correctly reproduces the curves of peak flood water heights for Nokoue lake. The return period of the largest flood experienced by Nokoue lake is 150 years. The Gumbel model will be used for the preliminary determination of the peak water heights related to floods in the Nokoue lake watershed. The main limitation of this work lies in the choice of probability distributions and the method of parameter estimation. Indeed, there is no universal and infallible method for choosing the distributions suited to different situations. However, this case study of the Nokoue lake watershed serves as a basis for all flood prevention structures in the Nokoue lake basin and the determination of alert thresholds.