Application of the Bivariate Exponentiated Gumbel Distribution for Extreme Rainfall Frequency Analysis in Contrasting Climates of Mexico

1. Introduction

[1] highlighted that variations in the hydrological cycle have profound impacts on human activities, including an increased risk of flooding as well as the occurrence and severity of droughts. Floods are among the most damaging natural hazards worldwide, and their impacts have intensified due to factors such as irregular human settlements near rivers, deforestation, and continuous land-use changes.

Mexico is especially vulnerable to hydrometeorological events that cause widespread damage across its territory. Between 2007 and 2020, heavy rains, floods, and tropical cyclones affected 27.8 million people, resulting in 1,173 fatalities and economic losses exceeding 24.6 billion USD [2]. Analyses of rainfall patterns in Mexico [3,4] indicate that floods caused by extreme rainfall events are likely to become more frequent and intense in the future. Given the potential consequences of underestimating or overestimating rainfall quantiles—ranging from hydraulic structure failure to unnecessarily inflated construction costs—reliable statistical tools are essential for estimating their magnitude and frequency.

Numerous at-site and regional rainfall frequency analyses have been conducted worldwide, and results consistently demonstrate that no single distribution universally fits annual maximum daily rainfall (AMDR) data. Although the Gumbel distribution (EVI) is widely applied, studies such as [5] in Egypt have shown that it is not always appropriate, while in other regions it has been identified as the most suitable model [6,7,8,9,10,11,12,13,14,15].

In Italy, [16] analyzed AMDR records from 297 gauging stations and found that the Fréchet distribution (EVIII) provided the best fit. Their study also revealed significant differences in return periods and return levels compared to results from EVI, reversed Weibull (EVII), Pareto, Lognormal, and Gamma distributions, echoing earlier findings by [17,18]. The EVI, EVII, and EVIII distributions are all special cases of the Generalized Extreme Value (GEV) distribution, which has often been identified as the best-fitting model for AMDR samples in multiple regions [19,20,21,22,23,24,25,26,27,28]. Similarly, the Log-Pearson Type III (LPIII) distribution has also proven effective in modeling AMDR series [29,30,31,32,33,34].

In recent years, the Exponentiated Gumbel (EG) distribution was proposed by [35], who demonstrated its improved performance in climate modeling compared to the EVI and GEV distributions. [36] applied EG, exponentiated Weibull (EW), exponentiated Fréchet (EF), and their mixed forms to AMDR series from 19 Mexican stations, showing that exponentiated and mixed exponentiated distributions provide flexible and reliable alternatives for modeling extreme rainfall. More recently, [37] compared the Gumbel and EG distributions, concluding—based on AIC and BIC criteria—that EG offered a more flexible fit for hydrological datasets, corroborating earlier findings by [35].

Regional frequency analysis has also emerged as a practical approach to reduce uncertainties associated with short or incomplete rainfall records at gauged sites. In this context, multivariate joint estimation models have proven valuable, as they enhance the estimation of marginal distribution parameters and improve regional at-site estimates of return levels by incorporating information from neighboring sites within homogeneous regions. Bivariate approaches have shown promise in flood frequency analysis [38,39,40,41,42,43,44,45,46,47].

In this study, a bivariate distribution with Exponentiated Gumbel marginals is proposed to improve the estimation of marginal parameters and corresponding quantiles. The performance of this model is compared with that of univariate probability distributions, namely GEV, EVI, and EG.

2. Materials and Methods

2.1. Study Area

Mexico is characterized by a summer–rainy and winter–dry regime, except in the northwest. Annual precipitation ranges from <500 mm in the north and northwest to >2,000 mm in the humid south and southeast. Two states that illustrate this climatic contrast are Coahuila, in northern Mexico, and Tabasco, in the southeast (Figure 1).

Coahuila (151,563 km²) has a semi-warm summer and cold winter climate, with a mean annual temperature of 20 °C (min. 4 °C, max. 30 °C). Rainfall is scarce, averaging ~400 mm yr⁻¹, concentrated in summer. Maximum daily rainfall varies between 11.3 and 453 mm. Recent extreme events have produced floods affecting Torreon, Monclova, and Piedras Negras.

Tabasco (25,267 km²) is predominantly warm-humid, with abundant summer rains (75.97%), year-round humid conditions (19.64%), and sub-humid summer rains (4.39%). The mean annual temperature is 27 °C (min. 18.5 °C, max. 36 °C). Rainfall occurs throughout the year, peaks from June to October, and averages 2,550 mm yr⁻¹. Maximum daily rainfall ranges from 38.8 to 816.9 mm.

The dataset consists of records from 106 rain gauge stations in Coahuila and 75 stations in Tabasco, obtained from [48]. For each station, annual maximum daily rainfall values were analyzed, with detailed station characteristics provided in Table 1 and Table 2.

2.2. Delineation of Homogeneous Regions

The joint parameter estimation model requires that all samples belong to the same homogeneous region. In this study, regions were delineated using the coefficient of variation (Cv) as the measure of comparability among stations.

In the first stage, three initial zones were identified by setting confidence limits defined as the mean Cv ± one standard deviation. At this stage, most stations were concentrated in region 2 (Figure 2), leaving only the stations in regions 1 and 3 distinctly separated.

In the second stage, the stations within region 2 were reclassified by recalculating new confidence limits and repeating the same procedure. This subdivision resulted in four regions (1, 2, 4, and 5) (Figure 3).

The process was repeated iteratively, each time focusing on the stations remaining in the central region (closer to the mean Cv value). In this case, region 3 was progressively subdivided until a sufficient distribution of stations across regions was achieved. After the final iteration, seven homogeneous regions were identified (Figure 4).

It is important to note that the total number of stages required depends on the number of stations analyzed, since the central group—located around the mean Cv—will ultimately determine how many subdivisions are necessary to achieve adequately homogeneous regions.

2.3. Univariate Distributions

a. GEV distribution

The cumulative distribution function is [49]:

$$F\left(x\right)=exp\left\{-{\left[1-\left({\displaystyle \frac{x-\epsilon }{\lambda }}\right)\beta \right]}^{{\displaystyle \frac{1}{\beta }}}\right\}$$

(1)

where $\epsilon ,\lambda ,\beta$ are the location, scale, and shape parameters, respectively.

And

$$f\left(x\right)={\displaystyle \frac{1}{\lambda }}exp\left\{-{\left[1-\left({\displaystyle \frac{x-\epsilon }{\lambda }}\right)\beta \right]}^{{\displaystyle \frac{1}{\beta }}}\right\}{\left[1-\left({\displaystyle \frac{x-\epsilon }{\lambda }}\right)\beta \right]}^{{\displaystyle \frac{1}{\beta }}-1}$$

(2)

The log-likelihood for the GEV distribution is:

$$\ln L\left(x;\epsilon ,\lambda ,\beta \right)=-n\ln \lambda +{{\displaystyle \sum }}_{i=1}^n\left\{-{\left[1-\left({\displaystyle \frac{x_i-\epsilon }{\lambda }}\right)\beta \right]}^{{\displaystyle \frac{1}{\beta }}}\right\}+\left({\displaystyle \frac{1}{\beta }}-1\right){{\displaystyle \sum }}_{i=1}^n\ln \left[1-\left({\displaystyle \frac{x_i-\epsilon }{\lambda }}\right)\beta \right]$$

(3)

The quantile ${\widehat{x}}_T$ corresponding to a given return period T is:

$${\widehat{x}}_T=\widehat{\epsilon }+\widehat{\lambda }\left\{{\displaystyle \frac{1-{\left[-\ln \left(1-\frac{1}{T}\right)\right]}^{\beta }}{\widehat{\beta }}}\right\}$$

(4)

b. EVI (Gumbel) distribution

$$F\left(x\right)=exp\left\{-exp\left[-\left({\displaystyle \frac{x-\nu }{\alpha }}\right)\right]\right\}$$

(5)

$$f\left(x\right)={\displaystyle \frac{1}{\alpha }}exp\left\{-exp\left[-\left({\displaystyle \frac{x-\nu }{\alpha }}\right)\right]\right\}exp\left[-\left({\displaystyle \frac{x-\nu }{\alpha }}\right)\right]$$

(6)

where $\nu ,\alpha$ are the location and scale parameters $\alpha >0$

$$\ln L\left(x;\nu ,\alpha \right)=-n\ln \alpha -{{\displaystyle \sum }}_{i=1}^{n}exp\left[-\left({\displaystyle \frac{x-\nu }{\alpha }}\right)\right]+{{\displaystyle \sum }}_{i=1}^n\left[-\left({\displaystyle \frac{x_i-\nu }{\alpha }}\right)\right]$$

(7)

$${\widehat{x}}_T=\widehat{\nu }-\widehat{\alpha }\ln \left[-\ln \left(1-{\displaystyle \frac{1}{T}}\right)\right]$$

(8)

c. Exponentiated Gumbel distribution (EG) [50].

$$F\left(x\right)=1-{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\right\}\right\}}^{\kappa }$$

(9)

$$f\left(x\right)={\displaystyle \frac{\kappa }{\sigma }}{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\right\}\right\}}^{\kappa -1}exp\left\{-exp\left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]\right\}exp\left[-\left({\displaystyle \frac{x-\mu }{\sigma }}\right)\right]$$

(10)

$$\begin{array}{l}\ln L\left(x;\mu ,\sigma ,\kappa \right)=n\ln \kappa -n\ln \sigma +\left(\kappa -1\right){{\displaystyle \sum }}_{i=1}^n\ln \left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)\right]\right\}\right\}-\\ {{\displaystyle \sum }}_{i=1}^{n}exp\left[-\left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)\right]+{{\displaystyle \sum }}_{i=1}^n\left[-\left({\displaystyle \frac{x_i-\mu }{\sigma }}\right)\right]\end{array}$$

(11)

$${\widehat{x}}_T=\widehat{\mu }-\widehat{\sigma }ln\left[-ln{\left(1-{\displaystyle \frac{1}{T}}\right)}^{{\displaystyle \frac{1}{\widehat{\kappa }}}}\right]$$

(12)

2.4. Biexponentiated Gumbel Distribution (BEG)

As already mentioned, multivariate extreme value distributions have shown to be a reliable option for fitting hydrological variables. [51] compared the mixed and logistic bivariate models and concluded that the logistic model (LM) is the best option in flood frequency analysis.

The LM is [52]:

$${\left[-\log F\left(x,y,m\right)\right]}^m={\left[-\log F\left(x\right)\right]}^m+{\left[-\log F\left(y\right)\right]}^m$$

(13)

where x and y denote the extreme events gauged in a couple of neighboring sites.

and

$$f\left(x,y,m\right)={\displaystyle \frac{{\partial }^{2}F\left(x,y,m\right)}{\partial x\partial y}};f\left(x,y,m\right)\ge 0$$

(14)

The bivariate likelihood function is:

$$L\left(x,y,\theta \right)={{\displaystyle \prod }}_{i=1}^{n}f\left(x_i,y_i,\theta \right)$$

(15)

The log-likelihood function is:

$$\ln L\left(x,y,\theta \right)={{\displaystyle \sum }}_{i=1}^n\ln f\left(x_i,y_i,\theta \right)$$

(16)

The expression (16) is used when both samples share a common length of record, however, the general form for considering different sizes is [53]:

$$\ln L\left(x,y,\theta \right)=I_1\left[{{\displaystyle \sum }}_{i=1}^{n_1}\ln f\left(r_i,\theta \right)\right]+I_2\left[{{\displaystyle \sum }}_{i=1}^{n_2}\ln f\left(x_i,y_i,\theta \right)\right]+I_3\left[{{\displaystyle \sum }}_{i=1}^{n_3}\ln f\left(s_i,\theta \right)\right]$$

(17)

Marginals F(x) and F(y) are proposed to be EG distributions. The BEG log-likelihood function to maximize in the parameter estimation procedure is:

(18)

where x and y are the variables in the common period (CP) with length $n_2$, r denotes to the variable x or y with length $n_1$before the CP, s is the variable x or y with length $n_3$ after the CP; $I_j=1$if $n_j>0\mathrm{or}{I}_j=0$if $n_j=0$ for j = 1, 2, 3; p =1 and q = 1 if r or s represent to variable x; p = 2 and q = 2 if r or s characterize to variable y, and

$$F\left(x\right)=1-{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{x_i-{\mu }_1}{{\sigma }_1}}\right)\right]\right\}\right\}}^{{\kappa }_1}$$

(19)

$$f\left(x\right)={\displaystyle \frac{{\kappa }_1}{{\sigma }_1}}{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{x_i-{\mu }_1}{{\sigma }_1}}\right)\right]\right\}\right\}}^{{\kappa }_1-1}exp\left\{-exp\left[-\left({\displaystyle \frac{x_i-{\mu }_1}{{\sigma }_1}}\right)\right]\right\}exp\left[-\left({\displaystyle \frac{x_i-{\mu }_1}{{\sigma }_1}}\right)\right]$$

(20)

$$F\left(y\right)=1-{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{y_i-{\mu }_2}{{\sigma }_2}}\right)\right]\right\}\right\}}^{{\kappa }_2}$$

(21)

$$f\left(y\right)={\displaystyle \frac{{\kappa }_2}{{\sigma }_2}}{\left\{1-exp\left\{-exp\left[-\left({\displaystyle \frac{y_i-{\mu }_2}{{\sigma }_2}}\right)\right]\right\}\right\}}^{{\kappa }_2-1}exp\left\{-exp\left[-\left({\displaystyle \frac{y_i-{\mu }_2}{{\sigma }_2}}\right)\right]\right\}exp\left[-\left({\displaystyle \frac{y_i-{\mu }_2}{{\sigma }_2}}\right)\right]$$

(22)

The Rosenbrock optimization algorithm for constrained variables [54] was selected for estimating the univariate and bivariate parameters by the direct maximization of equations (3), (7), (11), and (18).

2.5. Selection of Best Fit

The selection of fit between the Empirical and Theoretical distribution of the AMDR was based on the AICc and BIC goodness of fit tests.

The AIC was proposed by Akaike [55]:

$$\mathrm{AIC}=-2\ln L\left(\widehat{\theta }\right)+2p$$

(23)

while the BIC is [56]:

$$\mathrm{BIC}=-2\ln L\left(\widehat{\theta }\right)+p\ln \left(n\right)$$

(24)

where $\ln L\left(\widehat{\theta }\right)$ represents the log-likelihood of empirical distribution, p the number of maximum-likelihood estimates of the parameters, and n the length of record.

For small n, a corrected version of AIC is proposed as [57]:

$$\mathrm{AICc}=\mathrm{AIC}+{\displaystyle \frac{2p\left(p+1\right)}{n-p-1}}$$

(25)

Distribution having least value of AICc, or BIC is considered as best model.

2.6. Reliability of Estimated Quantiles

It is very important to evaluate whether the BEG distribution provides more accurate and reliable quantile estimates than traditional univariate approaches. This evaluation is essential in hydrological frequency analysis, since underestimation of quantiles can increase the risk of hydraulic structure failure, while overestimation may result in unnecessarily high construction costs. To assess reliability, quantile estimates obtained from the BEG model were compared with those from univariate distributions using two statistical criteria: bias (BIAS) and mean squared error (MSE). This framework allowed us to determine not only the accuracy of the estimated quantiles but also the extent to which the joint estimation procedure enhances the transfer of information across sites, thereby improving the robustness of extreme rainfall frequency analysis.

To formalize this evaluation, the reliability of estimated quantiles was expressed in terms of bias and mean squared error, which are defined as follows:

Let $\widehat{\psi }$ be the quantile to be computed by using the BEG distribution:

$$\mathrm{BIAS}=m\left(\widehat{\psi }\right)-\psi$$

(26)

$$\mathrm{MSE}=S^2\left(\widehat{\psi }\right)+{\left[m\left(\widehat{\psi }\right)-\psi \right]}^2$$

(27)

For the number of simulated samples “ $n_s$”

$$m\left(\widehat{\psi }\right)={\displaystyle \frac{1}{n_s}}{{\displaystyle \sum }}_{i=1}^{n_s}{\widehat{\psi }}_i$$

(28)

And

$$S^2\left(\widehat{\psi }\right)={\displaystyle \frac{1}{n_s}}{{\displaystyle \sum }}_{i=1}^{n_s}{\left[{\widehat{\psi }}_i-m\left(\widehat{\psi }\right)\right]}^2$$

(29)

When estimating quantiles is desirable to have unbiased and minimum MSE estimators.

3. Results

A procedure of data generation was performed to show if the quantiles obtained through of the bivariate joint estimation of parameters are more reliable with those obtained by its univariate counterpart.

For the EG distribution, data were generated using population parameters ${\mu }_1=15.8$, ${\sigma }_1=2.5$ and ${\kappa }_1=1.3$ with sample sizes n = 10, 20 and 50. A total of 1,000 simulated samples were considered for each n.

For the BEG distribution, quantiles were obtained by combining samples with sizes n₁-n₂: 10-10, 10-20, 20-20, 20-50, 50-50 and 50-100. Comparisons were performed for non-exceedance probabilities of 0.50, 0.80, 0.90. 095 0.98 and 0.99. The associated site has population parameters ${\mu }_2=18.7$, ${\sigma }_2=3.0$ and ${\kappa }_2=1.1$.

Results (Table 3 and Table 4) indicate that as $n_2$increased relative to $n_1$, both BIAS and MSE of the shorter series decreased. This demonstrates an effective transfer of information when parameters are jointly estimated, supporting the conclusion that quantiles computed using the BEG distribution are more reliable than those obtained from the univariate case.

The rain gauge stations in Coahuila and Tabasco were classified into homogeneous regions by applying the delineation procedure previously described. This methodological approach yielded nine regions in Coahuila (Table 5) and seven regions in Tabasco (Table 6), ensuring consistency with the adopted regionalization framework.

The analysis of the Coahuila dataset considered all possible pairwise station combinations within each homogeneous region to evaluate the dependence structure of extreme rainfall. Table 7 and Table 8 illustrate representative examples for selected stations, showing the BEG model combinations applied (Table 7) and the corresponding return level estimates (Table 8). When extended to the complete set of stations, the procedure enabled a systematic assessment of return levels using the BEG framework. For comparison, alternative univariate models—Gumbel, Generalized Extreme Value (GEV), and Exponentiated Gumbel—were also fitted to the same stations, and the resulting return levels were ranked according to goodness-of-fit statistics, including the Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) (Table 9). Notably, in 65% of the stations in Coahuila the best fit was achieved with the BEG distribution, highlighting its reliability for modeling extreme rainfall in this region.

A similar procedure was applied in Tabasco, where all possible station pairs within each homogeneous region were analyzed. Table 10 presents the parameters obtained for the best BEG combination at each station. Return levels derived from the BEG distribution were then compared with those obtained from alternative univariate models (Gumbel, GEV, and Exponentiated Gumbel), and the results were ranked according to statistical criteria (AIC and BIC) (Table 11). In this case, 81% of the stations achieved the best fit with the BEG distribution, confirming the robustness of this approach under the markedly different climatic conditions of southeastern Mexico.

Together, these results demonstrate that the BEG distribution consistently outperforms classical models, providing the best fit in both the arid to semi-arid conditions of Coahuila and the humid tropical environment of Tabasco. This consistency across regions and evaluation metrics underscores the versatility of the BEG framework for regional frequency analysis in contrasting hydroclimatic settings.

4. Discussion

The results obtained in Coahuila and Tabasco demonstrate the advantages of adopting flexible probability models such as the Bivariate Exponentiated Gumbel (BEG) distribution for extreme rainfall analysis. In both regions, the BEG framework provided the best fit for a majority of stations—65% in Coahuila and 81% in Tabasco—when compared to the classical Gumbel, Generalized Extreme Value (GEV), and Exponentiated Gumbel distributions. This superior performance, confirmed by goodness-of-fit criteria including AIC and BIC, highlights the ability of the BEG model to capture the dependence structure of extreme events more effectively than univariate alternatives.

The simulation experiments reinforce this conclusion. By jointly estimating parameters across paired samples, the BEG approach significantly reduced bias and mean squared error in quantile estimation, particularly as one sample length increased relative to the other. This result demonstrates the effective transfer of information between series, a feature that is especially relevant in hydrological contexts where records are often short, fragmented, or incomplete. In contrast, univariate models produced higher estimation errors under identical conditions, confirming the greater robustness of BEG.

An important aspect of these findings is the consistency of the BEG performance across contrasting climatic contexts. In Coahuila, a predominantly arid to semi-arid region, rainfall extremes are generally short-lived and spatially heterogeneous, which complicates regionalization and frequency analysis. Conversely, Tabasco is located in a humid tropical environment where rainfall extremes are often more spatially extensive and influenced by large-scale atmospheric dynamics. Despite these marked hydroclimatic differences, the BEG distribution consistently outperformed traditional models, underscoring its robustness and adaptability.

The systematic evaluation of all possible station pairs within homogeneous regions further strengthens the validity of the results. This approach ensures that spatial dependence is explicitly considered, rather than assuming independence between stations—a limitation common in univariate frequency analysis. The strong performance of the BEG distribution suggests that bivariate or multivariate models may be more appropriate for regions with high spatial variability, particularly where water resource planning and hydraulic design require accurate estimation of joint extremes.

From a practical perspective, the improved fit obtained with the BEG distribution has direct implications for risk management and infrastructure design. Underestimation of return levels, especially for long return periods, can lead to inadequate sizing of hydraulic structures and increased vulnerability to extreme events. By providing more reliable estimates, the BEG framework contributes to reducing uncertainty in hydrological design, supporting more resilient adaptation strategies in the face of climate variability and change.

Finally, these results align with recent studies emphasizing the importance of moving beyond stationary univariate models in hydrology. The incorporation of flexible, bivariate approaches not only improves statistical performance but also offers a more realistic representation of rainfall extremes, particularly in regions with complex climatic dynamics. The outcomes from Coahuila and Tabasco therefore provide empirical evidence supporting the broader adoption of BEG-based regional frequency analysis in Mexico and comparable hydroclimatic contexts.

5. Conclusions

This study applied the Bivariate Exponentiated Gumbel (BEG) distribution to extreme rainfall analysis in Coahuila and Tabasco, two Mexican states with contrasting climatic conditions. The main conclusions are as follows:

Across the stations analyzed, the BEG distribution provided the best fit in 65% of the cases in Coahuila and 81% in Tabasco, outperforming classical alternatives such as the Gumbel, Generalized Extreme Value (GEV), and Exponentiated Gumbel distributions. This performance was consistently confirmed by statistical indicators including AIC and BIC.

Simulation experiments showed that BEG reduced bias and mean squared error when estimating quantiles, especially for short or heterogeneous samples. The joint estimation of parameters allowed effective transfer of information, resulting in more reliable quantiles than those obtained from univariate approaches.

The BEG model demonstrated strong adaptability in both arid to semi-arid conditions (Coahuila) and humid tropical conditions (Tabasco), underscoring its versatility as a tool for regional frequency analysis of extreme rainfall.

The systematic evaluation of all possible pairwise station combinations within homogeneous regions captured spatial dependence more effectively than univariate models. This methodological strength enhances the reliability of rainfall return level estimates, particularly for long return periods.

The improved accuracy of return level estimates obtained with the BEG distribution can reduce underestimation of design events, thereby contributing to safer and more resilient hydraulic infrastructure, as well as supporting climate adaptation planning in water management.

In summary, the BEG distribution proves to be a reliable and flexible alternative for regional extreme rainfall analysis in Mexico. Its demonstrated robustness suggests that it can be extended to other regions with similar hydroclimatic variability, providing a valuable tool for both scientific research and applied hydrological practice. Future research should explore the application of the BEG model under non-stationary conditions, test its performance in multivariate settings that include variables such as streamflow or temperature, and assess its potential for integration into climate change impact studies.