Flood Frequency Analysis Using the Bivariate Logistic Model with Non-Stationary Gumbel and GEV Marginals

1. Introduction

Floods constitute one of the most critical natural hazards, producing extensive economic impacts, infrastructure deterioration, and considerable risks to human life. Reliable flood frequency analysis (FFA) is therefore essential, as it provides the basis for estimating the magnitude of extreme events given a return period (T) and for supporting the safe design and management of hydraulic structures and flood-risk mitigation strategies. In this context, [1] provided a brief review of univariate and multivariate approaches to FFA.

Building on this perspective, research has increasingly turned to multivariate methods, which allow for a more complete representation of flood phenomena. Multivariate approaches to FFA have been applied globally. These studies often consider flood peaks, volumes, and durations (e.g., [1,2,3,4,5,6]). While many analyses incorporate two or more hydrological variables to capture joint flood behavior, some studies have focused on a single variable across multiple stations within a multivariate framework, as in [7,8,9,10,11]. These multivariate frameworks improve the realism of flood representation and support more reliable design criteria and risk management strategies in complex hydrological contexts.

Despite these advances, a key limitation of both univariate and multivariate approaches lies in their reliance on the stationarity assumption, which according to [12] it implies that the probability density function of any variable does not change over time. However, a growing body of evidence shows that this assumption is often violated due to both climatic and non-climatic drivers, challenging the applicability of conventional statistical approaches [12,13,14]. In response, a substantial volume of literature has emerged on non-stationary FFA, with applications spanning different regions and methodological frameworks [13,15,16,17]. Most non-stationary flood frequency studies rely on extreme value distributions such as GEV, Gumbel, and related families, often combined with Mann-Kendall tests. A dominant approach has been to use time as the main covariate, consistently identifying linear trends across regions. Beyond time, large-scale climate indices such as ENSO, PDO, NAO, and global mean temperature have been widely incorporated, showing strong influence on flood frequency dynamics.

Extending FFA to bivariate non-stationary marginals represents a promising avenue to better capture the joint dynamics of flood-related variables under changing conditions. In this study, we address this gap by applying a bivariate logistic model with non-stationary Gumbel and Generalized Extreme Value (GEV) marginals to flood frequency analysis.

2. Materials and Methods

According to [18], at-site frequency analysis based on annual maxima has a long-standing tradition in hydrology and remains one of the most widely adopted methods for quantile estimation [19]. Theoretical principles support the use of the GEV distribution, which represents the limiting distribution of block maxima. In addition to the GEV distribution, this study also applies the Gumbel distribution in its stationary form, as part of the modeling framework.

Gumbel and GEV distributions:

$$F\left(x\right)=e^{-e^{-\left({\displaystyle \frac{x-v}{\propto }}\right)}},$$

(1)

$$F\left(x\right)=e^{-{\left[1-\left({\displaystyle \frac{x-v}{\propto }}\right)\beta \right]}^{{\displaystyle \frac{1}{\beta }}}}.$$

(2)

2.1. Bivariate Logistic Model

[20] emphasize that the bivariate logistic model is more adaptable than the mixed model, as it accommodates a broader range of correlation coefficients and dependency indices. This flexibility makes it a preferable choice for addressing frequency analysis problems in hydrology, as also recommended by [21].

The model is defined as follows:

$$F\left(x,y,m\right)=exp\left\{-{\left[{\left(-lnF\left(x\right)\right)}^m+{\left(-lnF\left(y\right)\right)}^m\right]}^{1/m}\right\},$$

(3)

where $F\left(x\right)$ and $F\left(y\right)$ represent the marginal distribution functions, and $m (m>1)$ is the bivariate association parameter. When $m=1$, the bivariate distribution simplifies to the case of independence: $F\left(x,y,m\right)=F\left(x\right)F\left(y\right)$.

Using notation 1 for the Gumbel distribution and 2 for the GEV distribution, the possible combinations of bivariate extreme value distributions (BEV) are denoted as: BEV11, BEV12, BEV21, and BEV22 [21]. The probability distribution and density function for each BEV combination are shown in Appendix A along with the return level functions.

2.2. Estimation of Bivariate Parameters

Following [20], the likelihood function for $n$ random variables is defined as the joint density of those variables and serves as a function of the parameters. Given a random sample $\left(x_1,y_1\right)$ from a bivariate density, the corresponding likelihood function is formulated accordingly.

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

(4)

Due to unequal record lengths in the analyzed samples (Figure 1 and Figure 2), a flexible formulation is required to accommodate all possible data arrangements. This formulation is based on the generalization proposed by [22].

The likelihood function corresponding to the maximum sample arrangement is derived accordingly.

$$L\left(x,y,{\theta }_i\right)={\left[{\prod }_{i=1}^{n_1}f\left(p_i,{\theta }_i^{(1 )}\right)\right]}^{I_1}·{\left[{\prod }_{i=1}^{n_2}f\left(x_i,y_i,{\theta }_i^{\left(2 \right)}\right)\right]}^{I_2}·{\left[{\prod }_{i=1}^{n_3}f\left(q_i,{\theta }_i^{\left(3 \right)}\right)\right]}^{I_3},$$

(5)

where

$n_1,n_3=$ univariate lengths of record before and after the common period $n_2$

$n_2=$ common bivariate period

$\mathrm{p}=$ variable with univariate record before the common period $n_2$

$\mathrm{q}=$ variable with univariate record after the common period $n_2$

$\left(x,y\right)=$ variable with bivariate record during the common period $n_2$

$I_i=$ indicator number such that $I_i$ = 1 if $n_i>0$ or $I_i$ = 0 if $n_i=0$

${\theta }_i$ = parameter vector

Since the maximum of a function and its logarithm occur at the same point, and the logarithmic form of the maximum likelihood equation (5) is easier to manipulate, the logarithmic likelihood function is employed.

$$lnL\left(x,y,{\theta }_i\right)=I_1\left[{\sum }_{i=1}^{n}lnf\left(p_i,{\theta }_i^{(1 )}\right)\right]+I_2\left[{\sum }_{i=1}^{n}lnf\left(x_i,y_i,{\theta }_i^{(2 )}\right)\right]+I_3\left[{\sum }_{i=1}^{n}lnf\left(q_i,{\theta }_i^{(3 )}\right)\right].$$

(6)

Maximum likelihood estimators for the parameters of bivariate extreme value distributions are obtained by maximizing equation (6). Due to the complexity of the functions associated with bivariate probability densities, analytical solutions via differential calculus are not feasible. Instead, a direct search optimization process is required.

Given the characteristics of the function to be optimized, Rosenbrock’s restricted multivariate nonlinear optimization algorithm is recommended [23,24], as suggested by [21].

2.3. Non-Stationary GEV Models and Covariates

[25] gives an example of an annual maximum series $X_t$ to which the GEV distribution can be applied, even though its record has changed linearly over the observation time, but in other respects, the distribution does not change. Using the notation $GEV\left(\nu ,\alpha ,\beta \right)$ to denote the GEV distribution with parameters $\nu ,\alpha$ and $\beta$, it follows that a suitable model for $X_t$, the annual maximum value in a year t, can be:

$$X_t~GEV\left(\nu \left(t\right),\alpha ,\beta \right), \mathrm{where} \nu \left(t\right)={\nu }_0+{\nu }_{1}t,$$

(7)

for parameters ν₀ and ν₁. In this way, variations over time in the observed process are modeled as a linear trend in the location parameter of the appropriate extreme value model, which in this case is the GEV distribution. The parameter ν₁ corresponds to the annual rate of change in the annual maximum.

According to [25], a different situation that may arise is that the extreme behavior of one series is related to that of another variable, called a covariate. For example, the values of a series of maximums may appear larger in years when the Southern Oscillation Index (SOI) is higher. This suggests the following model for $X_t$:

$$X_t~GEV\left(\nu \left(t\right),\alpha ,\beta \right), \mathrm{where} \nu \left(t\right)={\nu }_0+{\nu }_{1}SOI\left(t\right),$$

(8)

and SOI(t) denotes the Southern Oscillation Index in year t.

There is a structural unit in all these examples. In each case, the extreme value parameters can be written in the form:

$$\theta \left(t\right)=h\left(X^{T}B\right),$$

(9)

where $\theta$ denotes $\nu ,\alpha$, or $\beta$, $h$ is a specific function, $B$ is a vector of parameters, and $X$ is a model vector.

One advantage of maximum likelihood over other parameter estimation techniques is its adaptability to changes in model structure. Take, for example, a non-stationary GEV model to describe the distribution $X_t$ for $t=1,\dots ,n$:

$$X_t~GVE\left(\nu \left(t\right),\alpha \left(t\right),\beta \left(t\right)\right),$$

(10)

where each parameter has an expression in terms of a vector of parameters and covariates of the type (9). Denoting the complete parameter vector as B, the likelihood is simply:

$$L\left(B\right)={\prod }_{t=1}^{n}g\left(x_t;\nu \left(t\right),\alpha \left(t\right),\beta \left(t\right)\right),$$

(11)

where $g\left(x_{\mathrm{t}};\nu \left(t\right),\alpha \left(t\right),\beta \left(t\right)\right)$ denotes the GEV density function with parameters $\nu \left(t\right),\alpha \left(t\right),\beta \left(t\right)$ evaluated at $x_{\mathrm{t}}$.

According to [17], trends or changes in location parameters, such as the mean or median, are the most visible and easiest to detect. Trends and changes in the scale parameter are less visible and detectable than in the location parameter. Changes in the shape parameter are not considered in the literature because they are not present in all distributions and can completely change the class of distributions.

Therefore, in this study, the non-stationary Gumbel and GEV distributions were applied, incorporating time and the Pacific Decadal Oscillation (PDO) and the Southern Oscillation (SOI) climate indices as covariates in the scale parameter, formulated as $\nu \left(c\right)={\nu }_0+{\nu }_{1}c$. The selection of these covariates is based on their documented influence in the study region [16].

2.4. Choice of Non-Stationary Model

[25] presents an alternative method for quantifying uncertainty in maximum likelihood estimators using the deviation function, defined as:

$$D=2\left\{{\mathcal{l}}_1\left(M_1\right)-{\mathcal{l}}_0\left(M_0\right)\right\},$$

(12)

Here, ${\mathcal{l}}_0\left(M_0\right)$ and ${\mathcal{l}}_1\left(M_1\right)$ represent the maximized log-likelihoods under models M₀ and M₁, respectively. A high value of $D$ indicates that model M₁ explains significantly more variation in the data than M₀, while a low value suggests that increasing model complexity does not enhance explanatory power.

Model M₀ is rejected at a significance level α if $D>c_{\alpha }$, where $c_{\alpha }$ is the (1 - α) quantile of the ${\mathcal{X}}_k^2$ distribution, and k is the difference in dimensionality between models M₁ and M₀. This provides a formal criterion for preferring model M₁ based on the magnitude of D.

To reduce statistical uncertainty and design safer hydraulic structures, this study proposes analyzing annual maximum discharges using a bivariate logistic model with non-stationary marginals. Figure 3 presents a schematic representation of the methodology followed, highlighting the main steps from data collection to flood frequency modeling.

3. Results

Eight stations with at least 30 years of records, located in a homogeneous region of Sinaloa in northern Mexico, were selected to apply the bivariate logistic model with non-stationary marginals for flood frequency analysis. Data were obtained from the official database of the National Commission of Water (CONAGUA) [26]. Table 1 lists the available length of record periods and statistical characteristics of these stations.

Outlier detection methods (box plot, Grubbs test, and three times the standard deviation test) were applied and independence was verified using the Anderson, Mantilla, and Amigó test, as required for random variable distributions. Homogeneity tests (Helmert, Student’s t-test, Cramer, Pettit, Standard Normal, Buishand, and Von Neuman), trend tests (Spearman and Mann-Kendall). Results are shown in Table 2.

Station 10086 exhibited an increasing trend, confirmed by graphical analysis (Figure 4) and statistical tests.

After verifying independence and randomness, the Gumbel (1) and GEV (2) distributions were fitted, followed by four bivariate marginal combinations (3): BEV 11, BEV 22, BEV 33, and BEV 44. The best combination for each station was selected based on the maximum likelihood logarithmic function.

Next, parameters for the non-stationary Gumbel (4) and GEV (5) distributions were estimated, followed by non-stationary bivariate combinations (6). Stations with D > 3.84 were considered, and the one with the highest log-likelihood was selected for both univariate and bivariate non-stationary models.

Table 3 presents the methodology results for station 10086, where the best-fitting distribution is the non-stationary BEV 21 in time, associated with station 10137. Table 4 shows estimated design events for various return periods, and Figure 5 compares the non-stationary BEV 21 distribution with the GEV distribution.

Finally, Table 5 lists the best distribution for each station in the homogeneous region, based on the maximum likelihood logarithmic function.

4. Discussion

Traditional stationary univariate models assume that the statistical properties of extreme events remain constant over time. While simple to apply, they neglect temporal changes driven by climate variability, land use, or anthropogenic impacts. As a result, return levels derived from stationary univariate approaches often misrepresent current and future risks [27].

Stationary bivariate models improve on this by accounting for dependence between two variables (for example, peak flow in two gauged stations), leading to more accurate joint return period estimates than univariate analyses [28]. However, they still impose the limitation of time-invariant marginals, which can bias estimates if the underlying processes evolve.

Traditional assumptions of stationarity are increasingly invalidated by anthropogenic influences such as urbanization, deforestation, and climate change, which significantly alter discharge records. So, non-stationary univariate models allow distribution parameters to vary with time or covariates like PDO and SOI, capturing observed or expected trends [25,27]. Yet, they continue to overlook the dependence between variables, which is critical for compound flood risk.

The most comprehensive approach is the non-stationary bivariate model, which jointly models dependence while allowing the marginals to evolve over time or with climatic covariates [29]. This reduces bias in joint tail estimation, improves accuracy of return levels, and provides a framework for assessing risks under future climate or land-use scenarios. In practice, non-stationary bivariate models deliver the most realistic and policy-relevant results, since they reflect evolving hazards and their interactions, offering clear advantages over both stationary and univariate alternatives.

In this study, the superior performance of the non-stationary bivariate model, according to the log-likelihood function, confirms that jointly modeling dependence while allowing the marginals to evolve over time yields more realistic and reliable return level estimates than stationary or univariate alternatives.

4.1. Performance of GEV and Gumbel Marginals

The best results were obtained when using GEV and Gumbel marginals to model the maximum discharges gauged at the two sites. The Generalized Extreme Value (GEV) distribution, with its shape parameter, offers flexibility to represent heavy, light, or bounded tails, whereas the Gumbel distribution, as a special case of the GEV, provides a simpler yet robust model for hydrological extremes [25,30]. This combination of marginals balances flexibility and parsimony, ensuring stable and interpretable parameter estimates.

At both gauging sites, the fitted marginals captured the statistical properties of annual maximum discharges. The GEV distribution was particularly suitable for the site with higher variability and occasional extreme peaks, while the Gumbel distribution provided a better fit at the site with more regular, less heavy-tailed extremes. This differentiation reflects the hydrological characteristics of the catchments: in one, extreme runoff events are more pronounced and require a flexible tail representation; in the other, the rainfall–runoff response is more consistent, making the Gumbel model sufficient.

The site-specific choice of marginals ensures that the joint bivariate model faithfully represents local hydrological variability and inter-site dependence. Importantly, this translates directly into practice: return period estimates used for infrastructure design are more reliable when tailored to the discharge behavior of each site, reducing the risk of under-design at the more variable site, and avoiding unnecessary over-design at the more stable site.

4.2. Case Study: Station 10086

The analysis of station 10086 revealed a statistically significant increasing trend, detected by both the Spearman rank correlation test and the Mann-Kendall test [31,32]. This trend is reflected in the positive value of the location parameter ($v_{11}$) of the non-stationary model. Figure 5 illustrates the importance of explicitly incorporating non-stationarity when modeling extremes.

For example, under stationary assumptions, the GEV distribution estimates a 50-year return flow of 772.66 m³/s. In contrast, the non-stationary BEV21 model projects a 50-year return flow of 861.00 m³/s for the year 2075 scenario. This difference implies that a discharge of 772.66 m³/s, considered a 50-year event under stationary assumptions, would be expected approximately every 30 years by 2075 when non-stationarity is accounted for. Such a misrepresentation has direct implications for hydraulic design, as ignoring trends may lead to underestimation of design floods and compromise the safety and resilience of hydraulic infrastructure [29].

4.3. Practical and Policy Implications

The adoption of non-stationary bivariate models has important implications for both engineering design and water policy. For engineers, it means that hydraulic structures—such as dams, levees, and urban drainage systems—can be designed with greater confidence and safety margins, reflecting realistic compound risks under changing conditions. For policymakers, non-stationary bivariate approaches support adaptive management strategies, allowing the testing of “what-if” scenarios under different climate or socio-economic pathways. This provides a sound basis for long-term investment decisions, disaster risk reduction strategies, and compliance with national and international climate adaptation frameworks [33].

Ultimately, the progression from stationary to non-stationary, and from univariate to bivariate analysis, represents a critical methodological advance in flood frequency analysis. It ensures that design standards and water policies evolve in step with changing climatic and hydrological realities, providing a more resilient and cost-effective basis for flood risk management.

5. Conclusions

This study confirms that non-stationary bivariate models provide the most reliable estimates of return levels, since they capture both the dependence between variables and the temporal evolution of extremes. Compared with stationary or univariate approaches, this framework reduces bias, improves accuracy, and supports scenario-based risk assessment under changing climatic and land-use conditions.

The use of GEV and Gumbel marginals yielded the best overall performance for the maximum discharges at the two gauging sites. The GEV distribution was particularly suited to the site with greater variability and occasional extremes, while the Gumbel distribution offered a parsimonious yet robust fit at the more stable site.

Importantly, the case of station 10086 illustrates how trend detection alters return level estimates. A discharge of 772.66 m³/s, considered a 50-year event under stationary assumptions, is projected by the non-stationary model to occur every ~30 years by 2075. Ignoring such non-stationarity could lead to systematic underestimation of design floods and compromise the safety of hydraulic infrastructure.

These findings underscore the importance of adopting flexible, non-stationary approaches for both engineering design and water policy. By tailoring flood frequency analysis to evolving hydrological regimes, decision-makers can avoid under- or over-design, improve safety margins, and align planning with climate adaptation frameworks.

Appendix A.1. Probability Distribution and Density Functions for the BEV11 Distribution

$$F\left(x_i,y_i;v_1,{\propto }_1,v_2,{\propto }_2,m\right)=exp\left\{-{\left[e^{-m\left({\displaystyle \frac{x-v_1}{{\propto }_1}}\right)}+e^{-m\left({\displaystyle \frac{y-v_2}{{\propto }_2}}\right)}\right]}^{1/m}\right\}$$

$$\begin{array}{c} f\left(x_i,y_i;v_1,{\propto }_1,v_2,{\propto }_2,m\right)=\left[{\displaystyle \frac{e^{-m\left(\frac{x-v_1}{{\propto }_1}\right)}·e^{-m\left(\frac{y-v_2}{{\propto }_2}\right)}}{{\propto }_1{\propto }_2}}\right]·{\left[e^{-m\left({\displaystyle \frac{x-v_1}{{\propto }_1}}\right)}+e^{-m\left({\displaystyle \frac{y-v_2}{{\propto }_2}}\right)}\right]}^{{\displaystyle \frac{1}{m}}-2}·\\ exp\left\{-{\left[e^{-m\left({\displaystyle \frac{x-v_1}{{\propto }_1}}\right)}+e^{-m\left({\displaystyle \frac{y-v_2}{{\propto }_2}}\right)}\right]}^{{\displaystyle \frac{1}{m}}}\right\}· \left\{\left(m-1\right)+{\left[e^{-m\left({\displaystyle \frac{x-v_1}{{\propto }_1}}\right)}+e^{-m\left({\displaystyle \frac{y-v_2}{{\propto }_2}}\right)}\right]}^{{\displaystyle \frac{1}{m}}}\right\}\end{array}$$

Appendix A.4. Log-Likelihood Function for the BEV11 Distribution

$$A_1=-\left({\displaystyle \frac{p_i-v_{\mathrm{p}0}-v_{\mathrm{p}1}c}{{\propto }_P}}\right),A_2=\mathrm{e}\mathrm{x}\mathrm{p}(A_1 )$$

$$B_1=-\left({\displaystyle \frac{q_i-v_{\mathrm{q}0}-v_{\mathrm{q}1}c}{{\propto }_q}}\right), B_2=\mathrm{e}\mathrm{x}\mathrm{p}(B_1 )$$

$$C_1=-\left({\displaystyle \frac{x_i-v_{10}-v_{11}c}{{\propto }_1}}\right), C_2=m{C}_1, { C}_3=\mathrm{e}\mathrm{x}\mathrm{p}(C_2 )$$

$$D_1=-\left({\displaystyle \frac{y_i-v_{20}-v_{21}c}{{\propto }_2}}\right), D_2=m{D}_1, D_3=\mathrm{e}\mathrm{x}\mathrm{p}(D_2 )$$

$$E_1=C_3+D_3, E_2={E_1}^{1/m}, E_3={E_1}^{\left[{\displaystyle \frac{1}{m}}-2\right]}$$

Appendix A.5. Log-Likelihood Function for the BEV12 Distribution

$$A_1=-\left({\displaystyle \frac{p_i-v_{\mathrm{p}0}-v_{\mathrm{p}1}c}{{\propto }_P}}\right),A_2=\mathrm{e}\mathrm{x}\mathrm{p}(A_1 )$$

$$B_1=-{\left(1-\left({\displaystyle \frac{q_i-v_{\mathrm{q}0}-v_{\mathrm{q}1}c}{{\propto }_q}}\right){\mathsf{\beta }}_q\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_q}}}, B_2=-\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{q_i-v_{\mathrm{q}0}-v_{\mathrm{q}1}c}{{\propto }_q}}\right){\mathsf{\beta }}_q\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_q}}-1}$$

$$C_1=-\left({\displaystyle \frac{x_i-v_{10}-v_{11}c}{{\propto }_1}}\right), C_2=m{C}_1, { C}_3=\mathrm{e}\mathrm{x}\mathrm{p}(C_2 )$$

$$D_2=\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{y_i-v_{20}-v_{21}c}{{\propto }_2}}\right){\mathsf{\beta }}_2\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_2}}-1}, D_3={\left(1-\left({\displaystyle \frac{y_i-v_{20}-v_{21}c}{{\propto }_2}}\right){\mathsf{\beta }}_2\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_2}}}$$

$$E_1=C_3+D_3, E_2={E_1}^{1/m}, E_3={E_1}^{\left[{\displaystyle \frac{1}{m}}-2\right]}$$

Appendix A.6. Log-Likelihood Function for the BEV12 Distribution

$$A_1=-{\left(1-\left({\displaystyle \frac{p_i-v_{\mathrm{p}0}-v_{\mathrm{p}1}c}{{\propto }_p}}\right){\mathsf{\beta }}_p\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_p}}}, A_2=-\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{p_i-v_{\mathrm{p}0}-v_{\mathrm{p}1}c}{{\propto }_p}}\right){\mathsf{\beta }}_p\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_p}}-1}$$

$$B_1=-{\left(1-\left({\displaystyle \frac{q_i-v_{\mathrm{q}0}-v_{\mathrm{q}1}c}{{\propto }_q}}\right){\mathsf{\beta }}_q\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_q}}}, B_2=-\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{q_i-v_{\mathrm{q}0}-v_{\mathrm{q}1}c}{{\propto }_q}}\right){\mathsf{\beta }}_q\right)}^{{\displaystyle \frac{1}{{\mathsf{\beta }}_q}}-1}$$

$$C_2=\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{x_i-v_{10}-v_{11}c}{{\propto }_1}}\right){\mathsf{\beta }}_1\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_1}}-1}, C_3={\left(1-\left({\displaystyle \frac{x_i-v_{10}-v_{11}c}{{\propto }_1}}\right){\mathsf{\beta }}_1\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_1}}}$$

$$D_2=\mathrm{l}\mathrm{n}{\left(1-\left({\displaystyle \frac{y_i-v_{20}-v_{21}c}{{\propto }_2}}\right){\mathsf{\beta }}_2\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_2}}-1}, D_3={\left(1-\left({\displaystyle \frac{y_i-v_{20}-v_{21}c}{{\propto }_2}}\right){\mathsf{\beta }}_2\right)}^{{\displaystyle \frac{m}{{\mathsf{\beta }}_2}}}$$

$$E_1=C_3+D_3, E_2={E_1}^{1/m}, E_3={E_1}^{\left[{\displaystyle \frac{1}{m}}-2\right]}$$

Appendix A.7. Return Levels for the Gumbel and GEV Distributions

Distribution	Function
Gumbel	${\widehat{X}}_T=\widehat{v}-\widehat{\propto }ln\left[-ln\left(F\left(x\right)\right)\right]$
GEV	${\widehat{X}}_T=\widehat{v}+{\displaystyle \frac{\widehat{\propto }}{\widehat{\beta }}}\left\{1-{\left[-ln\left(F\left(x\right)\right)\right]}^{\widehat{\beta }}\right\}$
Non-Stationary Gumbel	${\widehat{X}}_T={\widehat{v}}_0+{\widehat{v}}_{1}c-\widehat{\propto }ln\left[-ln\left(F\left(x\right)\right)\right]$
Non-Stationary GEV	${\widehat{X}}_T={\widehat{v}}_0+{\widehat{v}}_{1}c+{\displaystyle \frac{\widehat{\propto }}{\widehat{\beta }}}\left\{1-{\left[-ln\left(F\left(x\right)\right)\right]}^{\widehat{\beta }}\right\}$

where $F\left(x\right)=1-{\displaystyle \frac{1}{T}}$ stands for the non-exceedance probability and T represents the return period in years.