Mixed statistical and data mining analysis of river flow and catchment properties at regional scale

This paper presents an application of the L-moments and L-moment ratio diagrams (LMRD) to the analysis of hydrological data at regional (country) scale. Existing research focuses on two main areas of the analysis: statistical analysis using LMRD and regression analysis. Further research mixes both approaches applying regression analysis to L-moments. Another direction of the research is clustering of the climatic and physiographic catchment properties and its validation using LMRD. However, LMRD plots can be separately used as the clustering domain. It is proposed to decompose the features into some classes, and than present these results on the LMRD. Such plots constitute the source for the clustering. Obtained clusters are then validated against k-means clustering performed in the LMRD diagram domain. Results show that statistical L-moments analysis can be improved with data mining clustering algorithms. Such combination delivers a new perspective for the interpretation of the results. It is shown that clustering in the LMRD domain is consistent with the K-means clustering. It is anther argument showing that L-moments diagrams can be considered as a very powerful and informative tool for hydrologists enabling the comparison on the regional basis with respect to various catchment properties. The method is validated on data consisting of daily river flow data from 290 gauges covering entire Poland.


Introduction
Investigating the effects of catchment properties related to climate and physiography on river flow characteristics is a frequent area of research in hydrology. Such an analysis is mostly based on statistical approaches (Helliwell et al. 2007;Calver et al. 2009;Salinas et al. 2014;Singh 2017). Actually, the analysis can be conducted on the single gauge basis (Popat et al. 2020;Negi et al. 2021) or on the higher regional level (Smith et al. 2015;Kar et al. 2017). Regional analysis uses different approaches. Generally, there are two leading methodologies. Application of the method of L-moments in the regional frequency analysis of extreme events for river flows, apart from other natural sciences, is considered to be one possible approach (Hosking and Wallis 1993;Kjeldsen et al. 2017;Simková 2017). The second methodology uses regression analysis Mašiček and Toman, Palàt 2011;Lun et al. 2021). There is also research mixing both approaches (Di Baldassarre et al. 2006;Salinas et al. 2014), i.e. application of the regression analysis to correlate the L-moment diagram data.
Both approaches have their advantages and shortcuts. L-moments ratio diagrams (LMRD) origin from the traditional moment ratio diagrams (MRD) introduced by Karl Pearson in the beginning of the XIX century. The idea is to use and graphically represent the relationship between the first four moments of the distribution. There are various versions of MRD, but the diagram showing the relationship between the fourth and the third moment, i.e. kurtosis versus skewness is the most popular one. MRD of the variation factor versus skewness is the second popular representation.
MRDs can show empirical data or analytical theoretical curves. Theoretical relations can be evaluated for the majority of univariate distributions (Vargo et al. 2010). MRD can be used for several purposes, such as: 1. to quantify the distance between various univariate distributions, 2. to visualize distribution versatility within the range of its moments, 3. to select the best fit distribution to the empirical data, 4. to visualize and test the relationships between distribution families, 5. and to identify homogeneous process expressed by the empirical moments relations.
Historical formulation uses classical moments definitions as for normal Gaussian distribution. Introduction of the Lmoments by Hosking (1990) enabled to increase the approach applicability. L-moments offer better moments estimation properties, especially for small sample number and fat-tailed distributions (outliers). It has been shown that they might even outperform other estimation approaches, like the standard algorithm of moments, TL-or LQ-moments (Simková 2017). L-moments and respective L-moments diagrams are widely used in the extreme events analysis in climatology, hydrology (Katz et al. 2002;Maeda et al. 2013), astronomy (Podladchikova et al. 2003) or medicine (Louzada et al. 2016). LMRD is the most common tool supporting the identification of a suitable frequency distribution of empirical samples (Kjeldsen and Prosdocimi 2015). Typical diagram represents the Lkurtosis (s 4 ) versus L-skewness (s 3 ). Empirical data may be easily confronted with the theoretical plots for the selected candidate probabilistic density functions (PDFs). The proximity of the empirical data to the theoretical ones might play the role of a selection criterion to choose the best fit distribution (Peel et al. 2001). Apart from distribution fitting, L-moments diagrams are often utilized to compare various samples originating from different sources, like for instance river discharge, in search for their homogeneity. There are defined various discordance measures used to recognize sources, whose sample's L-moments are marked contrarily from the others (Hosking and Wallis 1993;Khan et al. 2017). The task associated with the homogeneity testing is realized in the 2-dimensional (2D) space spanned by L-skewness and L-kurtosis. Verification of the homogeneity hypothesis may be performed by the visual inspection or may be supported by the dedicated distance measures. L-Moment Ratio Diagrams deliver new perspective to the statistical analysis of various time series properties. These properties are described by the aggregated position of the time series on the diagram. Thus, they allow to perform regional analysis comparing different features associated with the data. This feature is the strong advantage of the approach. However, the method does not give any insight into the specific statistical properties. Regional analysis is limited to the general clustering of the observations without connection to the internal (mean, L-Cv, L-skewness and L-kurtosis) moments.
Regression analysis (Garmdareh et al. 2018;Desai and Ouarda 2021) is applied to address this issue. It aims at finding a quantitative relationship between certain catchment properties and respective flow time series L-moments. This approach is quite natural and promising. However, one has to take into account that standard regression using least squares (LS) has its limitations. It assumes that the process behind data is stationary, affected by Gaussian noise and without outliers (Hawkins 1980;Domański 2020). If these conditions are not met, the LS regression estimation will be biased (Rousseeuw and Leroy 1987;Huber and Ronchetti 2009). Moreover, LS regression assumes that the observations are relatively uniformly distributed over the domain number of observations and their number is high enough to satisfy convergence towards the limiting probabilistic density function (PDF). Unfortunately, hydrological data does not meet these assumptions. The process underlying statistical process are neither stationary, nor Gaussian. Outliers are frequent, what is visible in the tails fatness and the resulting need for the extreme analysis. The length of the time series is not long enough to assure the convergence. Therefore, obtained estimations can be susceptible to the bias.
Researchers use different approaches to minimize the effect, as for instance weighted moving average (WMA) (Salinas et al. 2014). The present study does not focus on that aspect, though sample regression analysis in different context is included. It is proposed to use robust regression performance index in form of the least median square (LMS) (Rousseeuw 1984).
Following above arguments considering LMRDs and regression analysis, the first approach (LMRD clustering) is selected as it opens a new research opportunity. This opportunity arises from the fact that LMRD plots enable various approaches to the data comparison in the 2D space of the moment ratio diagram. Such a statistical procedure resembles the 2D grouping tasks well known in the data mining research. One may find dozens of clustering algorithms that aim at data grouping, finding cluster centers or identifying the delimitation borders (Aggarwal and Reddy 2014;King 2015). Thus, the use of clustering methods in the grouping task in the L-moments diagrams data seems to be natural and obvious. As the grouping itself is not an ultimate focus of this research, classical k-means algorithm is used (Lloyd 1982), as it is one of the simplest and popular unsupervised machine learning methods.
Recent literature review shows some reports on clustering approaches used in the hydrological analysis. The Ward's data agglomeration method is a hierarchical approach and uses analysis of variance as the measure of the distance between clusters It has been used for catchment hydrological analysis (Cupak 2017;Cupak et al. 2017). A principal component analysis and a subsequent clustering of the principal components has been used to cluster catchment data (Jehn et al. 2020), while authors of Sharghi et al. (2018) compared four approaches of clustering, K-means, Ward, SOM and Wavelet-Entropy-based method in a similar task. K-means clustering has been applied in the analogous research in Dikbas et al. (2013), Aytaç (2020).
Step-wise cluster analysis hydrological model has been proposed in Wang et al. (2021). Actually, current research is frequently organized as follows: researchers select a clustering algorithm, do catchments' clustering by their features and finally observe the results in the LMRD. Proposed approach is opposite. First, LMRD are plotted and the K-means clustering is performed in the diagram 2D domain. Clustering results are compared with expert and statistical features classification.
Concluding, this study aims at the analysis of climatic and physiographic catchment properties in the domain of the L-moments ratio diagram based on the river daily flow data from 290 flow gauges from Polish rivers. Catchment properties are compared with the L-moments of the river flow data and the best fitted PDFs. The analysis is supported by the k-means clustering to identify leading properties and to allow the generalization of observations. The study starts with Sect. 2 presenting the considered geographical area and respective data. Regional clustering and analysis of properties is included in Sections 4, preceded by Sect. 3 describing applied methods and algorithms and followed by Sect. 5 with concluding remarks and open research items.
2 Study area and data Regional hydrological analysis is performed using daily river flow data ½m 3 =s provided by the Institute of Meteorology and Water Management -National Research Institute (IMGW-PIB) from 290 gauges located in Poland. The analysis uses available historical data (daily discharges) starting from November 1950 till the end of October 2019. The map showing the area of interest is presented in Fig. 6. Moreover, the analysis uses a set of catchment climatic and physiographic properties, i.e. the catchment area [km 2 ], mean annual precipitation [mm], mean annual of the minimum and maximum daily temperatures [deg. C], catchment slope [-], elevation [m asl], dominant land use classes (fractions), dominant soil hydrological groups (fractions), soil permeability classes (fractions) and fraction of lakes [-]. All the properties were calculated in ArcGIS software using readily available climatic and physiographic data for the respective catchment areas. Climatic data were derived from the G2DC-PLþ data set , and physiographic data were derived from the input maps used in the recently developed hydrological model SWAT covering the Polish territory (Marcinkowski et al. 2021).
At first, catchment properties are briefly analyzed, i.e. any possible effect of these features on the catchment flow is checked. Such an analysis is not straightforward, as there are at least are two challenges. First, the appropriate aggregate measure of the river flow variable has to be chosen. Three candidates are taken into account: mean value as the classical statistical shift estimator, median representing robust shift estimator and the maximum value. In order to account for the direct effect of catchment size on river flows, selected flow indices (mean, median and max) were normalized, i.e. divided by the respective catchment area. They are referred to as ''specific'' flows, e.g. a mean specific flow, etc., and their unit is Second issue is associated with the regression analysis itself. First of all, the number of observations is very low, which may bias any obtained estimation. Moreover, the underlying statistical process behind data cannot be assumed neither Gaussian, nor stationary. Thus, apart from classical mean shift estimator, the robust version is used (Rousseeuw 1984). Classical regression uses least squares (1) as the performance index, while the robust one -the LMS (2). The example of this effect is shown in sample Fig. 1.
The estimators differ. LS index takes into account even outlying observations, while the robust one concentrates of data inliers (majority). The decision, which one should be taken into consideration as the resulting index is not straightforward.
As expected, the effect of catchment area on selected specific flow indices evaluated by robust regression is close to 0 due to standardization process (Table 1). Two categories of catchment properties stand out from others: topography and climate. Catchment slope has the highest correlation with mean specific slope across all catchment properties, followed by mean annual precipitation and average elevation. The higher the slope, elevation and precipitation, the higher the mean, median and maximum specific flows. In contrast, average (min and max) temperature have negative correlation with mean and median specific flow, which can be explained by two facts: that higher elevation is associated not only with higher precipitation but also with lower temperature, and that lower temperature results in lower evapotranspiration. Other categories of catchment properties such as land cover, soils and hydrography exhibit less pronounced relationship with specific flow, particularly when assessed with robust regression. Linear regression shows that fraction of arable land has negative correlation with all specific flows, whereas fraction of forests, urban land and grassland has positive correlation. Catchments having a high fraction of impermeable soils have positive correlation with specific flow due to higher surface runoff. In general, the robust regression has significantly different results for maximum specific flows, perhaps because the outliers are most abundant there.
The values of each property are grouped into clusters. The number of clusters is set to 4 (as the results of the selection of quantiles: Q1, Q2 and Q3.). The division is realized in two ways: using hydrological expert knowledge and using statistical quantiles (Q1 ¼ 25%, Q2 ¼ 50% and Q3 ¼ 75%) for respective feature. Data are presented in Table 2. These properties are compared with the best fitted PDFs and L-moments properties.

Methods and algorithms
All the applied methods are specified and introduced in this section. Statistical approach covers used distributions, Q-Q plots, L-moments and moment ratio diagrams. Data mining approach is represented by the k-means clustering algorithm.

Distributions used for the extreme flow analysis
Analysis of the regional river flow data starts with the review of properties for each flow gauge. Hydrographs are plotted together with the histograms. Considered discharge time series are characterized by the periodicity and frequent extreme high value incidents. Histograms for such time series are highly asymmetric. One may find various probabilistic distribution functions appropriate in such a case. The idea of this research was to minimize unnecessary assumptions. Therefore the analysis takes into account seven different distributions being frequently used in the considered research context: Gamma (GAM), lognormal (LGN), exponential (EXP), Weibull (WEI), Generalized Extreme Value (GEV), generalized Pareto (GPD) and four-parameter Kappa (K4D). Respective probabilistic density functions and their parameters are shown in Table 3.
Each PDF should be fitted to the empirical data. The fitting is achieved using maximum likelihood estimation except the four parameter Kappa distribution, for which the L-moments approach is applied. Sample plot showing fitting, i.e. data histogram and fitted PDFs is sketched in Fig. 2.

Q-Q plots and the best fit PDF estimation
In general quantile-quantile (so called Q-Q) plots enable to compare different distributions by plotting their quantiles against each other. A point on the diagram corresponds to a certain quantile coming from an empirical and theoretical distribution. They may be used as the heuristics for testing ''goodness-of-fit'' between some empirical data and a theoretical PDF. It allows to check the assumption the empirical data are derived by the chosen distribution. The fitting is relatively simple. We need to find such a theoretical curve that appears to be the closest to the empirical points. The distance between both curves determines the ''goodness-of-fit'' index.
Formal PDF fitting is done using the maximum likelihood (MLE) approach. It gives reliable results and is simple. It is frequently used used for large datasets (as in our case) (Coles and Dixon 1999). The Q-Q plot is used to select the best fit function. The fitting measure is robust median absolute distance between empirical and theoretical data. It plays the role of the ''goodness of fit'' index supporting a decision about the best fitting distribution. Figure 3 shows a sample Q-Q plot for sample empirical data versus seven theoretical distributions. Permeability: low Relation is presented in form of the linear regression slope angle -higher angle absolute value means stronger impact

L-moments
L-moments have been introduced as linear combinations of order statistics (Hosking 1990). This approach significantly improves conventional methodology and definitions. It gives new characterization of the shape of a probability density function and allows to estimate the distribution parameters. Unlike product moments, L-moments deliver almost unbiased L-moments statistics, even for very small samples. Additionally they are less sensitive to the distribution's tails (Peel et al. 2001). Above properties suite them almost perfectly and allow to describe environmental data commonly characterized by a moderate or high skewness. Evaluation of L-moments is realized through the following procedure. At first, the data xð1Þ; . . .; xðNÞ f g , Nnumber of samples, ranked in ascending order from 1 to N. The sample L-moments (l 1 ; . . .; l 4 ), the sample L-skewness s 3 and L-kurtosis s 4 are evaluated as where Therefore, the data are summarized and described by the sample L-location (L-shift) l 1 , L-scale l 2 , L-Cv (dimensionless measure of variability) s 2 , L-skewness s 3 2 ðÀ1; 1Þ and L-kurtosis s 4 2 ðÀ 1 4 ; 1Þ. L-moments might be used to fit a distribution to a dataset. It is done through fitting of the sample empirical L-moments to the exact theoretical L-moments of the distribution. L-moments s 3 and s 4 may play the role of the goodness-of-fit measure. L-moments can be analytically calculated for different distributions. Theoretical relationships for univariate PDFs can be found in Hosking (1990), Hosking (1992), Kjeldsen et al. (2017).

Moment ratio diagrams
MRDs have been used in practical applications to answer questions about sample's statistical properties, such as theoretical distribution fitting, comparison of the distribution's shapes or PDFs classification according to certain categories (Bobee et al. 1993). The MRD is a graphical representation in a Cartesian coordinates of a pair of standardized moments. In general, there are two types of moment ratio diagrams: a graph of skewness versus kurtosis and a graph of a scale factor versus skewness for common univariate probability distributions (Vargo et al. 2010). The literature shows that the first one has gained significant popularity and is adopted in current research. Introduction of L-moments has been naturally used by the MRDs and the L-moment diagrams have been extensively analyzed. They form the graphical representation of a relationship between L-kurtosis and L-skewness. Figure 4 shows theoretical relationship for selected distributions. One may see that exponential distribution is represented as a single point, four parameters Kappa as a region, while the rest of them by polynomial curves. As  Left-bottom part of the L-moment diagram is occupied by datasets with low L-skewness s 3 and L-kurtosis s 4 . L-skewness closer to zero means that the PDF is more symmetrical. Simultaneously, lower L-kurtosis reflects lighter tails of the respective distribution. In contrary, topright diagram area depicts data, which are highly skewed with significantly fat tails. Therefore, positioning on the L-moment diagram allows fast data interpretation and classification.

Clustering
Data clustering algorithms belong to the machine learning area of science and are used to group data that lie close to each other and may be classified as more similar to each other in some sense, in contrast to other ones. The literature shows dozens of clustering articles. The one used herein is a k-means algorithm. It is one of the simplest unsupervised approaches as it makes inferences from data using only 1. Specify the number of clusters M. 2. Initialize centroids randomly selecting M points for the centroids without replacement. 3. Iterate until the assignment of points to clusters does not change. 4. Compute the performance index being the sum of the squared distance between points and centroids. 5. Assign each point to the closest data cluster. 6. Evaluate new centroids taking the average of all data belonging to each cluster.

Regional clustering results
The results of the regional clustering are presented in the below paragraphs, starting from the PDF goodness-of-fit analysis, followed by the clustering in the domain of L-moments diagrams.

PDF goodness-of-fit analysis
The analysis is grouped into formal steps. At first the goodness-of-fit assessment is done to assign the most appropriate distributions to each flow gauge. In contrast to the majority of reports, the Q-Q diagram approach is used. It allows better rationalization of the choice between fourparameter Kappa distribution and other PDFs, of which some belong to the same family (GEV, GPD, EXP). The results are presented in a graphical form. Figure 5 shows the bar plot with the PDF fitting results. Two versions are presented. In the first one all the distributions are used. As it is seen the overwhelming majority of gauges is well described by four-parameter Kappa (81%) and GEV (18%). Once K4D is excluded from the analysis, GEV distribution covers the majority of sites, however LGN, GAM and GPD PDFs still match quite a significant number of points. Nevertheless, the concluding remark is clear: four-parameter Kappa distribution should be considered as the best fitting probability density function in case of Polish daily river flow data.

L-moment diagram analysis
At first, L-moments are calculated for each flow time series. Next, they are plotted on the L-diagram (see Fig. 7). Additionally, empirical mean and median L-moments are evaluated and respective points added to the diagram. The plot confirms the fact that K4D distributions covers the

k-means clustering
Empirical points of the L-moments diagram can be treated as any two-dimensional data on a plane and as such can be subject to clustering. Figure 8 shows how the k-means algorithm allocates data to the four regions. The number of regions is equal to the number of classes defined for each catchment property as in Table 2. Expert classification is done into four classes. Similarly straightforward approach using basic quantiles Q1, Q2 and Q3 also gives classification into 4 groups.
The clustering plot depicts evaluated regions' clusters. As one can notice, the median of empirical points lies very close to the border between two middle clusters, which is as expected. For clarity, these clusters are called MLclusters (machine learning clusters). Machine learning clustering is further compared with expert based classification depending on catchment area categorization.

Effects of catchment properties
At that point all important elements are prepared and catchment properties analysis can be performed. Each property listed in Table 2 will be analyzed separately according to the same proposed methodology.
(1) Two L-moment plots are prepared for both types of feature cluster definitions, i.e. expert-and quantilebased. The points are classified according to the respective cluster membership. Expert clustering is called EX-clusters, while quantile-based QUclusters.
(2) Median point for each cluster is evaluated.
(3) Binary membership matching between each feature class point belonging either to EX-class or the QUclass, and selected ML-cluster is performed. According to that the number of each EX-or QU-cluster point and ML-clusters is evaluated showing how the feature data fit into clusters.

Catchment area
The analysis starts with the catchment area. Respective Lmoment diagrams are shown in Fig. 9. Medians for each cluster are fitting some pattern, i.e. flow gauges with the largest catchment area tend to group themselves in the leftbottom corner of the diagram, while those with small catchment area are grouped in the right-top part of the diagram. The meaning of these regions is as follows. Leftbottom part of the L-moment diagram depicts observations with low L-skewness s 3 and low L-kurtosis s 4 . L-skewness is closer to zero, which reflects the PDF to be near symmetrical. Lower L-kurtosis values depict lighter tails of the distribution.
Since L-kurtosis for normal distribution (supposed to have no fat-tails) equals to s 4 ¼ 0:1226, observed C4 category median s 4 ¼ 0:2357 means that catchment data belonging to C4 are less skewed and do not exhibit too fat tails.
Therefore, the plot brings forward the following observation. Large catchment area means more symmetrical and less tailed distribution. This fact can be simply explained by the large water accumulation and catchment area inertia. In contrast, sites with low catchment area can respond to the rainfall events more abruptly, which results in asymmetric PDF, fatter and longer tails, and finally in more persistent behavior. Finally, the ML-clusters matching analysis is done. Table 4 presents, how each catchment area category fits into clusters obtained using machine learning (ML-clusters). It is visible that the largest values, for both classifications (EX and QU) are positioned on the reverse diagonal, i.e. C1 matches ML4 and C4 ML1. Though this matching is not exact it confirms observations in a qualitative way.
An additional observation is that expert matching is not even, as for instance EX-class C1 has only 6 observations. That fact may also affect perception of obtained results.

Mean annual precipitation
The analysis continues with mean annual precipitation. Respective LMRDs are shown in Fig.10. Median values for each cluster are fitting opposite pattern to the catchment area. Flow gauges with the highest mean annual precipitation tend to group themselves in the top-right corner of the diagram, while smaller ones in the opposite side. Obtained relationship is not as explicit as in the previous case.
Large values of mean annual precipitation may cause larger flood events, skewing the distribution and generating the tail. Low precipitation has much lower probability for such extreme events.
Finally, the ML-clusters matching analysis is done. Table 5 presents the results. It is visible that the largest values, for both classifications (EX and QU) are on the diagonal, however not as clearly as in previous case. Probably getting better matching requires more exact realization of classification and/or clustering.

Average maximum and minimum temperatures
In the next step clustering by catchment-averaged annual means of daily minimum and maximum temperatures is investigated. Figure 11 shows respective diagrams for minimum and maximum temperatures and both classifications. The relationships differs between both temperatures.
Average minimum values follow visible pattern, as the lower minimum temperatures are more persistent. They exhibit fat tails and have more skewed distributions.  Average maximum temperatures are not so clear in interpretation. However, an interesting observation can be made. All classes are grouped close to each other, except 11; 12 ð Þ o C and in that region flow gauges exhibit properties closer to Gaussian, than for other temperature ranges. There is no obvious explanation to this fact, since relationship between temperature and flow regime is not straightforward but indirect. It might be that certain temperature ranges, such as ð11; 12Þ o C, are observed in some regions and other features of those regions directly affect flow levels. Definitely, this observation requires further attention. Table 6 presents the results for cluster matching. As it was already shown no pattern is visible as well for the maximum temperature, while minimum temperature slightly depicts diagonal pattern.   Figure 12 shows obtained LMRDs for the mean catchment elevation, while Table 7 respective cluster matching properties. Obtained results are in line with expectations. The higher the elevation, the more persistent, skewed and tailed the distribution is. It is natural that catchments with higher mean catchment elevation exhibit more frequent and higher rainfall events which translates into more flashy flow regime.

Mean elevation
Observations from the L-moment diagrams follow in a qualitative way the Table 7 clustering matching. Highlighted the most matching points lie within the diagonal of the table. Figure 13 shows obtained L-moment ratio diagram for the average slope. Observed results are less evident than for mean elevation. Probably catchment with low and medium slopes cannot be easily recognized on the diagrams, however high-slope catchments (typically in mountainous region) are clearly associated with distribution asymmetry and fat tails.

Average slope
Similarly to the LMRDs, the matching table is unclear. The matching seems to be fine, as the pattern is similar to the one observed with L-moments diagrams.

Land use
The analysis of the dominant land use is more complex, as it involves categorical rather than quantitative data. The following classes are used: arable land, urban, grassland and forests. Therefore, the analysis has to be performed in a different way. First, the results are evaluated similarly to previous categories. Finally, the L-moments diagram is plotted showing different types of land use, not the categories within the class (Table 8). Figure 14 shows L-moment diagrams for arable and urban land uses. Figure 15 presents grassland and forests.
Catchments are dived into four groups of different land use types in each case.
Observations confirm expectations and intuition. Low fraction of arable land use means more persistent behaviour. Possible explanation for that observation might be the result of the extensive drainage of agricultural land in Poland. An opposite situation occurs for urban land use fraction. Low urbanization allows more water retention and decreases data persistence, while higher fraction of urban area means more impervious surfaces and faster catchment response to rainfall. It is reflected by skewed distribution with fatter and longer tails.
Forests have clear impact. Firstly, all the categories are quite close to each other, so it is hard to distinguish statistical properties. L-moment diagram shows that the most natural (with the least probability of extreme incidents) is category C2 depicting the forest fraction of 20; 40 ð Þ%. Selection of the categories for grassland shows that categories C1 and C2 exhibit very close to each other, therefore meadows fraction of 0; 10 ð Þ% are characterized Cluster matching is presented in a single Table 9. Observation of these results and their comparison confirms qualitatively previous results. Figure 16 presents regional analysis of land use classes by means of L-moments diagram. Interestingly, urban land use class lies outside of the   (Table 10) is presented together with a combined L-moments diagram (see Fig. 17) showing relationships between all considered categories. Taking into account former experience with L-moment diagrams and matching properties the conclusions might be derived using the matching tables only. Soil hydrological groups A and B (typically sands and loamy sands), characterized by highest permeability, behave in a diagonal way, i.e. lower fraction of these soils (and thus a higher fraction of impermeable soils) results in higher site time series persistence, i.e. resulting in more extreme behavior.
L-moment diagram showing regional comparison between dominant soil hydrological group shows that group D does not dominate in any catchment area. It also confirms previous observations that soils with lower permeability class contribute to the higher probability of extreme events. Distributions for such catchment areas are characterized by higher skewness and kurtosis.

Soil permeability
Soil permeability is classified into five categories, i.e. very low, low, average, variable, high and very high. Similarly to the soil hydrological grouping, only Table 11 with clusters matching is presented. Observed results are in line with intuition. The lower soil permeability is, the higher persistence in distribution is noticed (higher skewness and kurtosis, i.e. fatter and longer tails are observed). In contrast high and very high soil permeability stimulate infiltration and thus protect against extreme events. Average permeability category allows to revert classification matching, which means that only high, very low and low permeable soils encourage persistence.

Fraction of lakes
Fraction of lakes as a single category is analyzed in a more simple way. Table 12 shows the comparison. Table 12 is characterized by counter diagonal relationship. Lower fraction of lakes very significantly moves the category towards ML3 and ML4, which are responsible for extreme properties (high kurtosis and skewness). As opposed, high fraction of lakes introduces accumulation and inertia into the process. Presence of lakes in a catchment results in a buffering effect on river flows, putting the PDF properties closer to the normal Gaussian distribution. This relation is also clearly confirmed by the L-moment diagrams (Fig. 18). This relationship is the strongest among all studied catchment properties. This classification is the most distinct considering the distance between the categories with respect to the difference in L-skewness and L-kurtosis. The range of variability in L-skewness s 3 is from 0.27 to 0.51, while the kurtosis is in the range of s 4 2 0:18; 0:35 ð Þ . It is noteworthy that the category C4 denoting the highest fraction of lakes results in the median being the closest to the normal distribution, out of all the considered features. Thus, presence of lakes offers the biggest positive impact protecting against river flow extremes.

Conclusions and further research
This work presents statistical regional daily river flow analysis. L-moment diagrams are used to investigate an impact of various catchment properties (catchment area, average precipitation, temperatures, slope, elevation, land use, dominant soil group, soil permeability and fraction of lakes) on the probability of flow extreme events. Decomposition of the catchment features into classes and their presentation using LMRD allows to perform clustering. Moreover, such classification is compared with direct use of the k-mean clustering.
Results show that statistical L-moments analysis can be improved with the use of data mining clustering algorithm. First of all, applied methodology works well and gives an insight into the nature of daily river flow time series for a given set of catchments. Obtained results are clear and confirm common hydrological knowledge. Such confirmation allows to use machine learning based clustering in the LMRD domain, apart from the existing approaches, as an alternative tool for the hydrological data analysis.
This work is interdisciplinary. Although the simplest clustering algorithm is applied, obtained results are in line with the hydrology-based categorization. Such results are promising. More research is required on better customization of the ML-based clustering algorithm. However, once it is fit and tuned, the hydrology analyst will get entirely new investigation methodology. Deeper coordination with the hydrological properties will allow not only to determine catchments that are most vulnerable to (a) (b) Fig. 18 L-moments diagrams with clustering by fraction of lakes extreme events, but also enable root cause analysis giving recommendations for protective actions. Sample regression analysis presented in Sect. 2 is still preliminary and it also requires further analysis,as initial investigation is promising.
Presented machine learning clustering approach is independent on data space dimension. In current research it is applied to the standard 2D L-moment ratio diagrams. However, there are no formal obstacles to conduct similar analysis in more dimension. Therefore it is planned to prepare 3D L-moment diagrams (L-Cv / L-skewness / Lkurtosis) and do clustering is such a space.
Presented results are even more valuable, in times of rapid climate change and increasingly frequent extremely dangerous flood events in Europe (Piniewski et al. 2017).