Autocorrelation Ratio as a Measure of Inertia for the Classification of Extreme Events

Abstract

One of the problems in modern hydrology concerns the estimation of design events at sites with scarce or null data. This is a challenge for developing countries because they do not have monitoring networks as extensive and reliable as in developed nations. This situation has caused states in the Latin American and Caribbean region, in particular, to rely on hydrological regionalization techniques. These procedures implement clustering algorithms in combination with aggregation rules and metric distances to generate homogeneous groups from which hydrological information can be transferred. In addition, it has been proven that the analysis of spatial variables is sensitive to the magnitudes of extreme events; therefore, a mathematical formulation that adopts this fact into consideration must be included. For this purpose, the autocorrelation distance of the daily rainfall data series is proposed as an estimator of temporal variability. The fit parameters of the mixed Poisson-exponential probability distribution are operated as estimators of spatial variability. These spatio-temporal conditions are combined to obtain a mathematical relation of the autocorrelation as a measure of inertia for the classification of extreme events. This procedure is applied to Hydrologic Region 10 in northwestern Mexico from daily rainfall records. This zone has already been explored in terms of its regional homogeneity, which allows validating the results obtained.

1. Introduction

Hydrologic Region 10 (RH10) located in northwestern Mexico is one of the most critical regions of the country because it is a vast agricultural productive area. In addition, one of the largest hydraulic works in the country is located within the region, the Huites dam, one of the principal sources of water supply [1]. In this area, there are frequent extreme rainfall events that cause flooding and serious damage to crops. However, equally extended periods of drought affect this region. These problems make it necessary to pay attention to the determination of regional parameters of rainfall distribution [2,3]. This region has been the location of a considerable number of studies because it is an area subject to extreme maximum and minimum phenomena such as floods and droughts [4,5]. Univariate, mixed and bivariate probability distributions are frequently used to estimate extreme events simultaneously [6,7,8].

However, to successfully characterize the hydrological cycle processes within the region, it is necessary to use regionalization techniques that take into account this particular condition of exhibiting extreme hydrometeorological phenomena [9,10,11]. Likewise, it is necessary to take into account the influence that topographic relief possesses on the hydrological response of the basins, since the Sierra Madre crosses RH10 [12]. As in all countries of the Latin American and Caribbean (LAC) region, regionalization techniques are used to estimate design events at sites with few or scarce data records [13]. These procedures employ clustering algorithms in combination with aggregation rules and metric distances to generate homogeneous groups, from which it is possible to obtain missing hydrological information [14,15,16]. The most recent regionalization was carried out by describing the precipitation regime by means of two parameters, whose numerical values can be associated with a physical meaning. These are the mean number of days with rain during the analysis period (λ) and the mean rainfall height (β) [17]. These parameters, which express a physical meaning, correspond to the numerical value of the parameters of the Poisson and exponential distributions, respectively (leak distribution) [17]. Gutierrez-Lopez et al. (2015) [1] demonstrated that by using these parameters as hydrological characteristics, combined with clustering techniques such as K-means, hydrologically homogeneous regions can be successfully composed [18].

Indeed, the use of probabilistic model parameters that retain physical significance and can also form homogeneous regions is singularly attractive in hydrologic modeling. On the other hand, it is possible to correlate the parameter (λ) with the occurrence of droughts [19] and the parameter (β) with extreme rainfall events [1,17,20]. This advantage of using the leakage distribution makes it possible to characterize the two extreme phenomena, extreme rainfall and droughts, simultaneously with the same probability distribution. However, three mathematical aspects must be previously solved, if maximum and minimum events are to be used, to perform an adequate grouping of hydrological homogeneous regions. These aspects correspond to the objectives of this work and are (i) to verify that the parameters of the mixed Poisson-exponential distribution come from a zero-truncated Gamma distribution and therefore can be operated for hydrological information transfer (regionalization); the later has been mentioned by [21], but the mathematical proof has not been shown; (ii) to obtain spatial clustering of extreme events, it is necessary to propose a metric distance that differs from the Euclidean distance, as a proximity index; this is possible if the distance can be obtained from a correlogram (autocorrelation distance) and a Gaussian variogram (CH-GLo); in this sense, Gutierrez-Lopez (2021) [22] has shown the CH-GLo represents the appropriate mathematical tool to describe the spatial variation in extreme events [22]; (iii) to suggest a modification to the inertia clustering method proposed by Ward in (1963) [23], incorporating the elements described in (i) and (ii) to consider the effect of extreme events simultaneously in the distribution function.

2. The Gamma Distribution from the Exponential Distribution

Definition 1.

The Gamma distribution results from the application of the distribution$x^{v-1}/\mathsf{\Gamma }\left(v\right)$on a simple increasing exponential distribution of the type $e^{-x}$ [24], i.e.,

$$\left[\frac{x^{v-1}}{\Gamma \left(v\right)},e^{-x}\right]\left(v\right)={{\displaystyle \int }}_0^x\frac{t^{v-1}}{\Gamma \left(v\right)}{e}^{-t}dt=\frac{1}{\Gamma \left(v\right)}\gamma \left(v,x\right)$$

(1)

The normality condition is:

$$\left[\frac{x^{v-1}}{\Gamma \left(v\right)},e^{-x}\right]\left(v\right)={{\displaystyle \int }}_0^{+\infty }\frac{x^{v-1}}{\Gamma \left(v\right)}{e}^{-x}dx=1$$

where

$$\rho \left(x,v\right)=\left\{\begin{array}{cc}\frac{x^{v-1}}{\Gamma \left(v\right)}{e}^{-x}& x>0\\ 0& x\le 0\end{array}\right.$$

Introducing a scaling factor $c>0$ we have

$$\rho \left(x,v\right)=\left\{\begin{array}{cc}\frac{c{\left(\frac{x}{c}\right)}^{v-1}}{\Gamma \left(v\right)}{e}^{-\frac{x}{c}}& x>0\\ 0& x\le 0\end{array}\right.$$

Whose characteristic function and Laplace transform are

$$f_{c\left(k,v,c\right)}=1/{\left(1-ick\right)}^v$$

Lemma 1.

The Chi distribution, $y=t^2/2, dy=tdt$

$$\begin{gathered} \left[\frac{x^{v-1}}{\Gamma \left(v\right)},\frac{1}{2^{\frac{v}{2}-1}}\frac{1}{\Gamma \left(\frac{v}{2}\right)}{e}^{-x^2/2}\right]\left(v\right)={{\displaystyle \int }}_0^x\frac{t^{v-1}}{\Gamma \left(v\right)}\left(\frac{1}{2^{\frac{v}{2}-1}}\frac{1}{\Gamma \left(\frac{v}{2}\right)}{e}^{-t^2/2}\right)dt \\ =\frac{1}{2^{\frac{v}{2}-1}}\frac{1}{\Gamma \left(\frac{v}{2}\right)\Gamma \left(v\right)}{{\displaystyle \int }}_0^{\frac{x^2}{2}}{\left[\left(2y^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\right]}^{v-1}{e}^{-y}\frac{dy}{{\left(2t\right)}^{\frac{1}{2}}} \\ =\frac{1}{2^{\frac{v}{2}-1}}\frac{1}{\Gamma \left(\frac{v}{2}\right)\Gamma \left(v\right)}{2}^{\frac{1}{2}v-\frac{1}{2}-\frac{1}{2}}{{\displaystyle \int }}_0^{\frac{x^2}{2}}\frac{{\left(y\right)}^{\frac{1}{2}v-\frac{1}{2}-\frac{1}{2}}}{\Gamma \left(v\right)}{e}^{-y}dy \\ =\frac{1}{\Gamma \left(\frac{v}{2}\right)\Gamma \left(v\right)}{{\displaystyle \int }}_0^{\frac{x^2}{2}}{y}^{\frac{1}{2}v-1}{e}^{-y}dy=\frac{1}{\Gamma \left(\frac{v}{2}\right)\Gamma \left(v\right)}\gamma \left(\frac{v}{2},\frac{x^2}{2}\right) \end{gathered}$$

Using a fractional derivative we obtain

$$\frac{x^{v-1}}{\Gamma \left(v\right)}=D_x^n\left(\frac{x^{v-1+n}}{\Gamma \left(v+n\right)}\right), \forall n\in \mathbb{N}$$

Lemma 2.

As a particular solution it can be written as

$${\left.\frac{x^{v-1}}{\Gamma \left(v\right)}\right|}_{v=0}={\delta }_0\left(x\right)$$

For $f\in {\mathrm{C}}_0^{\infty }\left({ℝ}_{+}\right)$ is defined as

$$\frac{x^{v-1}}{\Gamma \left(v\right)}f={{\displaystyle \int }}_0^x\frac{{\left(x-t\right)}^{v-1}}{\Gamma \left(v\right)}f\left(t\right)dt={\left.D^{-v}\right|}_0 f$$

In particular,

$${\left.D^{-v}\right|}_0 f=f$$

$${\left.D^{-1}\right|}_0 f={{\displaystyle \int }}_0^{x}f\left(t\right)dt$$

$${\left.D^{+1}\right|}_0 f=f^{\prime }$$

In general, ${\left.D^{-v_1}\right|}_0{\left.D^{-v_2}\right|}_0 f=D^{-(v_1+v_{2)}}f$.

When $f\in {\mathrm{C}}_0^{\infty }\left({ℝ}_{+}\right)$ then ${\left.D^{-v}\right|}_{0}f\in {\mathrm{C}}_0^{\infty }\left({ℝ}_{+}\right)$, $\forall v\in ℂ$ and the Fourier transform is $\int e^{-ikx}\frac{x^{v-1}}{\mathsf{\Gamma }\left(v\right)}dx=e^{\frac{1}{2}\pi iv}\frac{1}{{\left(k-i0\right)}^v}$ with $\frac{1}{{\left(k-i0\right)}^v}$ as marginal distribution.

Definition 2.

The characteristic distribution function of the Gaussian is

$$E\left(e^{ikx}\right)=\int e^{ikx}\frac{e^{-\frac{x^2}{2{\sigma }^2}}}{\sqrt{2\pi \sigma }}dx=\int e^{ikx}\frac{1}{\sqrt{2\pi \sigma }}{e}^{ikx-\frac{x^2}{2{\sigma }^2}}dx$$

Reducing

$$E\left(e^{ikx}\right)=e^{-kc+\frac{c^2}{2{\sigma }^2}}\int \frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{x^2}{2{\sigma }^2}-ix\left(\frac{c}{{\sigma }^2}-k\right)}dx$$

By going to $c={\sigma }^{2}k, \int e^{-\frac{x^2}{2{\sigma }^2}}dx=\sqrt{\pi 2{\sigma }^2}$

$$\begin{gathered} e^{-kc+\frac{c^2}{2{\sigma }^2}}\int \frac{1}{\sqrt{2\pi \sigma }}{e}^{-\frac{x^2}{2{\sigma }^2}}dx=e^{-k{\sigma }^{2}k+\frac{{\left({\sigma }^{2}k\right)}^2}{2{\sigma }^2}} \\ =e^{-{\sigma }^2{k}^2+\frac{{\sigma }^2{k}^2}{2}} \\ =e^{-\frac{1}{2}{\sigma }^2{k}^2} \end{gathered}$$

So that

$$f_{c\left(k\right)}=E\left(e^{ikx}\right)=\int e^{ikx}\frac{e^{-\frac{1}{2}{\left(\frac{x}{\sigma }\right)}^2}}{\sqrt{2\pi \sigma }}dx=e^{-\frac{1}{2}{\left(\sigma x\right)}^2}$$

In the same way, it is

$$\begin{gathered} E\left(e^{ik{x}^{\frac{1}{2}}}\right)=\int e^{ik{x}^{\frac{1}{2}}}\frac{1}{2^{\frac{\alpha }{2}}}\frac{x^{\frac{\alpha }{2}-1}}{\Gamma \left(\frac{\alpha }{2}\right)}{e}^{-\frac{x}{2}}dx \\ =\frac{1}{2^{\frac{\alpha }{2}}\Gamma \left(\frac{\alpha }{2}\right)}\int x^{\frac{\alpha }{2}-1}{e}^{ik{x}^{\frac{1}{2}}-\frac{x}{2}}dx \\ =\frac{1}{2^{\frac{\alpha }{2}}\Gamma \left(\frac{\alpha }{2}\right)}\int x^{\frac{\alpha }{2}-1}{e}^{-\frac{x}{2}\left(1-2ik{x}^{-\frac{1}{2}}\right)}dx \end{gathered}$$

Taking $x=x+ic$ then the expression $e^{ik{\left(x+ic\right)}^{\frac{1}{2}}-\frac{{\left(x+ic\right)}^2}{2{\sigma }^2}}$ can be written as

$$\begin{gathered} e^{x\left[ik{\left(x+ic\right)}^{\frac{1}{2}}-\left(\frac{x^2-c^2+i2xc}{2{\sigma }^2}\right)\right]}=e^{\left[ik{\left(x+ic\right)}^{\frac{1}{2}}-\frac{1}{2{\sigma }^2}{x}^2+\frac{1}{2{\sigma }^2}{c}^2-i\frac{1}{2{\sigma }^2}2xc\right]} \\ =e^{\left(+\frac{1}{2{\sigma }^2}{c}^2\right)}{e}^{\left(-\frac{1}{2{\sigma }^2}{x}^2-i\frac{1}{2{\sigma }^2}2xc\right)} \end{gathered}$$

$$E\left(e^{ik{x}^{\frac{1}{2}}}\right)=\int e^{ik{x}^{\frac{1}{2}}}\frac{e^{-\frac{x^2}{2{\sigma }^2}}}{\sqrt{2\pi \sigma }}dx$$

Similarly, when $x=x+ic$

$$\begin{gathered} e^{ik\left(x+ic\right)-\frac{{\left(x+ic\right)}^2}{2{\sigma }^2}}=e^{ikx-kc-\frac{x^2-c^2+i2xc}{2{\sigma }^2}} \\ =e^{\left(-kc+\frac{c^2}{2{\sigma }^2}\right)}{e}^{\left[-\frac{x^2}{2{\sigma }^2}-ix\left(\frac{c}{{\sigma }^2}-k\right)\right]} \end{gathered}$$

So, it is possible to write

$$\begin{gathered} =\frac{1}{\sqrt{2\pi \sigma }}\int e^{ik{x}^{\frac{1}{2}}-\frac{x^2}{2{\sigma }^2}}dx \\ =\frac{1}{\sqrt{2\pi \sigma }}\int e^{-\frac{x^2}{2{\sigma }^2}\left(1-2{\sigma }^{2}ik{x}^{-\frac{3}{2}}\right)}dx \\ =\frac{1}{\sqrt{2\pi \sigma }}\int e^{ik{\left(x+ic\right)}^{\frac{1}{2}}-\frac{{\left(x+ic\right)}^2}{2{\sigma }^2}}dx \end{gathered}$$

At this point, this distribution can be applied to an exponential distribution, as follows:

$$\frac{1}{2^{\frac{\alpha }{2}}}{{\displaystyle \int }}_0^x\frac{t^{\frac{\alpha }{2}-1}}{\Gamma \left(\frac{\alpha }{2}\right)}{e}^{-\frac{t}{2}} dt=\frac{1}{2^{\frac{\alpha }{2}}}\left[\frac{x^{\frac{\alpha }{2}-1}}{\Gamma \left(\frac{\alpha }{2}\right)},e^{-\frac{t}{2}}\right] \left(\frac{\alpha }{2}\right)$$

A similar result can be found in [25]. Now, employing the reduced variable $y=\frac{t}{2}$, it is obtained $dy=\frac{1}{2}dt$, and if it is used, the form of a function of the type $f\left(b,x\right)={{\displaystyle \int }}_0^x{t}^{b-1}{e}^{-t}dt$ then

$$\frac{1}{2^{\frac{\alpha }{2}}}{{\displaystyle \int }}_0^x\frac{t^{\frac{\alpha }{2}-1}}{\Gamma \left(\frac{\alpha }{2}\right)}{e}^{-\frac{t}{2}} dt=\frac{1}{2^{\frac{\alpha }{2}}}{{\displaystyle \int }}_0^{x/2}\frac{2^{-\left(\frac{\alpha }{2}-1\right)}{y}^{\frac{\alpha }{2}-1}}{\Gamma \left(\frac{\alpha }{2}\right)}{e}^{-y}2dy$$

(2)

$$=\frac{1}{\Gamma \left(\frac{\alpha }{2}\right)}f\left(\frac{\alpha }{2},\frac{x}{2}\right)$$

(3)

Thus, from the above process, the result is a truncated Gamma function.

3. Autocorrelation Length

Definition 3.

A normalized semi-variogram is defined as $\gamma {\left(h\right)}^{*}=\frac{\gamma \left(h\right)-{\gamma }_0\left(h\right)}{{\gamma }_{\infty }\left(h\right)-{\gamma }_0\left(h\right)}$ (Figure 1a). Because the variogram and the correlogram are two estimators of the variability in hydrological data [26,27], a standardized correlogram (Figure 1b) can also be defined as:

$$\rho {\left(h\right)}^{*}=1-\gamma {\left(h\right)}^{*}$$

(4)

Lemma 3.

The length of the autocorrelogram $L_c$can be approximated by an expression of the form $L_c={{\displaystyle \int }}_0^{\infty }\rho {\left(h\right)}^{*}dh$.

Lemma 4.

For the case of an exponential distribution according to Equation (1) we have:

$$\rho {\left(h\right)}^{*}=e^{-\left(\frac{h}{Le}\right)}$$

$$L_c={\left.{{\displaystyle \int }}_0^{\infty }{e}^{-\left(\frac{h}{Le}\right)}dh=-L_e{e}^{-\left(\frac{h}{Le}\right)}\right]}_0^{\infty }=-L_e\left[0,-1\right]\to L_c=L_e$$

The exponential variogram (exponential model) has been used for several years as the model that best characterizes the variability in extreme precipitation in mountainous regions [28].

Lemma 5.

For the case of a Gaussian distribution according to Equation (1) we have:

$$\rho {\left(h\right)}^{*}=e^{-{\left(\frac{h}{Lg}\right)}^2}$$

(5)

$$L_c={{\displaystyle \int }}_0^{\infty }{e}^{-{\left(\frac{h}{Lg}\right)}^2}dh=L_c{{\displaystyle \int }}_0^{\infty }{e}^{-{\left(\frac{h}{Lg}\right)}^2}dh=L_c\frac{\sqrt{\pi }}{2}$$

$$L_c=\frac{\sqrt{\pi }}{2}{L}_g\therefore L_g=\frac{2L_c}{\sqrt{\pi }}$$

(6)

which we are defining as autocorrelation ratio.

4. Autocorrelation Ratio as a Measure of Inertia

Expression (1) is a particular case of the Gamma–Laguerre function presented by [29]. The general expression has been used with great success as a method of analysis for extreme hydrological values [29]. Specifically, it has been proven that when the Gamma function multiplied by a series of orthogonal Laguerre-type polynomials is presented [30,31], the probabilistic modeling is appropriate and relevant since it uses higher order moments in the construction of the mathematical model. In particular, it is recommended that when the absolute value of the asymmetry coefficient of a series of maximum values is large, a Gamma–Laguerre GL(x) type function should be used and generally the number of moments to be used is between 3 and 4.

$$GL\left(x\right)=\frac{e^{-\frac{t}{2}} t^{\frac{\alpha }{2}-1}}{{\left(2\beta \right)}^{\frac{\alpha }{2}} \Gamma \left(\frac{\alpha }{2}\right)}$$

(7)

In expression (5), the parameters $\beta$ and $\alpha$ represent the scale and the number of degrees of freedom, respectively. From [30] and using expression (2) and simplifying it as a particular case of (3), we have:

$$\frac{e^{-\frac{t}{2}} t^{\frac{\alpha }{2}-1}}{{\left(2\beta \right)}^{\frac{\alpha }{2}} \Gamma \left(\frac{\alpha }{2}\right)}\approx \frac{\sqrt{\pi }{L}_g}{2 \Gamma \left(\frac{\alpha }{2}\right)}$$

Lemma 6.

This means that the Gamma distribution can effectively be used as a weight function to perform a distribution weighting, when talking about extremes. The autocorrelation ratio can also be introduced as a scaling parameter.

5. The Inertia Criterion for the Classification of Extremes

The classification of extreme events is fundamental to hydrological studies and represents one of the stages of regionalization [32,33]. Usually a clustering method is linked to calculations of inertia metrics in each of the stages of aggregation of elements. However, it has been shown that a single station with extreme value measurements has influence on the behavior of the whole region [22,34,35]. If we consider that $m_i$ y $m_j$ are the elements to be grouped in the dimensional space $\left(i,j\right)$, their momentum M is a function of the weights of the elements $m_{i,j}$ of the inertia they cause ${\theta }_{i,j}$ and of a representative distance of each element $L_{i,j}$. This can be written as $M=f\left\{m_i, m_j,{\theta }_i, {\theta }_j,L\left(i,j\right)\right\}$. The Minkowsky metric index is usually used, which is a multidimensional distance [36].

$$L\left(i,k\right)={\left[{\displaystyle \sum }_{j=1}^d{\left|m_{i,j}-m_{k,j}\right|}^r\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$r$}\right.}r\ge 1$$

When $r=2$ it is known as the Euclidean distance [37]. If it is used as a proximity estimator to the autocorrelation ratio $L_{i,j}=f\left(L_c\right)$, then it can be written as

$$M\left\{{\theta }_i, {\theta }_j\right\}=m_i{L}^r\left({\theta }_i, {\theta }_j\right)+m_j{L}^r\left({\theta }_i, {\theta }_j\right)$$

Now, if we use the autocorrelation ratio as a scalar number, we obtain the following:

$$L^r\left(\theta , {\theta }_i\right)\le \frac{m_i{\theta }_i+m_j{\theta }_j}{m_i+m_j}-{\theta }_i$$

$$\frac{m_i{\theta }_i+m_j{\theta }_j}{m_i+m_j}-{\theta }_i\ge \frac{m_j^r}{{\left(m_i+m_j\right)}^r}<{\theta }_j-{\theta }_i$$

$$L^r\left(\theta , {\theta }_i\right)=\frac{m_j^r}{{\left(m_i+m_j\right)}^r}{L}^r\left({\theta }_i, {\theta }_j\right)$$

$$M\left\{{\theta }_i, {\theta }_j\right\}=\frac{m_i{m}_j^r}{{\left(m_i+m_j\right)}^r}{L}^r\left({\theta }_i, {\theta }_j\right)+\frac{m_j{m}_i^r}{{\left(m_i+m_j\right)}^r}{L}^r\left({\theta }_i, {\theta }_j\right)$$

Finally, taking the autocorrelation ratio as a function, we have:

$$M\left\{{\theta }_i, {\theta }_j\right\}=\frac{m_i{m}_j}{m_i+m_j}{L}_c^r\left({\theta }_i, {\theta }_j\right)$$

$$M\left\{{\theta }_i, {\theta }_j\right\}=\frac{m_i{m}_j}{m_i+m_j}{L}_{i,j}\left({\theta }_i, {\theta }_j\right){\left(\frac{\sqrt{\pi }}{2}{L}_g\right)}^r$$

(8)

This expression is similar to the one proposed by Ward in (1963) [23] and is still in use today for the clustering of similar values. This variant allows a Gaussian behavior in the spatial distribution of extreme values to be considered.

6. Methodological Application and Results

RH10 has an area of about 80,000 km², is located in northwestern Mexico and is divided into nine watersheds. There are 93 meteorological stations with historical records (50 years on mean) of maximum daily rainfall. The methodology proposed in this work to delimit homogeneous regions, in areas subject to maximum and minimum extreme events, is detailed as follows:

• This realization began with the sampling at the 93 of rain gauge stations, giving the measurements. The probabilistic Poisson-Exponential model is fitted to the daily rainfall data sample. Next the parameters of the model are calculated: $\lambda$ as number of days with rain and mean rainfall per event as $\beta$ [1,17]. The proposal is that since these parameter values provide a physical meaning, they can then be used as weights for the aggregation algorithm, i.e., $\left(m_i, m_j\right)\approx \left({\lambda }_i,{\beta }_j\right)$;
• The spatial structure of each rainy event is characterized by its variogram $\gamma \left(h\right)$ [38]. Let $(s_i,s_j)$ denote the observed value at coordinates $s=\left(x,y\right)$; then the Robust Gaussian Variogram Estimator is estimated, where N(h) is the number pairs of points $\left(x_i, y_j\right)$ a distance away $h=\left|x_i-y_i\right|$. The scaled variogram is defined by [22] as

  $$\gamma \left(h\right)=3{\left(\frac{1}{2}\sum \sqrt{\left|s_i-s_j\right|}\right)}^2/{\left(0.457+\frac{0.494}{\left|N\left(h\right)\right|}+\frac{0.045}{{\left|N\left(h\right)\right|}^2}\right)}^{1/2}$$

  (9)
• From $\gamma \left(h\right),$ the value of the autocorrelation length $L_g=\frac{2L_c}{\sqrt{\pi }}$ is obtained using Equations (4)–(6). With this length, it is now possible to affect the values of the Euclidean distance matrix $L\left(i,k\right)\to r=2$. In this way we obtain the term $L_{i,j}\left({\theta }_i, {\theta }_j\right){\left(\frac{\sqrt{\pi }}{2}{L}_g\right)}^r$.
• Finally, Equation (8) is applied iteratively, until the hydrological homogeneous regions are obtained.

The aggregation proposed in this paper is based on using the value of the parameters for a mixed probability distribution as weight factors in the aggregation algorithm. Likewise, the autocorrelation length is used to take into account the effect of the spatial variation of extreme events within the distance matrix. The result of this procedure applied to RH10 is presented in Figure 2. The symbols in triangles in Figure 2 and Figure 3 represent the climatological stations; the colors relate to each of the groups composed.

7. Discussion

Delimitation of homogeneous regions has been typically based on the identification of geographically contiguous zones [40,41]. It has been supposed that the grouping of homogeneous zones is only related to the hydrological response of watersheds [42,43]. However, the current climate and spatial variability of exceptional events require mathematical methods to characterize such extreme variability [6,34]. A hydrological regionalization operating parameters with physical meaning as weighting factors seems to obtain advantages and allows prior information to be incorporated into the analysis [44]. In addition, a mixed probability distribution allows characterizing maximum and minimum extreme events, once the primary hypothesis proposed by [21] and verified in Section 2 of this paper is supported.

To verify the likelihood of the proposed methodology, it is vital to compare the delimitation obtained by other authors. Campos-Aranda (1994) [4] employed the Lagbein test carry out to hydrometric records to regionalize RH10, and did not find subgroups. At that time, in 2014 [39], they employed a principal component analysis (PCA) and hierarchical ascending method procedure. They recognized three regions at different altitudes, all parallel to the coast and parallel to the Sierra Madre (Figure 3). This result seems to be totally influenced by orography [9]. Studies have shown that the delimitation of homogeneous regions is in fact affected by physiography [45,46] and especially when a method such as hierarchical ascending is used [47]. Years later, a mapping of RH10 was presented from a Kriging of the parameters $\left({\lambda }_i,{\beta }_j\right)$; the spatial interpretation of this procedure allows for the identification of two homogeneous regions: a high-altitude area on the Sierra Madre and a low-altitude region on the coastal flat-plain [1]. The results of the delimitation for these studies are presented in detail in Table 1 and graphically in Figure 3.

The geographic subdivision of this critical region in Mexico has been investigated for years. This area comprises a portion of the Sierra Madre mountain chain, which supplies water to more than a quarter of Mexico’s total land area. Any change in the hydrological regime, be it a maximum or a minimum, can produce severe consequences for the development of this region. Floods are a consequence of extreme events such as hurricanes that affect RH10 every year. In the case of droughts, RH10 is peculiarly sensitive if we consider that 54% of the fields suitable for cultivation are located in northern Mexico. Unfortunately, this region only receives 7% of the total available water for the whole country [48].

Table 1. Techniques used for the delimitation of homogeneous regions.

Author	Procedure	Final Regionalization
Campos-Aranda 1994 [4]	Langbein test [49,50]	The whole region is homogeneous.
Arellano-Lara, and Escalante-Sandoval (2014) [39]	Principal component analysis [51] and hierarchical ascending clustering [52]	Three regions parallel to the coast and parallel to the Sierra Madre (Figure 3).
Gutierrez et al. (2015) [1]	Cartography of $\left({\lambda }_i,{\beta }_j\right)$ [53]	Two regions: one on the Sierra Madre and the other on the coastal flat-plain (Figure 3).
Current proposal (2022)	Equation (8)	Three regions: coastal flat-plain, mountainous region and a coastal desert (Figure 2).

Table 1. Techniques used for the delimitation of homogeneous regions.

Author	Procedure	Final Regionalization
Campos-Aranda 1994 [4]	Langbein test [49,50]	The whole region is homogeneous.
Arellano-Lara, and Escalante-Sandoval (2014) [39]	Principal component analysis [51] and hierarchical ascending clustering [52]	Three regions parallel to the coast and parallel to the Sierra Madre (Figure 3).
Gutierrez et al. (2015) [1]	Cartography of $\left({\lambda }_i,{\beta }_j\right)$ [53]	Two regions: one on the Sierra Madre and the other on the coastal flat-plain (Figure 3).
Current proposal (2022)	Equation (8)	Three regions: coastal flat-plain, mountainous region and a coastal desert (Figure 2).

8. Validation

How can the most appropriate regionalization be identified? The work presented by Drescroix et al. (2007) [54] to attempt to define the spatial and temporal boundaries of Hortonian (infiltration excess runoff) and Hewlettian (saturation excess overland flow) hydrological behavior found the most functional regionalization for northern Mexico. This study is based on data collected in four diverse areas of northern Mexico. The experimental sites were equipped with stream gauges and rain gauges. This study determined that the Sierra Madre is mostly characterized by a humid zone with Hewlettian hydrology. Another zone with Hortonian runoff completely dominates the semi-arid and arid areas. However, in certain areas saturation excess overland flow can appear due to the genuine possibility of the landscape infiltrating a significant proportion of rainwater, i.e., during hurricanes crossing the mountains. It has been shown that the relationship between a soil’s classification and the spatial parameters of a conceptual catchment scale is very important in the hydrological models [55]. This region would represent a coastal and mountainous area, identified as the transition zone to the Chihuahua Desert. This zone can be evidenced clearly in the regionalization work proposed therein (Figure 2). In detail, it can be remarked that each region maintains extremely specific characteristics, described below and in Appendix A.

Region I: This region is a transition zone to the Chihuahua desert. It is a region with a mean altitude of 400 m. As a transition zone, it maintains altitudes from 19 m to 2300 m. Near the northern deserts, the annual precipitation is only 500 mm. The mean value of $\lambda$ is 0.17 with a mean of 11 mm per event $\beta$. The autocorrelation ratio is between 1.6 and 126.92.

Region II: This region is a semi-arid zone, and a coastal desert with Hortonian runoff. It is a region with mean altitudes of 129 m above sea level. The coastal desert zone gains altitudes from 12 m to 800 m. The region includes an annual mean of 549 mm of rainfall per year. The average value of $\lambda$ is 0.16 with an expected 12.6 mm per event $\beta$. The autocorrelation ratio is between 1.1 and 21.5.

Region III: This region is a humid zone with Hewlettian hydrology at Sierra Madre. The region includes an altitude range of 1299 m. It is the high mountain zone with altitudes from 140 m to 2720 m. It is the rainiest region with an expected rainfall of 718 mm during the rainy season. The significant value of $\lambda$ is 0.25 with an expected 10.7 mm per event $\beta$. The autocorrelation ratio is between 1.8 and 7.2.

9. Conclusions

It was properly verified that the parameters of the mixed Poisson-exponential distribution come from a Gamma distribution truncated at zero. Therefore, they can be utilized for the hydrological information transfer of extremes. In this case, they are employed as weight factors in the aggregation algorithm to obtain hydrological homogeneous regions. Moreover, these parameters retain a physical meaning. Days of rainfall and mean rain per event create certainty in the weights operated in the aggregation algorithm. As a complement to the Euclidean distance, the autocorrelation ratio was added from a Gaussian variogram. The convergence in the clustering of regions is more efficient and faster when using Equation (6). The final clustering was compared with previous studies, obtaining similar and even more revealing results. The methodology was validated with the results from research that conducted experimental studies on runoff and soil properties in that region of northern Mexico. Therefore, it is considered a valid, acceptable and efficient proposal to delimit homogeneous regions subject to extreme events.

It is important to mention that the verification of hydrological homogeneous regions should be carried out with cross-validation [33]. It was not the objective of the current paper to perform that validation. Furthermore, this is a hydrological topic. The formulation of a mathematical expression to reach an aggregation of elements subject to extreme conditions is presented in this paper.