A Global-Scale Investigation of Stochastic Similarities in Marginal Distribution and Dependence Structure of Key Hydrological-Cycle Processes

: To seek stochastic analogies in key processes related to the hydrological cycle, an extended collection of several billions of records from hundred thousands of worldwide stations is used in this work. The examined processes are the near-surface hourly temperature, dew point, relative humidity, sea level pressure, and atmospheric wind speed, as well as the hourly/daily streamﬂow and precipitation. Through the use of robust stochastic metrics such as the K-moments and a second-order climacogram (i.e., variance of the averaged process vs. scale), it is found that several stochastic similarities exist in both the marginal structure, in terms of the ﬁrst four moments, and in the second-order dependence structure. Stochastic similarities are also detected among the examined processes, forming a speciﬁc hierarchy among their marginal and dependence structures, similar to the one in the hydrological cycle. Finally, similarities are also traced to the isotropic and nearly Gaussian turbulence, as analyzed through extensive lab recordings of grid turbulence and of turbulent buoyant jet along the axis, which resembles the turbulent shear and buoyant regime that dominates and drives the hydrological-cycle processes in the boundary layer. The results are found to be consistent with other studies in literature such as solar radiation, ocean waves, and evaporation, and they can be also justiﬁed by the principle of maximum entropy. Therefore, they allow for the development of a universal stochastic view of the hydrological-cycle under the Hurst–Kolmogorov dynamics, with marginal structures extending from nearly Gaussian to Pareto-type tail behavior, and with dependence structures exhibiting roughness (fractal) behavior at small scales, long-term persistence at large scales, and a transient behavior at intermediate scales.


Introduction
Stochastic approaches have been preferable over the deterministic ones when the identification and simulation of the complicated fluctuations observed in geophysical processes are the tasks at hand. In recent times, we characterize a process as complex if it is difficult to analyze or explain in a simple manner.
The origin of the word "complex" comes from Latin but has been re-borrowed from ancient Greek (originated from the verb "συµπλέκω") and is attributed to "a whole comprised of parts". It constitutes the Latin preposition "com" or "cum", which is related to the Greek preposition "συν" and is used, usually at the beginning of a word, to declare union, ensemble, etc., as well as the Latin verb "plectere", which comes from the Greek verb "πλέκω" meaning "weave", "twine", etc.
Climate dynamics is characterized by high complexity since it involves the spatiotemporal evolution of numerous geophysical variables (i.e., multivariate stochastic processes) interacting with each other in a nonlinear way, forming a hydrological cycle. Never-theless, even if we could determine a set of physical laws that describe in full detail the complexity of climate dynamics, it would be impossible to combine the equations for the purpose of predictability due to the existence of chaos, i.e., a nonpredictive sensitivity to initial conditions. For example, consider the analysis of Poincare [1] for the three-body problem, where chaotic behavior emerges from the equations of classical mechanics when studying the interacting gravitational forces between three bodies (e.g., planets). Similar results came into sight from Lorenz [2] while applying a simplified set of equations for the analysis of atmospheric dynamics. Lorenz came across the idea that nonlinear dynamic systems may have a finite limit of predictability, which for weather prediction was estimated to be around two weeks, even if the model is perfect and even if the initial conditions are known almost perfectly. Later on, numerous methodologies were initiated not for predicting the exact outcome of a nonlinear system, which as we already explained may be trivial, but rather for estimating the limits of this prediction through the alternative approach of stochastic analysis.
However, a detailed stochastic analysis of such complicated fluctuations would demand high availability of observational records at both small (e.g., hourly or even finer) and large scales (e.g., several decades)-a rare case in the hydrological-cycle processes.
In this work, we use robust stochastic metrics to increase the extracted information from time series in order to seek similarities among hydrological-cycle processes in terms of the second-order statistics. In particular, for the analysis of the marginal structure, we implement the central knowable moments (K-moments) that enable more reliable estimation from data since they substitute higher-order deviations from the mean by the use of the probability distribution function [3]. The second-order dependence structure (commonly expressed through the autocovariance function) is estimated through the climacogram (i.e., variance of the averaged process vs. scale) and other climacogram-based metrics (similar to the commonly applied autocorrelation function and power spectrum). They have all been shown to exhibit a smaller estimation bias in both large scales and in a variety of Markovian and long-term persistent (LTP) processes, as well as other advantages, compared to the aforementioned common metrics [4].
In seeking a broad and comprehensive investigation of the stochastic similarities among the hydrological-cycle processes, we implement the above metrics to a global-scale network of stations, as well as to the key processes of the hydrological cycle. To be specific, we study the near-surface hourly temperature, relative humidity, sea level pressure, and wind speed (extracted from the ISD database), as well as the hourly and daily streamflow (extracted from the USGS and CAMELS databases), hourly and daily precipitation (extracted from the HPD and the GHCN databases), and we also take into account other processes related to the hydrological cycle, such as hourly solar radiation, ocean waves and evaporation, as analyzed at global scale in other studies. In total, we extract and handle approximately 50 × 10 10 records from over 2 × 10 5 hydroclimatic global-scale stations on a daily, hourly, and subhourly resolution, with some dating even back to 1900. Specifically, we transform to hourly resolution the subhourly timeseries of air temperature, dew point, humidity, wind-speed, and sea level pressure (with around 15,000 active stations each). Similarly, we transform to hourly resolution over 10,000 timeseries of streamflow and 5,000 timeseries of precipitation located in the USA, and we merge them with over 600 USA streamflow and over 10 5 worldwide precipitation timeseries of daily resolution. It is noted that after quality control, only the 10% of records are finally selected for the analysis (Figure 1 and Table 1).  Finally, the microscopic processes driving and generating the hydrometeorological ones are governed by turbulent state. For example, the size of the drops, which is highly  Finally, the microscopic processes driving and generating the hydrometeorological ones are governed by turbulent state. For example, the size of the drops, which is highly linked to the form and intensity of precipitation events, is strongly affected by the turbulent state of small-scale atmospheric wind [5]. In addition, in a physical basis, the rain rate is found to be a function of the gradient level wind speed, the translational velocity of the tropical cyclone, the surface drag coefficient, and the average temperature and saturation ratio inside the tropical cyclone boundary layer [6]. Therefore, since the above key processes are dominated and driven by the turbulent shear and buoyancy in the nearly Gaussian atmospheric boundary layer [7][8][9], we investigate whether further stochastic similarities exist with respect to such turbulence time series, and specifically, to a huge grid-turbulence lab database (provided by the Johns Hopkins University) and to laboratory records of turbulent buoyant jets along the jet axis, both of which are found to be nearly Gaussian and isotropic.
The current study focuses on the Hurst-Kolmogorov dynamics, which characterizes a geophysical process exhibiting the Hurst phenomenon [10] (or else LTP), with an arbitrary marginal distribution function and dependence structure extending in a continuous manner from extremely fine to very large scales. The mathematical description of the Hurst phenomenon is attributed to Kolmogorov [11], who developed it while studying turbulence, and it corresponded only to a Gaussian process with a power-law decay autocorrelation function (e.g., [12]). To give credit to both contributing scientists and to incorporate alternate short-range dependence with an arbitrary marginal structure, Koutsoyiannis [13] named this general behavior as Hurst-Kolmogorov. In this work, similarities are found to exist in both (a) the marginal structure, through the mixed Hurst-Kolmogorov (e.g., Pareto-Burr-Feller for positively defined processes) probability density distribution functions, which, depending on the selected shape and scale parameters, is applicable to all the examined processes from (truncated) nearly Gaussian distributions to heavy-tail Pareto ones; and (b) the second-order dependence structure, through a Hurst-Kolmogorov generalized model, which is also expanded to include the observed curved behavior at the intermediate scales. These similarities can be well described within the framework of the maximum entropy and the Hurst-Kolmogorov dynamics (see definitions in [14][15][16][17][18]), which can be implemented through the method of moments (e.g., explicit models), through nonlinear transformation (e.g., copulas), or disaggregation (e.g., downscaling and pulse models) schemes of the Hurst-Kolmogorov dynamics, as summarized for the field of Hydrology [19,20] and beyond [3,4,[21][22][23][24][25][26][27][28][29][30][31], preserving both the marginal and dependence structures for a vast range of scales, including double periodic and intermittent behaviors (see the discussion in [18], and references therein).
In Section 2, we introduce the applied metrics and methodology, while in Section 3, we present several statistical characteristics of the selected stations of the global databases, as well as the results obtained from the analysis with a focus on the marginal and the dependence structures. Finally, in Sections 4 and 5, we summarize our findings, and we discuss how they may allow for the development of a uniting stochastic view of the hydrological-cycle processes under the Hurst-Kolmogorov dynamics, expanding from nearly Gaussian to Pareto-type tail behavior, and from fractal and intermittent behavior at small scales to LTP behavior at large scales.

Methodology
In this section, we define the central K-moments and the climacogram-based metrics that are used for the investigation of the stochastic similarities of the key hydrologicalcycle processes.

Dependence Structure Metrics
We focus on the dependence structure of second-order statistics (see the discussion in [3,4,26,32,33]), which is commonly estimated through the autocovariance function, c(h), where h is the continuous time lag in time units, or through the power spectrum, whose definition is based on the autocovariance function, i.e., s(w) := 2 ∞ −∞ c(h) cos(2πwh)dh, where w is the continuous frequency in reverse time units. A common estimator of the latter is the so-called periodogram, which is often used as a preferable estimator over the autocovariance classical estimator [34,35]. Since these two estimators are based only on the domain of lag and frequency, Beran [36] attempted to implement a similar method to the scale domain, by plotting in logarithmic axes the variance of the accumulated [37] or else aggregated [38] process vs. scale (a method often misnamed as "aggregated-variance"; however, variance is not aggregated but rather scale [39]). This attempt was inspired by the empirical work of Smith [40] in agricultural crops and by other similar works in other fields [41,42], but only to be later abandoned as a bad estimator of the longterm persistence mainly due to the large estimation bias [35,43]. The alternative estimator of the rescaled/range (R/S), which was introduced by Hurst [10] for the identification of long-term persistence, was also deemed to exhibit similar issues as the aggregatedvariance method [35,44,45] or even worse [19,46]. In addition, since a single name for this method did not exist (as, for example, for the periodogram or the correlogram methods), Koutsoyiannis [13] coined the term climacogram to emphasize the graphical representation and the link of the concept to scale (i.e., climax in Greek), so as not to be confused with the already established term of scale(o)gram. It is noted that the climacogram is explicitly linked to the autocovariance c(h), i.e., c(h) = 1 2 d 2 h 2 γ(h) /dh 2 , and thus, also to the power spectrum [32].
Finally, the identification of the stochastic structure in the scale domain is revisited much later and compared to the lag domain (i.e., through the correlogram estimator of the autocorrelation function) and the frequency domain (i.e., through the periodogram estimator of the power spectrum) but after theoretically defining the expressions, estimators, and bias expressions for all three metrics [20]. Specifically, after properly adjusting the climacogram for estimation bias ( [15,26,47,48]), and the autocovariance and power spectrum for estimation and discretization bias [4], and thus, after eliminating the main criticized limitations of the three classical estimators, Dimitriadis and Koutsoyiannis [4] show that the climacogram outperforms the other two estimators. Hence, a stochastic process with short-term fractal and long-term persistent behaviors should be preferred to be analyzed in the scale domain instead of in the lag and frequency domains.
The second-order climacogram, in terms of the variance, is defined as: where k is the scale in units of x, and similarly, the climacogram-based spectrum (CBS), or else called climacospectrum, is defined as follows [3]: The estimators of the above climacogram-based metrics can be expressed through the estimator of the climacogram, i.e., where κ = k/∆ is the dimensionless scale, ∆ is the time resolution of the process, n/κ is the integer part of n/κ, and x (κ) i is the i-th element of the averaged sample of the process at scale κ, i.e., Although the above two metrics contain the exact same information, their estimators exhibit a lower estimation bias in different range of scales. Particularly here, we use the climacogram for the identification of the dependence structure behavior at large scales (or else Hurst behavior; [10,37,43,91]), and the CBS for the investigation of small scales characterized by the fractal behavior [92][93][94][95], while for the intermediate scales, both metrics are used (e.g., see [75,76]).
Since stations may have different lengths and recording resolutions, we adjust the climacogram Var x (κ) i at resolution ∆ for the estimation bias as follows [26]: where γ(k) is the climacogram model, which for the examined processes is shown to be well described by a Hurst-Kolmogorov (HK) dynamics model, expressed as: where λ is the variance of the process; q is a scale parameter in units of the scale k that affects the behavior of the intermediate scales of the process; M is a dimensionless shape parameter that simulates the fractal behavior of the process at small scales affected by intermittency; and H is the Hurst parameter indicative of the long-term persistence of the process, i.e., for 0.5 < H < 1 the process exhibits LTP behavior, while for 0 < H < 0.5 an anti-persistent behavior, and for H = 0.5 a white-noise behavior.
Another model that is found to adequately simulate the observed convex-shaped behavior at the intermediate scales (see empirical evidence in all examined processes below) can be expressed as: where the second expression corresponds to a continuous-time Markov model. A more generalized model includes the replacement of the second Markovian model with an HK model in terms of the autocovariance function [4]. The standardized climacogram and CBS, which are used here for the identification of the dependence structure behavior, can be easily estimated by dividing them with the sample variance of the process, i.e.,γ(k)/γ(1) andζ(k)/γ(1), respectively. An in-depth analysis of other climacogram expressions for the HK dynamics can be found in [3].

Marginal Structure Metrics
The investigation of the marginal structure of the key hydrological-cycle processes is based here on the estimation of high-order moments, which can be indicative of the intermittence and tail behavior of the process (e.g., see several applications up to the 6th moment-order in [26]). Therefore, we focus our investigation on the similarities presented in the skewness and kurtosis coefficients, based on the unbiased estimators of the standardized central-moments (C-moments), as well as on the standardized L-moments, i.e., L-skewness and L-kurtosis, and the corresponding K-moments, which are described below.
Koutsoyiannis [3] has extensively illustrated that classical moments beyond the second order are unknowable, especially in the presence of LTP, in the sense that their estimation from data is not feasible. To tackle this issue, the so-called (knowable) K-moments were introduced, in which the higher-order estimations are implemented by the use of the probability distribution function. In particular, the central K-moments that are used for estimation here, are invariant to a shift of origin and they are homogeneous in terms of multiplication by a scalar, and thus, they enable more reliable estimations. The central K-moments are defined as: Koutsoyiannis [86] also introduced the hyper-central K-moments, which can be estimated directly from the central K-moments, but they present several advantages over the latter in model fitting and inference from data, and can be defined as: for p ≥ q, where p and q reflect the order of the estimation, x is the stochastic process of interest with marginal probability distribution F(x) (it is noted that underlined quantities denote random variables), and µ is the process mean. The hyper-central K-skewness and K-kurtosis that are used here for illustrating the behavior of the marginal structure can be expressed through the central K-moments, respectively, as: and A sample estimator of the central K-moments for q = 2, which is used here and is approximately adjusted for bias, is derived as: where n is the length of the sample,μ = x i /n is the estimator of the mean, and x (i:n) is the sample of the process rearranged in ascending order. For comparison, we show results from the classical unbiased estimators of the skewness and kurtosis coefficients, denoted here as C-skewness and C-kurtosis, respectively. Moreover, we show results for the classical L-moments by Hosking [96], but they are estimated here through the hyper-central K-moments for comparison, respectively, as: and where λ 2 , λ 3 , and λ 4 , are the well-known L-moments of the 2nd, 3rd, and 4th order, respectively. We note that the commonly used L-kurtosis, i.e., λ 4 /λ 2 , is linearly related to the modified version of the L-kurtosis shown in the latter equation.
The investigation here focuses on the marginal moments instead of the probability distribution itself due to the large estimation bias entailed in the estimators of the latter. Instead of directly fitting a distribution function to data, one may estimate a small number of moments from data equal to the number of parameters in the target probability distribution function, and then estimate the parameters of the distribution from the values of the moments. A flexible probability distribution density function that can simulate processes with a lower threshold, such as the ones examined here, is the truncated-mixed Pareto-Burr-Feller (PBF) distribution, as named after (a) the engineer Pareto, who discovered the family of power-type distributions while working on the size distribution of incomes in a society [97,98]; (b) to Burr [99] who identified and analyzed this distribution family (but without giving a justification), first proposed as an algebraic form by Bierens de Haan; and (c) to Feller [100] who linked it to the Beta function and distribution through a linear power transformation, which was further analyzed and summarized by Arnold and Press [101]. The truncated-mixed PBF distribution is used here to test whether its limits on the skewness-kurtosis frame plot can adequately simulate the observed values estimated through the K-moments, by including a nonzero probability at zero and a lower-truncation value x > d. The mixed PBF can be expressed as follows (e.g., [3]): where ζ and ξ are shape parameters, λ is a scale parameter, and P 1 = 1 − P{x = 0}. Then, the truncated mixed PBF distribution is constructed by dividing the above expressions by F(x u ) − F(x l ), where x u and x l are the upper and lower limits of the truncation. Specifically, x l = 0 for temperature (if expressed in Kelvin degrees), dew point, humidity, wind speed, pressure, streamflow, and precipitation, and x ul = 1 for humidity, while for the dew point x u is set at the temperature value, and for the rest of the processes at potential infinity.
The PBF distribution can also be viewed from the perspective of the generalized HK model defined in the previous section, i.e., f ( is a scale parameter, b and c are shape parameters, d is a position parameter, and λ can be estimated so that f (x) = 1. Similarly to the PBF, this expression contains a variety of distributions from nearly Gaussian (e.g., for ( where µ and σ are the mean and standard deviation of the process)such as the temperature, atmospheric pressure, or grid turbulence and turbulent buoyant jet processes-to the Pareto-tail distribution (e.g., for (|x − d|/a) 2b >> 1, a = p 2 p 1 /(p 1 +1) , c = (1 − p 1 )/2, and λ = p 1 , where p 1 and p 2 are the shape and scale parameters of the Pareto distribution)-such as the wind speed, streamflow, precipitation, and ocean wave processes-while it can be easily truncated between absolute zero and temperature for the dew point, between solar radiation on the ground and on the top of the atmosphere for the clearness index, and between 0 and 1 for the relative humidity processes.

Global-Scale Data Extraction and Processing
The selected global-scale databases for the analysis after quality control are the ISD database [102][103][104][105][106][107][108][109][110][111] for near-surface hourly temperature, dew point, sea level pressure, and wind speed; the CAMELS database [112][113][114][115] for the daily streamflow; the USGS database [116] for the hourly streamflow; the GHCN database [117][118][119][120][121] for the daily precipitation; and the HPD database (https://www.ncei.noaa.gov/, lastly accessed on 15 December 2020) for the hourly precipitation. From the contained stations in each database, we selected stations with more than 5 years of full records (see Table 1 for more information on the selected stations, and Figure 1 for the visualization of their locations), and we discarded recordings that have quality flags. In addition, for other streamflow-related processes, we analyzed and discussed one of the largest daily discharge time series of the river Po [122] and the largest annual time series of the stage at the river Nile (e.g., [123]).
The relative humidity is estimated from the hourly temperature and dew point values as [124]: where U is the relative humidity (%); T and T d are the near-surface temperature and the dew point expressed in Kelvin degrees, respectively; and T 0 = 273.16 K is the triple point of water.
To mitigate the effect that the periodicity of hydrological-cycle processes, prominent both in the diurnal and seasonal cycles (e.g., [125][126][127][128][129][130][131][132]), exerts on their modelling, we apply a double standardization on the processes with hourly resolution and a seasonal standardization on the ones with daily resolution. In particular, we subtract the mean from each periodicity cycle, and we divide with its standard deviation. In this way, although the effect of the double periodicity cannot be entirely removed-since this is only possible for Gaussian-distributed variables-we diminish the propagation of its effect on the dependence structure along the subsequent months and hours. Moreover, to increase the available records in each cycle and the estimation accuracy, and without much affecting the dependence structure, we standardize each process based on a 6 hourly × 3 monthly resolution of the cycles, i.e., the diurnal periodicity is divided into 24/6 = 4 cycles, and the seasonal periodicity is divided also in 12/3 = 4 cycles.
In addition, the marginal and dependence structures of other key hydrological-cycle processes were included in the analysis, and the results from previous studies of globalscale stations or from the longest available samples worldwide were discussed. Specifically, these processes are the potential evapotranspiration [133,134], the solar radiation clearness index [62], and the ocean wind-wave height and period processes [135].
Finally, for the investigation of the stochastic similarities between the above key hydrological-cycle processes and the isotropic and nearly Gaussian turbulence, which resemble the turbulent shear and buoyancy in the atmospheric boundary layer, we show the results from a massive grid-turbulence database that includes 40 time series, each with 36 × 10 6 records of longitudinal wind velocity along the flow direction, all measured by X-wire probes placed downstream of the grid and in different positions, and with a sampling temporal resolution of 25 µs [136]. To shift from the spatial to the temporal domain [137], we apply a standardization to all time series by subtracting their mean and dividing by their standard deviation (see more information and results in [76]). In this manner, we may directly estimate the expected marginal and temporal dependence structure by combining the estimations from all the time series, approximately as if the same experiment was performed multiple times in the same position. For the buoyancy behavior, we discuss the results from several studies of Papanicolaou and List [138,139] and Dimitriadis et al. [18,140,141], where more than 10 time series of horizontal and vertical positively buoyant thermal jets of temperature concentration, recorded with the laserinduced-fluorescence technique, and with a 30 ms resolution, various nozzle diameters, discharges, initial and ambient temperature, and of more than 10 4 sample length each, were analyzed.

Results
The investigation of the marginal structure of the key hydrological-cycle processes is based here on the estimation of moments, and especially of high-order, which are indicative of the tail behavior. In this section, we present the results from the global-scale analysis, in terms of the second-order statistics, and particularly, the plots of the mean vs. standard deviation, of the skewness vs. kurtosis estimated from the C-, L-, K-moments, and of the climacogram and the climacospectrum (mean, and the 5% and 95% quantiles). Note that for the estimation of the model parameters adjusting for bias, we use up to the 10% of the maximum scales of the climacogram (and thus, 5% of the CBS), in order not to use variance estimates from a sample length less than 10% of the original length [4].
In Figures 2-8, which show the climacogram and climacospectrum for each process, we observed similarities in the shape of the dependence structure spanning from strong correlations at the small scales to a power-type behavior at large scales, with a convex-shape at the intermediate scales. al. [18,140,141], where more than 10 time series of horizontal and vertical positively buoyant thermal jets of temperature concentration, recorded with the laser-induced-fluorescence technique, and with a 30 ms resolution, various nozzle diameters, discharges, initial and ambient temperature, and of more than 10 4 sample length each, were analyzed.

Results
The investigation of the marginal structure of the key hydrological-cycle processes is based here on the estimation of moments, and especially of high-order, which are indicative of the tail behavior. In this section, we present the results from the global-scale analysis, in terms of the second-order statistics, and particularly, the plots of the mean vs. standard deviation, of the skewness vs. kurtosis estimated from the C-, L-, K-moments, and of the climacogram and the climacospectrum (mean, and the 5% and 95% quantiles). Note that for the estimation of the model parameters adjusting for bias, we use up to the 10% of the maximum scales of the climacogram (and thus, 5% of the CBS), in order not to use variance estimates from a sample length less than 10% of the original length [4].
In Figures 2-8, which show the climacogram and climacospectrum for each process, we observed similarities in the shape of the dependence structure spanning from strong correlations at the small scales to a power-type behavior at large scales, with a convexshape at the intermediate scales.               Although similarities may also be traced in the skewness-kurtosis plots, the more robust estimation of the L-skewness vs. L-kurtosis (Figures 9 and 10) allows for more powerful comparisons. Still, the more advanced K-skewness and K-kurtosis estimations enable an even clearer view of the detected pattern in Figures 9 and 11. Specifically, the K-skewness increases from the dew point to the humidity, to the temperature, sea level pressure, and grid turbulence, followed by a weaker increase of the K-kurtosis, while for the wind speed, streamflow, and precipitation processes, both the K-skewness and K-kurtosis highly increase. In addition, we observe how all sample K-skewness and K-kurtosis fall into the range of the mixed PBF distribution only by the analysis through K-moments. Finally, in Figures 12 and 13, we observed that the dew point, relative humidity, temperature, and sea level pressure, along with grid turbulence, exhibit a stronger structure at the intermediate and large scales, followed by the wind speed, streamflow, and precipitation.
Therefore, this hierarchy of the key hydrological-cycle processes may be described by a transformation of the (truncated) nearly Gaussian processes (i.e., dew point, humidity, temperature, and sea level pressure) with a stronger LTP dependence at large scales (H > 0. 75) to the Pareto-type processes (i.e., wind speed, streamflow, and precipitation) with a weaker dependence LTP at larger scales (H ≤ 0.75). Interestingly, the same hierarchy was observed in a simplified cycle of the energy exchange among hydrological processes trough turbulent mixing, with the dew point, relative humidity, temperature, and sea level pressure, feeding the wind speed, while triggering precipitation, whose energy is then temporally stored in rivers and soil, and finally returned to the former processes through evaporation.
For the rest of the processes included in the current analysis, it was found that the clearness index of the solar radiation index [62] and the ocean waves (height and period; [135]), as analyzed from the global databases mentioned above, exhibit a marginal structure similar to the sea level pressure and wind speed, respectively. In addition, they exhibit an HK dependence structure with H ≈ 0.8 and q ≈ 30 h, and H ≈ 0.9 and q ≈ 10 h, for the wave height and period, respectively, and H = 0.83 and q = 2 h for the clearness index of the solar radiation, while both have a fractal behavior with roughness (M < 0.5). The standardized temperature concentration of the turbulent buoyant jet along the axis [141], is found to exhibit similar behavior as in the grid turbulence with a nearly Gaussian marginal structure, and a rough (M < 0.5, at small scales) and strong LTP (H > 0.6, at the jet-like area and H > 0.9 at the plume-like area) dependence structure. Finally, the evapotranspiration process was also found to exhibit a weak LTP behavior (H ≈ 0.6) and a similar marginal structure to temperature and solar radiation [133,134]. Although similarities may also be traced in the skewness-kurtosis plots, the more robust estimation of the L-skewness vs. L-kurtosis (Figures 9 and 10) allows for more powerful comparisons. Still, the more advanced K-skewness and K-kurtosis estimations enable an even clearer view of the detected pattern in Figures 9 and 11. Specifically, the Kskewness increases from the dew point to the humidity, to the temperature, sea level pressure, and grid turbulence, followed by a weaker increase of the K-kurtosis, while for the wind speed, streamflow, and precipitation processes, both the K-skewness and K-kurtosis highly increase. In addition, we observe how all sample K-skewness and K-kurtosis fall into the range of the mixed PBF distribution only by the analysis through K-moments. Finally, in Figures 12 and 13, we observed that the dew point, relative humidity, temperature, and sea level pressure, along with grid turbulence, exhibit a stronger structure at the intermediate and large scales, followed by the wind speed, streamflow, and precipitation.
Therefore, this hierarchy of the key hydrological-cycle processes may be described by a transformation of the (truncated) nearly Gaussian processes (i.e., dew point, humidity, temperature, and sea level pressure) with a stronger LTP dependence at large scales (H > 0.75) to the Pareto-type processes (i.e., wind speed, streamflow, and precipitation) with a weaker dependence LTP at larger scales (H ≤ 0.75). Interestingly, the same hierarchy was observed in a simplified cycle of the energy exchange among hydrological processes trough turbulent mixing, with the dew point, relative humidity, temperature, and sea level pressure, feeding the wind speed, while triggering precipitation, whose energy is then temporally stored in rivers and soil, and finally returned to the former processes through evaporation. Figure 9. L-skewness vs. L-kurtosis, and K-skewness vs. K-kurtosis estimated through the hyper-central K-moments, for the key hydrological-cycle and the grid-turbulence processes.        For the rest of the processes included in the current analysis, it was found that the clearness index of the solar radiation index [62] and the ocean waves (height and period; [135]), as analyzed from the global databases mentioned above, exhibit a marginal structure similar to the sea level pressure and wind speed, respectively. In addition, they exhibit an HK dependence structure with H ≈ 0.8 and q ≈ 30 h, and H ≈ 0.9 and q ≈ 10 h, for the wave height and period, respectively, and H = 0.83 and q = 2 h for the clearness index of the solar radiation, while both have a fractal behavior with roughness (M < 0.5). The standardized temperature concentration of the turbulent buoyant jet along the axis [141], is found to exhibit similar behavior as in the grid turbulence with a nearly Gaussian marginal structure, and a rough (M < 0.5, at small scales) and strong LTP (H > 0.6, at the jetlike area and H > 0.9 at the plume-like area) dependence structure. Finally, the evapotranspiration process was also found to exhibit a weak LTP behavior (H ≈ 0.6) and a similar marginal structure to temperature and solar radiation [133,134].   For the rest of the processes included in the current analysis, it was found that the clearness index of the solar radiation index [62] and the ocean waves (height and period; [135]), as analyzed from the global databases mentioned above, exhibit a marginal structure similar to the sea level pressure and wind speed, respectively. In addition, they exhibit an HK dependence structure with H ≈ 0.8 and q ≈ 30 h, and H ≈ 0.9 and q ≈ 10 h, for the wave height and period, respectively, and H = 0.83 and q = 2 h for the clearness index of the solar radiation, while both have a fractal behavior with roughness (M < 0.5). The standardized temperature concentration of the turbulent buoyant jet along the axis [141], is found to exhibit similar behavior as in the grid turbulence with a nearly Gaussian marginal structure, and a rough (M < 0.5, at small scales) and strong LTP (H > 0.6, at the jetlike area and H > 0.9 at the plume-like area) dependence structure. Finally, the evapotranspiration process was also found to exhibit a weak LTP behavior (H ≈ 0.6) and a similar marginal structure to temperature and solar radiation [133,134]. Finally, the summary statistics from the global-scale analysis can be seen in Tables 2  and 3, where the mean from the C-, L-, and K-moments as well as the scale, fractal, Hurst, and scale parameter of the expected second-order dependence structure, respectively, were estimated and are discussed in more detail at the next section. It is noted that the merge between the time series of a different resolution, such as in streamflow and precipitation, was performed after adjusting for bias (see more details from a similar global-scale analysis of fine-resolution and large-scale temperature and wind in [70]).

Discussion
The investigation of the uncertainty in the hydrological cycle is an important scientific field, as recognized by the International Association of Hydrological Sciences (IAHS) by launching the Panta Rhei research initiative for the Scientific Decade 2013-2022 [142].
An overall conclusion from this study is that the Hurst parameter is estimated significantly above 0.5 for all processes (Table 2 and Figures 12 and 13), indicating that the observed uncertainty and climatic variability in the hydrological cycle may be caused by the presence of the long-term persistent (LTP) behavior (see also discussion in [143,144]). This is consistent with the universality of LTP behavior as confirmed in various other studies (than the ones mentioned in Section 2.1) and statistical attributes in literature (for a review see [18,[145][146][147]), such as in global-scale key hydrometeorological processes [3,18], trend analysis [148] and extremes [149,150], precipitation [151][152][153][154][155], streamflow [155][156][157], turbulent jets [158,159] and grid turbulence [160,161].
The intermediate scale behavior of the dependence structure (Table 2 and Figure 12), which was found to be consistent with the K41 law of Kolmogorov [162,163] for all processes, corresponds to a fractal parameter of M = 1/3 (i.e., roughness behavior of M < 0.5). For even smaller scales, the fractal parameter was estimated even lower than 0.5 or close-tozero for all the examined processes. However, for the more robust estimation of the fractal parameter, additional data and of higher resolution are required. For example, in [70], where high resolution samples (of 10 Hz) are applied to the temperature and wind speed, the fractal parameter was estimated in both processes as M ≈ 1/3. It is noted that in the multifractal analysis, the change of the dependence structure is viewed as a scale break, and similar positions of the break is identified in the intermediate scales of streamflow and precipitation [164][165][166][167][168][169].
The high uncertainty of geophysical dynamics is linked to the power-law type of the marginal distribution as well as of the dependence structure through empirical evidence and physical justification. Although the above and other studies have focused on one or a limited number of processes, in this study we analyze several key hydrological-cycle processes, namely near-surface temperature, dew point, humidity, sea level pressure, atmospheric wind speed, streamflow, precipitation, as well as other processes from previous studies, such as solar radiation clearness index and ocean waves, where we trace stochastic similarities in their marginal and dependence structures. Moreover, we find similar stochastic structures in turbulent shear and buoyancy processes, as studied through laboratory records of grid-turbulent wind speeds and temperature concentrations of buoyant turbulent jets along the axis.
Specifically, a hierarchy emerges (a) for the marginal structure visualized through the skewness-kurtosis plot, both estimated through the K-moments, while similar but weaker empirical conclusions can be derived by the plot of the C-moments and L-moments; and (b) for the dependence structure visualized through the climacogram and climacospectrum. This hierarchy, related to the hydrological cycle, starts with temperature, dew point, relative humidity, solar radiation index, evaporation, and sea level pressure, all of which exhibit a stronger skewness over kurtosis absolute ratio than the turbulent processes, wind speed, and ocean waves. All the latter processes exhibit a stronger LTP behavior in the dependence structure, whereas streamflow and precipitation present a weaker skewness over kurtosis absolute ratio and an LTP behavior.
Interestingly, the same hierarchy is observed in the energy exchange among processes in the hydrological cycle. It starts with the solar radiation, temperature, dew point (or equivalently, relative humidity), and sea level pressure. These feed the wind speed and ocean waves through nearly Gaussian and isotropic turbulent mixing in the boundary layer. Next, they trigger precipitation, which then moves to streams and soil, and is finally returned to the former processes through evaporation.

Conclusions
The major innovation of this study is the uniting view of the key hydrological-cycle processes through the analysis of several billions of observations from hundred thousands of stations by robust statistical metrics of (a) the K-moments, for the estimation of the marginal structure of the first four moments, and of (b) the climacogram, for the estimation of the second-order dependence structure. The key examined hydrological-cycle processes are the near-surface temperature, dew point, humidity, sea level pressure, atmospheric wind speed, streamflow, and precipitation, as well as other processes from previous studies, such as shear and buoyant turbulent processes analyzed through small-scale laboratory experiments, and solar radiation, and ocean waves. The main traced stochastic similarities are as follows: (1) A hierarchy related to the hydrological cycle was identified with the dew point, temperature, relative humidity, solar radiation, and sea level pressure all exhibiting a lower skewness over kurtosis absolute ratio than the turbulent processes, wind speed, and ocean waves, and with a stronger long-term persistence (LTP) behavior in the dependence structure (H > 0.75), followed by streamflow and precipitation, both of which exhibit a smaller skewness-kurtosis absolute ratio and a weaker LTP behavior (H ≤ 0.75). (2) All the examined processes can be adequately simulated by the truncated mixed-PBF distribution, adjusting for probability dry and lower (or upper) truncation, in terms of the first four moments, and ranging from (truncated) nearly Gaussian to Pareto-type tails. (3) As the sample size increases, different records of the same process from several locations converge to a smaller area of the nondimensionalized statistics (skewnesskurtosis), indicating a common marginal behavior. (4) All the examined hydrological-cycle processes exhibit a similar dependence structure that extends from the fractal behavior with roughness (M < 0.5) located at the smallintermittent scales to the LTP behavior at large scales (H > 0.5), while both indicate large uncertainty and high climatic variability.