Assessment of Multifractal Fingerprints of Reference Evapotranspiration Based on Multivariate Empirical Mode Decomposition

Abstract

This study analyzed the multifractal characteristics of daily reference evapotranspiration (ET_o) time series of the Tabriz and Urmia stations of northwestern Iran and its cross-correlation with five other meteorological variables. The results of multifractal detrended fluctuation analysis (MFDFA) of ET_o, temperature, pressure, relative humidity solar radiation, and wind velocity showed that all the time series of both stations exhibited multifractality and long memory persistence with higher persistence and complexity in the datasets of Urmia station. Then, a multivariate empirical mode decomposition (MEMD)-(MFDFA) coupled framework was proposed to identify the dominant modes suitable for the forecasting of the different variables. The examination of reconstructed time series consistently displayed an increase in persistence and multifractality. The cross-correlation between different candidate variables and ET_o was examined using a recently proposed multifractal cross-correlation analysis (MFCCA) method. The results showed that in each pair-wise cross-correlation analysis, the joint persistence is approximately half of the persistence of an individual time series, reinforcing the universality in the fractal cross-correlation analysis. The cross-correlation properties displayed diverse patterns in different pair-wise combinations of cross-correlation analysis despite the similarity of patterns among the data of the two stations.

1. Introduction

Evapotranspiration (ET) is one of the prominent components of a hydrologic cycle, and its accurate prediction is vital in irrigation and water resources management. The rate of ET from a reference crop of 0.12 m height, bulk surface resistance of 70 sm⁻¹, and albedo of 0.23, assuming a surplus of soil moisture, is called reference evapotranspiration (ET_o) [1], is perhaps the most tangible form of evapotranspiration in the analysis of a hydrologic cycle. ET_o is often controlled by numerous meteorological variables, and the inherent complexities of such variables often force researchers to search for diverse methods, ranging from statistical to machine learning, for accurate modeling [2,3]. Hydro-climatic time series often possess nonlinear and nonstationary characteristics with the presence of trends, shifts, and abrupt fluctuations, which are also exhibited in the ET_o time series. The characterization of a time series and the examination of its long-term dependency is crucial for the accurate modeling of the process [4]. Hurst proposed the rescaled range (R/S) analysis as one of the initial contributions to the persistence of a time series [4]. Mantelbrot [5] presented the theory of fractality, which determines that on subdividing a geometrical object into parts, each part will be a downscaled copy of the object. Fractality is closely linked with persistence, and a considerable number of algorithms have been proposed for estimating the fractal behavior of the time series. The Fourier spectral method [6], wavelet-based analysis [7], arbitrary order Hilbert spectral analysis (AOHSA) [8], structure function [9,10], and visibility graph analysis [11,12] are among them. Peng et al. [13] proposed a robust technique, namely detrended fluctuation analysis (DFA), for analyzing the fractality based on a detrending operation. Kantelhardt et al. [14] propounded the multifractal (MF) extension of DFA (so-called MFDFA). The multifractal parameters can be considered as ‘fingerprints’ of the series, as they are useful in modeling [15].

In the past two decades, numerous studies on the multifractal analysis of hydro-climatic data were performed using MFDFA [16,17,18,19,20,21,22]. Multifractal properties were also detected in meteorological time series, such as temperature, wind velocity, solar radiation, and relative humidity, in different parts of the world [23,24,25,26,27,28]. Baranowski et al. [29] analyzed the scaling properties of the time series of temperature, wind speed, relative humidity, rainfall, and solar radiation comprising over 11,000 points from stations in Germany, Finland, Spain, and Poland. Long-range correlations were the reasons for the multifractality of the majority of the series, except for rainfall. Krzyszczak et al. [27] used MFDFA to analyze the multifractal properties of the meteorological datasets of four stations in Poland and Bulgaria located in varying climatic zones. They have analyzed daily time series comprising more than 11,000 observations of temperature, rainfall, relative humidity, and wind speed. It was revealed that the multifractal characteristics of rainfall are quite different from those of others. Entire time series were partitioned at about 2001–2002, and the multifractal characteristics were analyzed to evaluate the differences in dynamics. The multifractal spectrum was found to be susceptible to climatic shift, and it was especially apparent for asymmetry, which changed from being right to left skewed, implying more extreme events in the recent past.

Evapotranspiration is a complex process that happens as a visible natural outcome of the concurrent interactions of different meteorological variables. The modeling of any hydrologic variable is a challenging problem for a hydrologist, for which either a time-series approach, a cause–effect approach, or a combined approach involving the two can be followed. Several methods have been proposed in the past for modeling reference-crop ET_o including mass transfer, water budget, pan evaporation, radiation, and temperature-based methods, as well as the recent practices of data-driven methods [30,31]. Recently, Sreedevi et al. [32] investigated the teleconnection between the meteorological variables and reference ET_o of Urmia and Tabriz of northwestern Iran in a multi-scale perspective. The region is hydro-climatologically sensitive and of ecological importance because of the presence of Lake Urmia, and it is important to investigate the impact of different climate regimes on the modeling of hydrological variables. Identifying the most appropriate predictor variables and understanding the joint association between the influential variables with Et_o are among the key steps in the more popular cause–effect approach of modeling evapotranspiration. According to the authors’ knowledge, even though many studies have been proposed for evaluating the fractal properties of hydrological variables employing the MFDFA procedure, the study on the multifractal characterization of evapotranspiration is scarce in the literature [2,20,33]. Ariza-Villaverde et al. [34] used joint multifractal analysis for examining the association between temperature, relative humidity, and ET_o of Cordoba city in Spain. Multifractal characterization of daily ET_o in different climate zones of China was performed by Zhan et al. [35]. Gómez-Gómez et al. [36] analyzed the cross-correlation between the meteorological variables and ET_o of the Guadalquivir River Valley using the multifractal detrended cross-correlation method.

The objectives of this paper include (i) analyzing the multifractality daily ET_o time series and five other climate variables (temperature, pressure, relative humidity, wind velocity, and solar radiation) of two stations located in northwestern Iran, (ii) proposing an MEMD-MFDFA method for the identification of the dominant modes for the prediction of candidate variables, and (iii) investigating the cross-correlation between ET_o and the five climate variables using multifractal cross-correlation analysis (MFCCA).

2. Study Area and Observation Data

In this study, temperature (T) (°C), pressure (P) (kPa), wind velocity (U) (ms⁻¹), relative humidity (RH) (%), solar radiation (SR) (kJ m⁻²), and ET_o (mm) at daily time scales obtained from the two stations at Tabriz and Urmia cities were considered for the multifractal analysis (June 1992 to October 2018 for Tabriz and June 1992 to October 2016 for Urmia). The positions of the stations are given in Figure 1. Tabriz city, which has a mean altitude of 1350 m, is located between latitudes 38°8′ N and 46°15′ E at the intersection of the Aji and Quri rivers. Tabriz has a rainfall of about 378 mm annually; it experiences a dry to semi-hot summer and a mild climate in spring. The mean temperature is nearly 14 °C per year with a rate of ET around 1740 mm/year. Urmia is located at 37°34′ N and 44°58′ E in northwestern Iran at an average altitude of 1332 m. Urmia experiences scanty summer rainfall and heavy rainfall in winter. In most years, freezing occurs over four months. The annual rainfall of the station is nearly 310 mm [37]. Time-series plots of the data for Tabriz and Urmia used in this study are shown in Figure 2 and Figure 3.

The quick reception of solar irradiance results in the rise in temperature and the subsequent increase in ET. Therefore, Tabriz records a higher maximum pan evaporation and maximum temperature (15 mm/day and 35 °C, respectively) than Urmia (11 mm/day and 34 °C). However, because of Lake Urmia, evaporation will be higher in that region, leading to higher maximum and average daily rainfall (3.8 mm and 0.7 mm, respectively) than in Tabriz (3 mm and 0.6 mm). The statistical properties of the hydro-meteorological data of the Tabriz and Urmia stations are stated in Table 1.

3. Materials and Methods

3.1. Multifractal Detrended Fluctuation Analysis (MFDFA)

MFDFA is a well-established algorithm for detecting the scaling behaviors and multifractality of nonlinear and nonstationary datasets. This method involves (i) the computation of a profile by subtracting each value in the time series from its mean, (ii) the division of the profiles into segments of different lengths (so-called scale, s) and the polynomial fitting for each of them, (iii) estimating the local trend by taking the difference of the segment from the fitted polynomial, for each moment order (q), by the aggregation and averaging of operations, and (v) performing a logarithmic fit between the fluctuation function to obtain the generalized Hurst exponents (GHE) for different moment orders, and this value for moment order q = 2 is quite similar to the classical Hurst exponent, which can provide information on the persistence of the time series. Then, based on the plot between GHE and the q-order, one can comment on the fractality of the time series, i.e., whether it is monofractal or multifractal. Subsequently, the multifractal spectrum can be derived, which can provide useful insights into the multifractal strength of the series. More details on the algorithm are provided in Appendix A.

3.2. Multifractal Cross-Correlation Analysis (MFCCA)

The mutual association between different hydrologic and climatic variables was traditionally analyzed using a correlation coefficient (R). However, the R value may be misleading if outliers are present and it is only deciphering the linear association between the variables. Podobnik and Stanley [38] propounded the detrended cross-correlation analysis (DCCA) to analyze the power-law cross-correlations between the two signals. Extension of DCCA to a multifractal domain was made by researchers via the proposal of multifractal DCCA [39], multifractal detrending moving-average cross-correlation analysis (MFXDMA) [40], and multifractal cross-correlation analysis (MFCCA) [41]. These improvisations have been applied to time series from different scientific fields [2,20,33,42,43]. The steps of MFCCA are as follows:

• Consider two time-series signals, x_i and y_i (i = 1,2,…,N), to find their profiles:

  (1)$Xp(j)={\displaystyle \sum _{i=1}^j\left[x_i-x_m\right]}$

  and

  (2)$Yp(j)={\displaystyle \sum _{i=1}^j\left[y_i-y_m\right]}$

  where i = 1,2,……., N; x_m and y_m are the mean of the respective signals.

• Divide the series into N_s independent windows, both in progressive and retrograde orders (hence, 2Ns), to avoid any exclusion of data points at the head and tail ends.

• For each window, compute the local trends by fitting an m-order polynomial:

  (3)${ff}_{XY}{}^2(v,s)=\left\{\frac{1}{s}\sum _{k=i}^s[\left(X((v-1)+k)-p_{X,\nu }{}^m\left(k\right)\right)\times \left(Y((v-1)+k)-p_{Y,\nu }{}^m\left(k\right)\right)]\right\}$

• Estimate the detrended covariance:

  (4)$F^q{}_{XY}(s)=\frac{1}{2N_s}{\displaystyle \sum _{\upsilon =0}^{2N_s-1}sign\left[f{f}_{XY}{}^2(\upsilon ,s)\right]}{\left|f{f}_{XY}{}^2(\upsilon ,s)\right|}^{q/2}$

  where F^q_XY(s) is the fluctuation function (FFn), bearing the relation:

  F^q_XY(s) ~ s^λ(q)(5)

  where s is the segment size.

The exponent λ(q) is similar to the GHE in MFDFA. Unlike DCCA, it is a method that also accounts for the sign of FFn of different statistical moment orders. The modifications of Equations (3) and (4) for backward computations and q = 0 can be found in the literature [44].

The coefficient ρ_XY is coined as the ratio between the covariance function $F_{XY}$ and the variance functions $F_X$ and $F_Y$ after detrending [45]:

(6)${\rho }_{XY}=\frac{F^q{}_{XY}}{\sqrt{F^q{}_X}\sqrt{F^q{}_Y}}$

ρ_XY range [−1,1], (±0.66 to ±1) indicates high, (±0.33 to ±0.66) is moderate, and (zero to ±0.33) is low [46]. The MFCCA analysis can provide a better physical insight, as it indicates the scale-specific correlation between two variables, which computes the covariance function of two time series along with the signs [47,48].

3.3. Multivariate Empirical Mode Decomposition (MEMD)

Empirical mode decomposition (EMD) is a data-based decomposition technique provided by Nordan E. Huang in 1998, which decomposes a signal into enough modes with definite periodicity and is based on spline fitting through the extremes (crest and trough points) of the signal. The multivariate EMD propounded by Rehman and Mandic [49] is an extended version of EMD, which can decompose more than three signals concurrently by tracking the common scales in them. Here, multiple envelopes are developed by finding the projections of inputs along various directions in an n-D space. Let $V(t)=\left\{V_1(t),V_2(t)\dots V_n(t)\right\}$ be the m variable signals and $D{V}^{{\theta }_d}=\left\{D_1{}^d,D_2{}^d,\dots ,D{n}^d\right\}$ denote the direction vectors (DVs) defined by angles ${\theta }_d=\left\{{\theta }_1{}^d,{\theta }_2{}^d,\dots {\theta }_{n-1}{}^d\right\}$ in direction set DV (d = 1,2,3,…D), where D is the number of directions. It is noteworthy to state that instead of there being oscillatory modes in EMD, rotational modes are generated in MEMD.

The mean zero rotation components (called intrinsic mode functions or IMFs) can be obtained by the following steps:

• Generating appropriate sets of DVs by sampling on a (n − 1) unit hypersphere;

• Computing the projections $p^{{\theta }_d}(t)$ of V(t) along the DVs $D{V}^{{\theta }_d}$ for all d;

• Finding the time instants $t_i{}^{{\theta }_d}$ of the maxima of projections for all d;

• Interpolating [$t_i{}^{{\theta }_d},D(t_i{}^{{\theta }_d})$] to get the surface $E^{{\theta }_d(t)}$ for all d;

• Finding the mean of surfaces:

  (7)$S(t)=\frac{1}{K}{\displaystyle \sum _{k=1}^K{E}^{{\theta }_d}(t)}$

• Extracting R(t) = V(t) − S(t). If R(t) satisfies the termination criteria, repeat from (1) onward upon (V(t) − R(t)), otherwise repeat from (2) upon the reminder R(t).

This method has many advantages over conventional EMD, as it provides cleaner decomposition and spectral representation of the signals, ensures an equal number of components for different candidate series by recognizing the scales commonly present, and is easier to implement [49,50,51]. The Hammersley sampling scheme is suitable for generating DVs, and the Cauchy-type termination criteria is suitable for implementing the algorithm [51]. MEMD is used as a data-processing tool in a hybrid modeling framework for modeling hydrological variables including evapotranspiration [52,53,54]. This study proposes the use of MEMD instead of conventional EMD, for identifying the dominant modes, which can improve the predictability of the time series [55].

The MEMD algorithm can decompose any complex time series into a definite number of IMFs of specific frequency and end residual. In general, the modes are arranged in descending order of their frequencies, with the mode of highest frequency IMF at the top. In the practical hydrological signals, there is a good chance that some of the modes, like high-frequency IMFs, may contain noise, which may affect the performance in predictions. Hence, the elimination of the noise-intruded component is very important in improving the predictability of the hydrological signals.

3.4. MEMD-MFDFA Framework for Reconstruction

An MFDFA-based approach is proposed for identifying the less-contributing signals. The different steps of this framework are presented in Figure 4.

Most of the EMD-based de-noising methods in the past considered the highest frequency mode only (IMF1) as noise [55], and in this procedure, we consider both Gaussian white noise and anti-persistent processes by fixing a threshold H of 0.5.

4. Results and Discussion

In this section, the results of the MFDFA analysis performed on different time-series data collected from the Tabriz and Urmia stations are provided first. Then, the results of the reconstructed time series using the proposed MEMD-MFDFA framework are presented, followed by the results of the MFDFA analysis of the reconstructed data. In the final subsection, the results of the cross-correlation analysis between different variables and ET_o are presented.

4.1. MFDFA of ET_o and Meteorological Variables

In the MF-DFA framework, the presence of seasonality in a time series is associated with potential crossovers in the log–log plots of F_q(s) vs. s [28]. This characteristic can significantly impact the results and their interpretation, primarily because of the existence of different scaling regimes in the multifractal properties. An effective approach to removing seasonality is to use seasonal-trend decomposition based on Loess (STL) [56]. STL applies a filtering technique to decompose a time series into its seasonal, trend, and remainder components. In this study, the STL method is employed to remove the seasonal and trend components from the six hydro-meteorological time-series signals of the Tabriz and Urmia stations. The remainder component is then used as the input for performing the MF-DFA analysis. Thus, the MFDFA is invoked for investigating the multifractal properties of the different time series, considering the moment orders varying between −4 and 4 to avoid any possible bias [57]. The appropriate scaling range was chosen, which varies from 8 days to 900 days (~L/10, where L is the signal length), and was found to be depicting the scaling regions in most of the series, except wind velocity. The wind velocity series of both stations are very complex in the characteristics that possess scaling ranges from 8 days to 256 days. The fluctuation function plots of different variables are given in Figure 5 and Figure 6, respectively, for Tabriz and Urmia. Similarly, the GHE plots, mass exponent (Renyi exponent) plots, and multifractal spectra (singularity spectra) are provided in Figure 7. The mass exponent plots of the different variables of the two stations display good similarity.

From the GHE plots in Figure 7, it is noted that there exists a q-dependent relation between the scaling exponents of the different hydro-meteorological variables at both stations, by which we infer the multifractality of the time series. One can notice that the scaling exponents of the different variables of Urmia station are more than those of Tabriz station, even though the GHE plots of the temperatures and pressures of the stations are much closer. The Renyi exponent plots of all the variables for both stations are very close to each other. The multifractal spectra of the temperatures and pressures of both stations display good similarity. The spectra of the remaining variables of Urmia station display a rightward shift when compared to those of Tabriz station. The prominent multifractal properties such as Hurst exponent (H), spectral width (W), spread of GHE (∆h(q)), spread of singularity spectrum (∆f(α)), asymmetry index (R), and singularity exponent (α0) are computed for all-time series of both stations. The values of these properties are summarized in Table 2.

From Table 2, it is noted that all hydro-meteorological variables of both stations preserve long memory persistence (H > 0.5), with the lowest values for pressure and wind velocity, by which we infer the difficulties in the predictability of these series. The solar radiation time series shows the highest persistence and, as a result, the highest predictability, which is an expected outcome. Since the Hurst exponent measures the dependence of present values on past values, H is considered as an indicator of the predictability of the series. In other words, if a time series has persistence (as indicated by a Hurst exponent value greater than 0.5), it is also predictable [58]. Research conducted by various scholars has shown that the term “persistence” can be used as a criterion for assessing predictability, and the concepts of memory and predictability, as measures of persistence, are analogous [59]. Some of the researchers defined the predictability index [60] as $PI=2\left|D-1.5\right|$, where D is the fractal dimension defined as D = 2 − H [6]. The values of PI and the fractal dimension computed from H are given in Table 3. Table 3 indicates that the PI values of wind velocity and pressure are the lowest for both stations.

Few researchers have showcased the use of H as a predictability measure [58,59,60,61,62]. According to some of these researchers, for a time series to be considered predictable, not only the complete time series but also different sub-series should have H-values significantly higher than 0.5 [58]. In the presence of self-similar behavior, apart from the full-length time series, if the sub-series are also maintaining good persistence, they are considered to be predicable. To confirm this, an analysis was conducted by dividing all-time series into four sub-series with a length of L/4. This approach ensures that enough data points are preserved and helps to avoid any potential inconsistencies in the results during MFDFA application. Table 4 presents the results of the persistence of various sub-series for the different variables.

From Table 4, it is evident that both the complete time series and different sub-series of various variables exhibit persistent behavior, with none of the cases falling below 0.5. This clear indication of reasonably good predictability in the considered variables reinforces the findings of this study. This indicates that the series is predictable. For both stations, the spectral width, ∆h(q), ∆f(α), of the ET_o time series is largest when compared with the values of the respective properties of the remaining variables, by which we infer that the highest multifractality is for the ET_o time series when compared to the other variables. The values of the singularity exponent matches well with those of the Hurst exponent for all variables (giving an overall correlation of about 0.99 for both stations), supporting the association between the two dominant multifractal properties of the temperature datasets. The singularity spectra of all-time series (of both stations) are found to be left skewed (positive value of R), which indicates a high probability of high fluctuations.

4.2. Cause of Multifractality

Determining the reason for multifractality is one of the important steps in multifractal studies. Two of the major reasons behind multifractality are (i) long-range correlations and (ii) the fatness (spreadness) of the probability density function (PDF). In this work, the shuffling processes were followed to find the origins of the multifractality. The shuffling process destroys any long-range correlation in the series but will retain the same probability distribution. To quantify the effect of the broadness of the PDF, a surrogate (SU) series developed from the original series is used. If the multifractal character is due to a correlation property, the h(q) of the shuffled (SH) series will be reduced to 0.5. If it is due to the spreadness of the PDF, the h(q) obtained for the surrogate series will show q-independency [63]. If both the long-range temporal correlations and broadness of the PDFs are leading to multifractality, both the SH and SU series will show lower multifractal strength than the original series. The details of random shuffling and the steps for generating SU series by the iterative amplitude-adjusted Fourier transform (IAAFT) algorithm can be found in the literature [64]. The comparison of the GHE plots of the original, SH, and SU series will help to find the origins of the multifractal behavior. The GHE plots of the original, SH, and SU series of all datasets of the Tabriz and Urmia stations are presented in Figure 8.

Overall, from Figure 8, it is noted that, for all variables, there is an overall reduction in h(q) values of the SH and SU series, which implies that the multifractal behavior can be attributed to both the long-range correlations and fatness of the PDF. It is noticed that multifractality is not fully destroyed by the SU scheme. Also, the shuffling operation has reduced the GHE values to close to 0.5, exhibiting the dominant role of the temporal correlation in the multifractality of the series. The statistics of GHE values of the original, SR, and SU series, along with the ∆h(q) values, are presented in Table 5.Table 5 shows that the H values of the shuffled series have been brought down to a value close to 0.5 with a small SD (<0.03). Thus, it can be concluded that the multifractal behaviors of the different series are dominated by temporal correlations.

4.3. MEMD-MFDFA Approach for Fractality Detection

The proposed MEMD-MFDFA framework presented in Section 3.4 was invoked for time-series reconstruction for eliminating the nondominated (short memory modes) along with the Gaussian white noise expected in the high-frequency modes. The MEMD algorithm was applied to the multivariate dataset comprising six variables including evapotranspiration. The MEMD parameters of maximum threshold, minimum threshold, and fraction of decomposition were set as 0.75, 0.075, and 0.5, respectively, considering 64 direction vectors and following the guidelines given in past studies [49,50,51,65]. The decomposition resulted in 11 IMFs and one final residue for all six-time series. Subsequently, the MFDFA of each mode of each variable was performed, and the scaling exponent (value for q = 2) was extracted. The threshold value of the scaling exponent was fixed at 0.5, considering the classic categorization of long/short memory persistence of a time series [15]. The modes withvalues of scaling exponents falling below 0.5 were considered to be short memory and, therefore, were eliminated during the reconstruction process. The values of the scaling exponents of the rotatory modes of the different time series for the Tabriz and Urmia stations are presented in Figure 9 and Figure 10, respectively. For all the time series of Urmia station, the first five modes were found to be of short memory and were not considered for reconstruction. Among the different time series of Tabriz station, the temperature and ETo series for the first five modes were found to be of short memory, while for relative humidity and solar radiation, the first six modes were found to be of short memory. For the time series of pressure, the first four modes were found to be of short memory, while for the wind velocity time series of Tabriz station, modes 2 to 5 were found to be short memory. After eliminating the modes of short memory persistence, each time series was reconstructed and MFDFA was performed. The results of MFDFA for the reconstructed time series are presented in Table 6.

A comparison of the multifractal properties of the original series and reconstructed series (Table 2 and Table 4) show that the persistence and predictability of the different variables increased after de-noising and reconstruction.

4.4. MFCCA of Meteorological Variables with Reference Evapotranspiration

The MFCCA algorithm was then employed for finding the relationship between the different reconstructed meteorological variables with the ET_o of both stations. The persistence of individual series, joint persistence, correlations at different scales, and multifractal spectra were determined for each case. The multifractal spectra and correlation plots of different pairs of variables (such as T-ET_o, P-ET_o, RH-ET_o, U-ET_o, and SR-ET_o) for Tabriz station are presented in Figure 11, and the same for Urmia station are presented in Figure 12.

From Figure 11 and Figure 12, it is noticed that the shape of the multifractal spectra of different links are different for the two stations, while the spectra of the solar radiation and ET_o of Tabriz station shows good similarity. However, the patterns of correlations for different links display similarity, even though the magnitudes of correlation differ. At all-time scales, the association of T, U, and SR with ET_o are positive, while RH-ET_o correlations are negative for both stations. It may be noted that P-ET_o association for Tabriz station is negative at all time scales, while at some intra-annual scales, the association is weakly positive for Urmia station. Thus, MFCCA is successful in capturing such scale-dependent changes in the nature of the associations between different candidate hydro-climatic variables with ET_o.

The cross-correlation coefficient (indicated as ρ_CCA in this paper) corresponding to 90 days is termed a seasonal cross-correlation coefficient and that corresponding to 365 days is termed an annual cross-correlation coefficient. The seasonal, annual, and overall correlations, along with the scaling exponent, spectral width (W), and R for different variables with ET_o, are provided in Table 7.

From Table 7, it is noted that the combined scaling exponent (λxy) is nearly half of the scaling exponent (λx and λy) of the individual series. This is in agreement with the universal property on the persistence of time series by DCCA, proven by the researchers considering the sunspot–streamflow linkages. The seasonal and annual correlations with ET_o were found to be positive for all variables except RH and P. In all links, the seasonal scale correlations were weak, while the annual scale correlations were found to be strong at both stations. The strongest annual scale correlations and overall correlations were noted in the SR-ET_o link and T-ET_o link, which is a quite anticipated result on considering the climatology of northwestern Iran. High multifractality was noted in the ET_o series of both stations (spectral widths of 0.809 and 0.942 for Tabriz and Urmia, respectively). The multifractal spectra obtained for all combinations were right skewed, which was also indicated by the negative value of the asymmetric index. This implies that the extreme events play a prominent role in the temporal structure of the time series [27], indicating the high probability of low fluctuations. This study performed the multifractal analysis of daily reference ET_o time series along with relevant meteorological variables affecting the process of evapotranspiration. This study was performed for the two stations of Tabriz and Urmia in northwestern Iran. The climatic processes of the region are controlled by the elevation and presence of Lake Urmia, along with the large-scale global climatic oscillations [32]. Thus, the scale-dependent relation of the meteorological variables of the region with ETo are controlled by these factors. Lake–atmosphere interactions create a complex coupled system, in which feedback currents between them can shape and modify the local/global climate conditions on Earth. As a hypersaline lake, Lake Urmia has impacts on its local atmosphere (also on the global climate in the largescale), of which the most important ones can be mentioned: like other great lakes, Lake Urmia (a) serves as a natural modulator of temperature during the seasons, resulting in cooler summers and warmer winters, (b) contributes to cloud formation, which helps adjust the Earth’s radiation equilibrium and leads to a series of intricate feedback mechanisms in atmospheric circulation, ultimately increasing precipitation and snowfall in the region, (c) reduces, by increasing the cloud cover, the incoming atmospheric radiation to the Earth’s surface, and (d) plays a significant role in wind formation by establishing a temperature gradient (the area surrounding a lake is consistently warmer than the lake itself, as water reflects most of the solar radiation back into the atmosphere, while land surfaces absorb more atmospheric radiation). In addition (e), as the water level drops, it results in a decrease in the local relative humidity and the lake breezes, as well as an increase in the temperature of and evaporation from the surrounding terrain, leading to salty sandstorms. The presence of airborne dust, including in salt and dust storms, poses risks of salinity for agriculture, adverse weather conditions for human life, and atmospheric dimming. Dimmed weather reduces incoming radiation and acts in a manner comparable to the greenhouse effect [66,67,68,69,70,71,72].

Researchers in the region are conducting in-depth studies on hydroclimatic issues. One specific task in this context is to determine which province (left or right) contributes more to the problem from both natural and anthropogenic perspectives. The interactions among the various parameters may have played a role in the fluctuation and nonlinear dynamic characteristics of the climatic variables and ET_o of the region. This made differences and increased complexity, which is more evident in the series involving wind speed and surface pressure at Urmia station. These insights into the nonlinear dynamic patterns are considered valuable for enhancing the modeling of ET_o and, consequently, for improving irrigation-water management. Furthermore, the recent changes in hydro-climatic extremes may be affecting the data characteristics, and to study their impact on fractality, it is recommended that more recent and extended datasets be utilized. Employing datasets covering a longer time span can aid in capturing the impact of evolving climate conditions on the ET_o in any region, given that the methodology presented is applicable across different contexts.

The differences in multifractal properties may provide useful insights for developing prediction models of ET_o, which are expected to display improved accuracy. The proposed MEMD-MFDFA framework could be used as a potential tool for signal reconstruction, and the knowledge of short or long memory persistence can be useful for improving the predictability of the series. The mutual association of different reconstructed meteorological variables with evapotranspiration was examined in a multifractal perspective through the MFCCA, a recently proposed sign-preserved variant in the family of cross-correlation algorithms. However, more experiments need to be solicited using the methods used in the study by considering long and continuous hydro-meteorological datasets of different climatic conditions. The knowledge gained from the study is to be extended for the prediction of ET_o or related variables, integrating with the potential of hybrid machine learning or multifractal algorithms, in which, first, a decomposition-based data preprocessing can be undertaken for capturing scale-based associations. Subsequently, the scale-dependent modeling can be followed, and integration can be undertaken using machine learning methods or multifractal models [52]. Also, the presented MEMD framework has potential not only as a reconstruction tool but also opens the scope for performing scale-dependent multifractal modeling, and studies in these directions are in the pipeline.

5. Conclusions

This study investigated the multifractal characteristics of hydro-meteorological datasets of the Tabriz and Urmia stations of northwestern Iran using the MFDFA and MFCCA methods. MFDFA analysis showed that the daily reference evapotranspiration (ET_o) and five other meteorological variables (temperature, pressure, relative humidity_, wind velocity, and solar radiation) of both stations have multifractality and long memory persistence. The ET_o time series of both stations exhibited the highest multifractality, while pressure and wind velocity exhibited the lowest persistence. The proposed MEMD-MFDFA coupled framework is a promising approach for reconstructing a hydro-meteorological time series. MFDFA analysis of the reconstructed signals highlighted an increase in persistence and multifractality behaviors compared to the original time series. The knowledge on short or long memory persistence gathered from the proposed framework can be used for improving the predictability of the series. The application of the MFCCA method could capture the scale-dependent association of different meteorological parameters with ET_o. It was found that the joint persistence is nearly half of the persistence of the individual time series, which reinforced the universal property in the multifractal cross-correlation studies. The seasonal and annual cross-correlation properties displayed diversity in the strengths of different pair-wise combinations of cross-correlation analysis, but they displayed similar patterns for both stations. The knowledge gained from the study is helpful for developing decomposition-aided hybrid modeling frameworks for ET_o prediction, in which MEMD can be used as a decomposition paradigm, and multifractal models can be used for the predication of its components at different scales.

Author Contributions

Conceptualization, A.S. (Adarsh Sankaran) and T.P.; methodology, A.S. (Adarsh Sankaran) and T.P.; software, A.S. (Adarsh Sankaran); formal analysis, S.V. and A.S. (Aayisha Salim); investigation, S.V. and A.S. (Aayisha Salim); data curation, V.N.; writing—original draft preparation, A.S. (Adarsh Sankaran); writing—review and editing, A.S. (Adarsh Sankaran), T.P. and V.N.; visualization, S.V. and A.S. (Aayisha Salim); supervision, A.S. (Adarsh Sankaran) and T.P. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data cannot be shared in public platforms. The data used in the study will be transferred for academic purposes upon genuine requests to authors after receiving permission from the Iranian meteorological agencies.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Locations of Tabriz and Urmia in northwestern Iran.

Figure 2. Time-series plots of hydro-meteorological variables at Tabriz station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration (ETo).

Figure 3. Time-series plots of hydro-meteorological variables at Urmia station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration (ETo).

Figure 4. Flowchart of MEMD-MFDFA method for signal reconstruction.

Figure 5. Fluctuation function plots of different variables at Tabriz station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration.

Figure 6. Fluctuation function plots of different variables at Urmia station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration.

Figure 7. Multifractal characteristics of hydro-meteorological variables at Tabriz and Urmia stations for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration. The upper panels show the GHE plots, middle panels show the Renyi exponent plots, and the lower panels show the singularity spectrum.

Figure 8. Generalized Hurst exponentplots of original, SHand SU time series of different hydro-meteorological variables at Tabriz and Urmia stations for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration. The upper panels show the plots of Tabriz station, and the lower panels show the plots of Urmia station.

Figure 9. Plots of scaling exponents (for q = 2) of the modes of different variables at Tabriz station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration.

Figure 10. Plots of scaling exponents of the modes of different variables at Urmia station for (a) temperature, (b) pressure, (c) relative humidity, (d) wind velocity, (e) solar radiation, and (f) reference evapotranspiration.

Figure 11. Multifractal spectra and cross-correlation plots obtained by MFCCA of different links at Tabriz station for (a) temperature–ETo, (b) pressure–ETo, (c) relative humidity–ETo, (d) wind velocity–ETo, and (e) solar radiation–ETo.

Figure 12. Multifractal spectra and cross-correlation plots obtained by MFCCA of different links at Urmia station for (a) temperature–ETo, (b) pressure–ETo, (c) relative humidity–ETo, (d) wind velocity–ETo, and (e) solar radiation–ETo.

Table 1

Statistical properties of hydro-meteorological variables of the Tabriz and Urmia stations. SD and CV respectively show the standard deviation and the coefficient of variation.

Station	Variable	From	To	Statistical Property
				Maximum	Minimum	Mean	SD	CV (%)
Tabriz	Temperature (°C)	1/6/1992	26/10/2018	33.20	−15.000	13.34	10.38	77.81
	Pressure (kPa)	1/6/1992	26/10/2018	878.44	850.95	864.40	4.32	0.50
	Relative Humidity (%)	1/6/1992	26/10/2018	95.875	10.5	50.57	17.683	34.96
	Wind Velocity (ms−1)	1/6/1992	26/10/2018	9.500	0.00	3.384	1.60	47.28
	Solar Radiation (kJ m−2)	1/6/1992	26/10/2018	3479	43.00	1698.0	792.22	46.65
	Evapotranspiration (mm)	1/6/1992	26/10/2018	21.20	−0.185	6.1598	4.88935	79.37
Urmia	Temperature (°C)	1/6/1992	10/11/2018	30.00	−13.00	11.88	9.42	79.29
	Pressure (kPa)	1/6/1992	26/10/2018	882.00	853.24	867.19	4.47	0.515
	Relative Humidity (%)	1/6/1992	26/10/2018	99.50	20.13	57.99	15.69	27.06
	Wind Velocity (ms−1)	1/6/1992	26/10/2018	7.63	0.00	2.13	1.11	52.11
	Solar Radiation (kJ m−2)	1/6/1992	26/10/2018	3540.00	16.00	1774.16	805.30	45.39
	Evapotranspiration (mm)	1/6/1992	26/10/2018	13.80	0.00	4.13	3.39	82.08

Table 2

Multifractal properties of different hydro-meteorological variables of the Tabriz and Urmia stations.

Station	Variable	H	W	R	∆f(α)	∆h(q)	α0
Tabriz	Temperature (T)	0.818	0.298	0.296	0.192	0.149	0.845
	Pressure (P)	0.601	0.366	0.237	0.170	0.207	0.646
	Relative Humidity (RH)	0.754	0.213	0.495	0.239	0.107	0.771
	Wind Velocity (WV)	0.607	0.430	0.346	0.318	0.230	0.651
	Solar Radiation (SR)	0.845	0.438	0.445	0.340	0.244	0.883
	Reference Evapotranspiration (ET0)	0.746	0.735	0.487	0.600	0.430	0.810
Urmia	Temperature (T)	0.821	0.331	0.341	0.231	0.182	0.856
	Pressure (P)	0.619	0.374	0.157	0.097	0.213	0.666
	Relative Humidity (RH)	0.785	0.276	0.247	0.137	0.146	0.814
	Wind Velocity (WV)	0.659	0.274	0.146	0.075	0.166	0.700
	Solar Radiation (SR)	0.960	0.451	0.073	0.030	0.240	0.991
	Reference Evapotranspiration (ET0)	0.857	0.927	0.686	0.870	0.607	0.919

Table 3

Fractal dimension and predictability index for the different variables.

Station	Variable	D	PI
Tabriz	T	1.182	0.636
	P	1.399	0.202
	RH	1.246	0.508
	WV	1.393	0.214
	SR	1.155	0.690
	ET0	1.254	0.492
Urmia	T	1.179	0.642
	P	1.381	0.238
	RH	1.215	0.570
	WV	1.341	0.318
	SR	1.040	0.920
	ET0	1.143	0.714

Table 4

Persistence of the total time series and the different sub-series for all variables (CL: complete length; FH: first half; SH: second half; FQ: first quarter; SQ: second quarter; TQ: third quarter; FQ: fourth quarter; FRQ: first three quarter; LTQ: last three quarter).

Station	Variable	CL	FH	SH	MH	FQ	SQ	TQ	FRQ	FTQ	LTQ
Tabriz	T	0.818	0.800	0.867	0.828	0.856	0.782	0.904	0.876	0.813	0.829
	P	0.601	0.639	0.644	0.629	0.720	0.678	0.682	0.664	0.617	0.608
	RH	0.754	0.717	0.815	0.785	0.740	0.736	0.871	0.815	0.758	0.779
	WV	0.607	0.639	0.582	0.591	0.622	0.580	0.605	0.556	0.643	0.573
	SR	0.845	0.908	0.815	0.908	0.888	0.878	0.893	0.690	0.895	0.848
	ET0	0.746	0.763	0.768	0.783	0.743	0.741	0.797	0.793	0.760	0.769
Urmia	T	0.821	0.804	0.883	0.844	0.874	0.798	0.910	0.889	0.822	0.837
	P	0.619	0.643	0.661	0.638	0.732	0.650	0.702	0.683	0.630	0.628
	RH	0.785	0.786	0.826	0.802	0.797	0.759	0.843	0.837	0.797	0.808
	WV	0.659	0.624	0.722	0.610	0.681	0.534	0.617	0.750	0.679	0.691
	SR	0.960	0.994	0.967	0.991	0.922	0.999	0.990	0.877	0.996	0.977
	ET0	0.857	0.860	0.868	0.868	0.814	0.782	0.833	0.850	0.871	0.871

Table 5

Statistics of the generalized Hurst exponent of the original, SH, and SU series of the different datasets of the Tabriz and Urmia stations. SD is the standard deviation.

Variable	Property	Tabriz Station			Urmia Station
		Original Series	Shuffled Series	SurrogateSeries	OriginalSeries	ShuffledSeries	SurrogateSeries
Temperature	Mean	0.856	0.501	0.790	0.871	0.506	0.855
	SD	0.045	0.010	0.042	0.057	0.008	0.035
	∆h(q)	0.149	0.033	0.140	0.182	0.027	0.121
Pressure	Mean	0.659	0.486	0.631	0.678	0.511	0.632
	SD	0.065	0.008	0.033	0.067	0.012	0.048
	∆h(q)	0.207	0.028	0.105	0.213	0.042	0.157
Relative Humidity	Mean	0.782	0.515	0.753	0.823	0.507	0.786
	SD	0.033	0.008	0.031	0.045	0.003	0.016
	∆h(q)	0.107	0.028	0.100	0.146	0.008	0.056
Wind Velocity	Mean	0.649	0.535	0.616	1.021	0.520	0.999
	SD	0.059	0.010	0.018	0.073	0.011	0.012
	∆h(q)	0.184	0.033	0.058	0.240	0.037	0.041
Solar Radiation	Mean	0.912	0.506	0.856	1.033	0.507	0.907
	SD	0.078	0.016	0.003	0.210	0.013	0.003
	∆h(q)	0.244	0.055	0.000	0.607	0.047	0.010
ReferenceEvapotranspiration	Mean	0.864	0.517	0.764	0.708	0.494	0.716
	SD	0.138	0.009	0.013	0.054	0.017	0.010
	∆h(q)	0.430	0.033	0.047	0.166	0.057	0.034

Table 6

Multifractal properties of reconstructed hydro-meteorological variables of Tabriz and Urmia stations.

Station	Variable	H	W	R	∆f(α)	∆h(q)	α0
Tabriz	Temperature	0.875	0.666	0.522	0.537	0.407	0.933
	Pressure	0.662	0.584	0.345	0.340	0.343	0.726
	Relative Humidity	0.943	0.756	0.471	0.510	0.480	0.919
	Wind Velocity	0.884	0.322	0.196	0.102	0.168	0.914
	Solar Radiation	0.922	0.992	0.629	0.578	0.369	0.966
	Reference ET0	0.865	0.942	0.598	0.735	0.631	0.947
Urmia	Temperature	0.910	0.642	0.536	0.503	0.403	0.967
	Pressure	0.871	0.897	0.506	0.601	0.577	0.952
	Relative Humidity	0.947	0.811	0.485	0.544	0.518	0.927
	Wind Velocity	0.998	0.471	0.279	0.218	0.253	0.997
	Solar Radiation	0.987	0.825	0.522	0.529	0.519	0.948
	Reference ET0	0.909	0.871	0.556	0.642	0.582	0.993

Table 7

Results of MFCCA between different meteorological variables with ET_o: x refers to the meteorological variable; y refers to ET_o;λ refers to scaling exponent;ρ refers to cross-correlation.

Station	Link	Scaling Exponent			Cross-Correlation Coefficient			Spectral Width			Asymmetry Index
		λx	λy	λxy	ρs	ρa	ρo	Wx	Wy	Wxy	Rx	Ry	Rxy
Tabriz	T-ETo	0.990	0.959	0.975	0.655	0.973	0.908	0.304	0.809	0.502	−0.014	−0.584	−0.483
	P-ETo	0.693	0.959	0.826	−0.134	−0.715	−0.416	0.274	0.809	0.375	−0.055	−0.584	−0.448
	RH-ETo	0.838	0.959	0.899	−0.476	−0.935	−0.748	0.240	0.809	0.477	−0.055	−0.584	−0.402
	U-ETo	0.771	0.959	0.865	0.271	0.856	0.474	0.190	0.809	0.432	−0.273	−0.584	−0.505
	SR-ETo	0.991	0.959	0.980	0.349	0.930	0.758	0.831	0.809	0.486	−0.523	−0.584	−0.435
Urmia	T-ETo	0.905	0.980	0.943	0.431	0.947	0.848	0.380	0.942	0.452	−0.039	−0.651	−0.476
	P-ETo	0.706	0.980	0.843	−0.048	−0.688	−0.422	0.272	0.942	0.484	−0.084	−0.651	−0.501
	RH-ETo	0.858	0.980	0.919	−0.338	−0.889	−0.636	0.229	0.942	0.554	−0.112	−0.651	−0.470
	U-ETo	0.922	0.980	0.951	0.077	0.473	0.155	0.252	0.942	0.458	−0.052	−0.651	−0.574
	SR-ETo	0.989	0.980	0.985	0.274	0.928	0.739	0.667	0.942	0.582	−0.393	−0.651	−0.517

Appendix A

The schematic view of the MFDFA algorithm is provided in Figure A1.

The Hurst exponent (h(q = 2)) may decipher the persistence of the signal, while the deduced results, such as the Renyi exponent (τ(q)) and multifractal spectrum, may help to decode the multifractal behavior of the signal. The multifractal exponents can be computed as [17]:(A1)$$\tau (q)=qh(q)-1$$ (A2)$$\alpha =\frac{d\tau (q)}{dq}$$ and (A3)$$f(\alpha )=q\alpha -\tau (q)$$

The f(α) vs. α plot is called the multifractal spectrum and the support of the spectrum.

($\Delta \alpha ={\alpha }_{\max }-{\alpha }_{\min }$) shows the multifractal strength of the time series. The wider the spectrum, the higher will be the strength of the multifractality [18]. For a multifractal signal, the geometry of the plot will be a downward parabola that converges to a point for a monofractal series. The asymmetry index (R) is deduced from the singularity spectrum to indicate the fluctuation frequency:(A4)$$R=\frac{\Delta {\alpha }_L-\Delta {\alpha }_R}{\Delta {\alpha }_L+\Delta {\alpha }_R}$$ where $\Delta {\alpha }_R={\alpha }_{\max }-{\alpha }_0$ and $\Delta {\alpha }_L={\alpha }_0-{\alpha }_{\min }$ are the widths of the right and left limbs of the spectrum, respectively, and ${\alpha }_0$ is the singularity exponent for q = 0 (known as the Holder exponent), $R\in [-1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}1]$, and the positive value of this index represents the right-hand deviation of the multifractal spectrum, indicating a high probability of high fluctuations. The negative value of R shows that the series is characterized by localized low fluctuations.

Figure A1. Procedure of MFDFA.

References

1. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56. Fao Rome; 1998; 300, D05109.

2. Adarsh, S.; Nityanjaly, L.J.; Sarang, R.; Pham, Q.B.; Ali, M.; Nandhineekrishna, P. Multifractal Characterization and Cross correlations of Reference Evapotranspiration Time Series of India. Eur. Phys. J. Spec. Top.; 2021; 230, pp. 3845-3859. [DOI: https://dx.doi.org/10.1140/epjs/s11734-021-00325-4]

3. Goyal, P.; Kumar, S.; Sharda, R. A review of the Artificial Intelligence (AI) based techniques for estimating reference evapotranspiration: Current trends and future perspectives. Comput. Electron. Agric.; 2023; 209, 107836. [DOI: https://dx.doi.org/10.1016/j.compag.2023.107836]

4. Hurst, H.E. Long-term storage capacity of reservoirs. Trans. Am. Soc. Civ. Eng.; 1951; 116, pp. 770-808. [DOI: https://dx.doi.org/10.1061/TACEAT.0006518]

5. Mandelbrot, B. The Fractal Geometry of Nature; WH Freeman Publishers: New York, NY, USA, 1982.

6. Pandey, G.; Lovejoy, S.; Schertzer, D. Multifractal analysis of daily river flows including extremes for basins five to two million square kilometers, one day to 75 years. J. Hydrol.; 1998; 208, pp. 62-81. [DOI: https://dx.doi.org/10.1016/S0022-1694(98)00148-6]

7. Kantelhardt, J.W.; Rybski, D.; Zschiegner, S.A.; Braun, P.; Koscielny-Bunde, E.; Livina, V.; Havlin, S.; Bunde, A. Multifractality of river runoff and precipitation: Comparison of fluctuation analysis and wavelet methods. Phys. A Stat. Mech. Its Appl.; 2003; 330, pp. 240-245. [DOI: https://dx.doi.org/10.1016/j.physa.2003.08.019]

8. Huang, Y.; Schmitt, F.G.; Lu, Z.; Liu, Y. Analysis of daily river flow fluctuations using empirical mode decomposition and arbitrary order Hilbert spectral analysis. J. Hydrol.; 2009; 373, pp. 103-111. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2009.04.015]

9. Schmitt, F.; Schertzer, D.; Lovejoy, S. Multifractal analysis of foreign exchange data. Appl. Stoch. Models Data Anal.; 1999; 15, pp. 29-53. [DOI: https://dx.doi.org/10.1002/(SICI)1099-0747(199903)15:1<29::AID-ASM357>3.0.CO;2-Z]

10. Plocoste, T.; Calif, R.; Jacoby-Koaly, S. Temporal multiscaling characteristics of particulate matter PM10 and ground-level ozone O₃ concentrations in Caribbean region. Atmos. Environ.; 2017; 169, pp. 22-35. [DOI: https://dx.doi.org/10.1016/j.atmosenv.2017.08.068]

11. Carmona-Cabezas, R.; Ariza-Villaverde, A.B.; de Ravé, E.G.; Jiménez-Hornero, F.J. Visibility graphs of ground-level ozone time series: A multifractal analysis. Sci. Total Environ.; 2021; 661, pp. 138-147. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2019.01.147]

12. Plocoste, T.; Carmona-Cabezas, R.; Jiménez-Hornero, F.J.; de Ravé, E.G.; Calif, R. Multifractal characterisation of particulate matter (PM10) time series in the Caribbean basin using visibility graphs. Atmos. Pollut. Res.; 2021; 12, pp. 100-110. [DOI: https://dx.doi.org/10.1016/j.apr.2020.08.027]

13. Peng, C.K.; Buldyrev, S.V.; Havlin, S.; Simons, M.; Stanley, H.E.; Goldberger, A.L. Mosaic organization of DNA nucleotides. Phys. Rev. E; 1994; 49, 1685. [DOI: https://dx.doi.org/10.1103/PhysRevE.49.1685]

14. Kantelhardt, J.W.; Zschiegner, S.A.; Koscielny-Bunde, E.; Havlin, S.; Bunde, A.; Stanley, H.E. Multifractal detrended fluctuation analysis of non-stationary time series. Phys. A Stat. Mech. Its Appl.; 2002; 316, pp. 87-114. [DOI: https://dx.doi.org/10.1016/S0378-4371(02)01383-3]

15. Kantelhardt, J.W.; Koscielny-Bunde, E.; Rybski, D.; Braun, P.; Bunde, A.; Havlin, S. Long-term persistence and multifractality of precipitation and river runoff records. J. Geophys. Res. Atmos.; 2006; 111, [DOI: https://dx.doi.org/10.1029/2005JD005881]

16. Koscielny-Bunde, E.; Kantelhart, J.W.; Braun, P.; Bunde, A.; Havlin, S. Long-term persistence and multifractality of river runoff records: Detrended fluctuation studies. J. Hydrol.; 2006; 322, pp. 120-137. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2005.03.004]

17. Li, E.; Mu, X.; Zhao, G.; Gao, P. Multifractal detrended fluctuation analysis of streamflow in the Yellow River Basin, China. Water; 2015; 7, pp. 1670-1686. [DOI: https://dx.doi.org/10.3390/w7041670]

18. Hou, W.; Feng, G.; Yan, P.; Li, S. Multifractal analysis of the drought area in seven large regions of China from 1961 to 2012. Meteorol. Atmos. Phys.; 2018; 130, pp. 459-471. [DOI: https://dx.doi.org/10.1007/s00703-017-0530-0]

19. Adarsh, S.; Dharan, D.S.; Nandhu, A.R.; Anand Vishnu, B.; Mohan, V.K.; Watorek, M. Multifractal Description of Streamflow and Suspended Sediment Concentration Data from Indian River Basins. Acta Geophys.; 2020; 68, pp. 519-535. [DOI: https://dx.doi.org/10.1007/s11600-020-00407-2]

20. Adarsh, S.; Krzyszczak, J.; Baranowski, P.; Archana, D.S.; Nandhinee krishna, P.; Nityanjaly, L.J.; Vandana, T.; Ali, M. Multifractal Cross Correlation Analysis of Agro-Meteorological Datasets (Including Reference Evapotranspiration) of California. Atmosphere; 2020; 11, 1116. [DOI: https://dx.doi.org/10.3390/atmos11101116]

21. Martínez, J.L.M.; Segovia-Domínguez, I.; Rodríguez, I.Q.; Horta-Rangel, F.A.; Sosa-Gómez, G. A modified Multifractal Detrended Fluctuation Analysis (MFDFA) approach for multifractal analysis of precipitation. Phys. A Stat. Mech. Its Appl.; 2021; 565, 125611. [DOI: https://dx.doi.org/10.1016/j.physa.2020.125611]

22. Zhan, C.; Liang, C.; Zhao, L.; Jiang, S.; Niu, K.; Zhang, Y. Multifractal characteristics of multiscale drought in the Yellow River Basin, China. Phys. A Stat. Mech. Its Appl.; 2023; 609, 128305. [DOI: https://dx.doi.org/10.1016/j.physa.2022.128305]

23. Lin, G.; Chen, X.; Fu, Z. Temporal–spatial diversities of long-range correlation for relative humidity over China. Phys. A Stat. Mech. Its Appl.; 2007; 383, pp. 585-594. [DOI: https://dx.doi.org/10.1016/j.physa.2007.04.059]

24. Yuan, N.; Fu, Z.; Mao, J. Different multi-fractal behaviors of diurnal temperature range over the north and the south of China. Theor. Appl. Climatol.; 2013; 112, pp. 673-682. [DOI: https://dx.doi.org/10.1007/s00704-012-0762-3]

25. Burgueño, A.; Lana, X.; Serra, C.; Martínez, M.D. Daily extreme temperature multifractals in Catalonia (NESpain). Phys. Lett. A; 2014; 378, pp. 874-885. [DOI: https://dx.doi.org/10.1016/j.physleta.2014.01.033]

26. Kalamaras, N.; Philippopoulos, K.; Deligiorgi, D.; Tzanis, C.G.; Karvounis, G. Multifractal scaling properties of daily air temperature time series. Chaos Solitons Fractals; 2017; 98, pp. 38-43. [DOI: https://dx.doi.org/10.1016/j.chaos.2017.03.003]

27. Krzyszczak, J.; Baranowski, P.; Zubik, M.; Kazandjiev, V.; Georgieva, V.; Sławiński, C.; Siwek, K.; Kozyra, J.; Nieróbca, A. Multifractal characterization and comparison of meteorological time series from two climatic zones. Theor. Appl. Climatol.; 2019; 137, pp. 1811-1824. [DOI: https://dx.doi.org/10.1007/s00704-018-2705-0]

28. Plocoste, T.; Pavón-Domínguez, P. Multifractal detrended cross-correlation analysis of wind speed and solar radiation. Chaos Interdiscip. J. Nonlinear Sci.; 2020; 30, 113109. [DOI: https://dx.doi.org/10.1063/5.0026354]

29. Baranowski, P.; Krzyszczak, J.; Slawinski, C.; Hoffmann, H.; Kozyra, J.; Nieróbca, A.; Siwek, K.; Gluza, A. Multifractal analysis of meteorological time series to assess climate impacts. Clim. Res.; 2015; 65, pp. 39-52. [DOI: https://dx.doi.org/10.3354/cr01321]

30. Nourani, V.; Elkiran, G.; Abdullah, J. Multi-station artificial intelligence-based ensemble modeling of reference evapotranspiration using pan evaporation measurements. J. Hydrol.; 2019; 577, 123958. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2019.123958]

31. Nourani, V.; Elkiran, G.; Abdullah, J. Multi-step ahead modeling of reference evapotranspiration using a multi-model approach. J. Hydrol.; 2020; 581, 124434. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2019.124434]

32. Sreedevi, V.; Adarsh, S.; Nourani, V. Multiscale coherence analysis of reference evapotranspiration of north-western Iran using wavelet transform. J. Water Clim. Chang.; 2022; 13, pp. 505-521. [DOI: https://dx.doi.org/10.2166/wcc.2021.379]

33. Suman, A.; Archana, D.S.; Nayak, A.K.; Adarsh, S.; Biswal, B. Unveiling the Climatic origin of Streamflow persistence through Multifractal Analysis of Hydrometeorological datasets of India. Hydrol. Sci. J.; 2023; 68, pp. 290-306. [DOI: https://dx.doi.org/10.1080/02626667.2022.2157726]

34. Ariza-Villaverde, A.B.; Pavón-Domínguez, P.; Carmona-Cabezas, R.; Gutierrez de Rave, E.; Jimenez-Honero, F.J. Joint multifractal analysis of air temperature, relative humidity and reference evapotranspiration in the middle zone of the Guadalquivir river valley. Agric. For. Meteorol.; 2019; 278, 107657. [DOI: https://dx.doi.org/10.1016/j.agrformet.2019.107657]

35. Cun, Z.; Liang, C.; Zhao, L.; Zhang, Y.; Cheng, L.; Jiang, S.; Xing, L. Multifractal characteristics analysis of daily reference evapotranspiration in different climate zones of China. Phys. A Stat. Mech. Its Appl.; 2021; 583, 126273. [DOI: https://dx.doi.org/10.1016/j.physa.2021.126273]

36. Gómez-Gómez, J.; Ariza-Villaverde, A.B.; Gutiérrez de Ravé, E.; Jiménez-Hornero, F.J. Relationships between Reference Evapotranspiration and Meteorological Variables in the Middle Zone of the Guadalquivir River Valley Explained by Multifractal Detrended Cross-Correlation Analysis. Fractal Fract.; 2023; 7, 54. [DOI: https://dx.doi.org/10.3390/fractalfract7010054]

37. Nourani, V.; Fard, M.S. Sensitivity analysis of the artificial neural network outputs in simulation of the evaporation process at different climatologic regimes. Adv. Eng. Softw.; 2016; 47, pp. 127-146. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2011.12.014]

38. Podobnik, B.; Stanley, H.E. Detrended cross-correlation analysis: A new method for analyzing two nonstationary time series. Phys. Rev. Lett.; 2008; 100, 084102. [DOI: https://dx.doi.org/10.1103/PhysRevLett.100.084102]

39. Zhou, W.-X. Multifractal detrended cross-correlation analysis for two non-stationary signals. Phys. Rev. E; 2008; 77, 066211. [DOI: https://dx.doi.org/10.1103/PhysRevE.77.066211]

40. Jiang, Z.-Q.; Zhou, W.-X. Multifractal detrending moving-average cross-correlation analysis. Phys. Rev. E; 2011; 84, 016106. [DOI: https://dx.doi.org/10.1103/PhysRevE.84.016106]

41. Oświȩcimka, P.; Drożdż, S.; Forczek, M.; Jadach, S.; Kwapień, J. Detrended cross-correlation analysis consistently extended to multifractality. Phys. Rev. E; 2014; 89, 023305. [DOI: https://dx.doi.org/10.1103/PhysRevE.89.023305] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25353603]

42. Hajian, S.; Movahed, M.S. Multifractal detrended cross-correlation analysis of sunspot numbers and river flow fluctuations. Phys. A Stat. Mech. Its Appl.; 2010; 389, pp. 4942-4954. [DOI: https://dx.doi.org/10.1016/j.physa.2010.06.025]

43. Dey, P.; Mujumdar, P.P. Multiscale evolution of persistence of rainfall and streamflow. Adv. Water Resour.; 2018; 121, pp. 285-303. [DOI: https://dx.doi.org/10.1016/j.advwatres.2018.08.018]

44. Maghyereh, A.; Abdoh, H.; Wątorek, M. The impact of COVID-19 pandemic on the dynamic correlations between gold and U.S. equities: Evidence from multifractal cross-correlation analysis. Qual. Quant.; 2023; 57, pp. 1889-1903. [DOI: https://dx.doi.org/10.1007/s11135-022-01404-x]

45. Kwapień, J.; Oświęcimka, P.; Drożdż, S. Detrended fluctuation analysis made flexible to detect range of cross-correlated fluctuations. Phys. Rev. E; 2015; 92, 052815. [DOI: https://dx.doi.org/10.1103/PhysRevE.92.052815]

46. Brito, A.A.; Santos, F.R.; de Castro, A.P.N.; da Cunha Lima, A.T.; Zebende, G.F.; da Cunha Lima, I.C. Cross-correlation in a turbulent flow: Analysis of the velocity field using the σ DCCA coefficient. Eur. Lett.; 2018; 123, 20011. [DOI: https://dx.doi.org/10.1209/0295-5075/123/20011]

47. Wątorek, M.; Drożdż, S.; Kwapień, J.; Minati, L.; Oświȩcimka, P.; Stanuszek, M. Multiscale characteristics of the emerging global cryptocurrency market. Phys. Rep.; 2021; 901, pp. 1-82. [DOI: https://dx.doi.org/10.1016/j.physrep.2020.10.005]

48. Wątorek, M.; Drożdż, S.; Oświęcimka, P.; Stanuszek, M. Multifractal cross-correlations between the world oil and other financial markets in 2012–2017. Energy Econ.; 2019; 81, pp. 874-885. [DOI: https://dx.doi.org/10.1016/j.eneco.2019.05.015]

49. Rehman, N.; Mandic, D.P. Multivariate empirical mode decomposition. Proc. R. Soc. A; 2011; 466, pp. 1291-1302. [DOI: https://dx.doi.org/10.1098/rspa.2009.0502]

50. Hu, W.; Si, B.-C. Soil water prediction based on its scale-specific control using multivariate empirical mode decomposition. Geoderma; 2013; 193–194, pp. 180-188. [DOI: https://dx.doi.org/10.1016/j.geoderma.2012.10.021]

51. Huang, G.; Su, Y.; Kareem, A.; Liao, H. Time-frequency analysis of non-stationary process based on multivariate empirical mode decomposition. J. Eng. Mech.; 2016; 142, 04015065. [DOI: https://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0000975]

52. Adarsh, S.; Sulaiman, S.; Murshida, K.K.; Nooramol, P. Scale-dependent prediction of reference evapotranspiration based on Multivariate Empirical Mode Decomposition. Ain Shams Eng. J.; 2018; 9, pp. 1839-1848. [DOI: https://dx.doi.org/10.1016/j.asej.2016.10.014]

53. Adarsh, S.; Janga Reddy, M. Multiscale characterization and prediction of reservoir inflows using MEMD-SLR Coupled Approach. J. Hydrol. Eng.; 2019; 24, 04018059. [DOI: https://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0001732]

54. Ali, M.; Deo, R.C.; Maraseni, T.; Downs, N.J. Improving SPI-derived drought forecasts incorporating synoptic-scale climate indices in multi-phase multivariate empirical mode decomposition model hybridized with simulated annealing and kernel ridge regression algorithms. J. Hydrol.; 2019; 576, pp. 164-184. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2019.06.032]

55. Agana, N.A.; Homaifar, A. EMD-based predictive deep belief network for time series prediction: An application to drought forecasting. Hydrology; 2018; 5, 18. [DOI: https://dx.doi.org/10.3390/hydrology5010018]

56. Cleveland, R.B.; Cleveland, W.S. STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Off. Stat.; 1990; 6, pp. 3-33.

57. Drożdż, S.; Minati, L.; Oświȩcimka, P.; Stanuszek, M.; Wątorek, M. Signatures of the Crypto-Currency Market Decoupling from the Forex. Future Internet; 2019; 11, 154. [DOI: https://dx.doi.org/10.3390/fi11070154]

58. Chandrasekaran, C.; Poomalai, P.; Saminathan, B.; Suthanthiravel, S.; Sundaram, K.; Hakkim, F.F.A. An investigation on the relationship between the Hurst exponent and the predictability of a rainfall time series. Meteorol. Appl.; 2019; 16, pp. 511-519. [DOI: https://dx.doi.org/10.1002/met.1784]

59. Orth, R.; Seneviratne, S.I. Predictability of soil moisture and streamflow on subseasonal timescales: A case study. J. Geophys. Res. Atmos.; 2013; 118, pp. 963-979. [DOI: https://dx.doi.org/10.1002/jgrd.50846]

60. Rangarajan, G.; Sant, D.A. A climate predictability index and its applications. Geophys. Res. Lett.; 1997; 24, pp. 1239-1242. [DOI: https://dx.doi.org/10.1029/97GL01058]

61. Peyghami, M.R.; Khanduzi, R. Predictability and forecasting automotive price based on a hybrid train algorithm of MLP neural network. Neural Comput. Appl.; 2012; 21, pp. 125-132. [DOI: https://dx.doi.org/10.1007/s00521-011-0605-7]

62. Khalili, N.; Khodashenas, S.R.; Davary, K.; Baygi, M.M.; Karimaldini, F. Prediction of rainfall using artificial neural networks for synoptic station of Mashhad: A case study. Arab. J. Geosc.; 2016; 9, 624. [DOI: https://dx.doi.org/10.1007/s12517-016-2633-1]

63. Movahed, M.S.; Jafari, G.R.; Ghasemi, F.; Rahvar, S.; Tabar, M.R.R. Multifractal detrended fluctuation analysis of sunspot time series. J. Stat. Mech. Theory Exp.; 2006; 2006, P02003. [DOI: https://dx.doi.org/10.1088/1742-5468/2006/02/P02003]

64. Keylock, C.J. Hypothesis testing for nonlinear phenomena in the geosciences using synthetic, surrogate data. Earth Space Sci.; 2019; 6, pp. 41-58. [DOI: https://dx.doi.org/10.1029/2018EA000435]

65. Rilling, G.; Flandrin, P.; Goncalves, P. On empirical mode decomposition and its algorithms. Proceedings of the IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing NSIP-03; Grado, Italy, 8–11 June 2003; Volume 3, pp. 1-5.

66. Gerken, T.; Biermann, T.; Babel, W.; Herzog, M.; Ma, Y.; Foken, T.; Graf, H.F. A modelling investigation into lake-breeze development and convection triggering in the Nam Co Lake basin, Tibetan Plateau. Theor. Appl. Climatol.; 2014; 117, pp. 149-167. [DOI: https://dx.doi.org/10.1007/s00704-013-0987-9]

67. Sun, X.; Xie, L.; Semazzi, F.H.; Liu, B. A numerical investigationof the precipitation over Lake Victoria basin using a coupledatmosphere-lake limited-area model. Adv. Meteorol.; 2014; 2014, 960924. [DOI: https://dx.doi.org/10.1155/2014/960924]

68. Dehghanipour, A.H.; Moshir Panahi, D.; Mousavi, H.; Kalantari, Z.; Tajrishy, M. Effects of Water Level Decline in Lake Urmia, Iran, on Local Climate Conditions. Water; 2020; 12, 2153. [DOI: https://dx.doi.org/10.3390/w12082153]

69. Nourani, V.; Tootoonchi, R.; Andaryani, S. Investigation of climate, land cover and lake level pattern changes and interactions using remotely sensed data and wavelet analysis. Ecol. Inform.; 2021; 64, 101330. [DOI: https://dx.doi.org/10.1016/j.ecoinf.2021.101330]

70. Nourani, V.; Sayyah-Fard, M.; Kantoush, S.A.; Bharambe, K.P.; Sumi, T.; Saber, M. Optimization-based prediction uncertainty qualification of climatic parameters. J. Hydrometeorol.; 2023; [DOI: https://dx.doi.org/10.1175/JHM-D-23-0043.1]

71. Plocoste, T.; Calif, R.; Euphrasie-Clotilde, L.; Brute, F.N. Investigation of local correlations between particulate matter (PM10) and air temperature in the Caribbean basin using Ensemble Empirical Mode Decomposition. Atmos. Pollut. Res.; 2020; 11, pp. 1692-1704. [DOI: https://dx.doi.org/10.1016/j.apr.2020.06.031]

72. Plocoste, T.; Calif, R. Is there a causal relationship between Particulate Matter (PM10) and air Temperature data? An analysis based on the Liang–Kleeman information transfer theory. Atmos. Pollut. Res.; 2021; 12, 101177. [DOI: https://dx.doi.org/10.1016/j.apr.2021.101177]