key: cord-0078418-tvxa8798
authors: Millán, Humberto; Macías, Idalberto; Valderá, Nathalí
title: Multifractality of the standardized precipitation index: influence of pan evaporation and virtual temperature-based potential evapotranspiration
date: 2022-05-22
journal: Meteorol Atmos Phys
DOI: 10.1007/s00703-022-00894-6
sha: cc6c101cc24446e0a7cc1403ea9e49a185f85e4d
doc_id: 78418
cord_uid: tvxa8798

Hydrological extremes experience an increase in some regions and a decrease in other zones. The objectives of the present work were (i) to introduce Class A pan evaporation data and virtual temperature-based potential evapotranspiration (PET(v) hereafter) into the Standardized Precipitation–Evapotranspiration Index (SPEIm hereafter) computation and (ii) to describe small and large fluctuations of SPI and SPEIm through multifractal detrended fluctuation (MF-DFA) and multifractal detrended cross-correlation (MF-DCCA) analyses. We used 40 years data (1974–2013) of monthly rainfall (P), mean, minimum and maximum air temperature, pan evaporation (E), relative humidity (RH) and relative sun brightness (RSB). Meteorological variables were collected from Puyo meteorological station, Pastaza Province, Ecuador. SPI time series for 1 and 6 months timescales were determined following two approaches. We computed SPI values using precipitation as the only input variable. Additionally, we incorporated pan evaporation and virtual temperature-based potential evapotranspiration into the standard SPEI computation (SPEIm). The SPEIm revealed some differences as compared with the classical SPI methodology. Five out of fifteen Asymmetry Index (AI) values were positive (0.095 ≤ AI ≤ 0.419).This indicates the relevance of high fluctuations at different time scales. Joint multifractal spectra between SPI (1,6)/SPEIm(1,6) versus RH and RSB rendered negative AI values which suggests the importance of low fluctuations at shorter time scales. The DCCA cross-correlation coefficient allows one to identify those time scales where SPI and SPEIm are influenced by other meteorological variables. Long-term correlation and sub-Gaussian behaviour of meteorological variables (apart from air temperature) are the main causes of multifractal structures.

Global warming is exerting a major impact on developing countries with agricultural, economical, and societal vulnerabilities. In particular Ecuador still needs more drought investigations including its three geographical regions (Campozano et al. 2020) . Climate is a nonlinear complex system that is out of thermodynamic equilibrium (Ghil and Lucarini 2020) where many scales of turbulence can occur (Ghil 2019) . Climate-related extreme events (e.g. floods, droughts, hurricanes, heat/cold waves, etc.) are effects of equilibrium deviation of the climatic system (Kantz et al. 2004 ). The World Meteorological Organization has recognized the economic and humanitarian impact of droughts in many vulnerable parts around the world since 2015 (World Meteorological Organization 2019). The impacts of climate change in general and droughts in particular on agriculture and global food security are well documented in the literature (World Meteorological Organization 2003; Croppenstedt et al. 2018; Zhang et al. 2020; Sharifi 2021) .

The SPI (McKee et al. 1993) has been proposed by the World Meteorological Organization for describing the intensity and frequency of drought events (Hayes et al. 2011 ). However, some limitations of the SPI have been recognized. For example Mishra and Singh (2010) have identified the length of the precipitation record and the nature of the Probability Distribution Function as limiting factors. The World Meteorological Organization (2012) identified the use of only precipitation as the main weaknesses of the SPI. The SPI computation assumes that precipitation and other meteorological variables are stationary signals without significant trends (Vicente-Serrano et al. 2010) . That statement has been debatable as some variables (e.g. temperature, relative humidity and evaporation) can influence moisture dynamics within the present global warming scenario. For example, Abramopoulos et al. (1988) and Rebetez et al. (2006) have shown that an increase in temperature also increases drought characteristics as duration and severity. According to their specific characteristics, drought events have been classified into meteorological, agricultural, hydrological and socio-economical droughts. In general, some drought characteristics as intensity and duration can be computed from the consecutive number of months with SPI below a threshold value. Tong et al. (2018) have pointed out that most drought types begin with a meteorological drought. On the other hand the air temperature-based PET has been taken into account with the implementation of the Standardized Precipitation Evapotranspiration Index (SPEI hereafter) (Vicente-Serrano et al. 2010; Beguería et al. 2014) . Ma et al. (2013) also found some potential limitations of the SPEI methodology regarding the use of precipitation and temperature as input variables. Those authors incorporated soil water capacity to the SPEI computation. They developed a new standardized water balance-derived index (SWBI). According to the authors, the SWBI seems to be more robust than the SPEI. In an effort for improving drought indexes, some authors have proposed other modifications to the SPEI or SPI methodologies. For example, Beguería et al. (2014) suggested that actual evapotranspiration (ET a ) would be a better estimator of water deficit than precipitation. On the other hand the use of virtual temperature instead of air temperature with drought indicators could be also an interesting avenue to explore in the future. Virtual temperature depends on air humidity which in turn impacts convection. Valuable theoretical and practical information on the topic can be found in Monteiro and Torlaschi (2007) . Therefore, virtual temperature could be an important variable for tropical (humid) environments where relative humidity is over 85% along the year. The effect of virtual temperature correction on small Convective Available Potential Energy (CAPE) has been investigated by Doswell and Rasmussen (1994) . Those authors recommended that CAPE computations should include virtual temperature correction.

MF-DFA (Kantelhardt et al. 2002) and MF-DCCA (Zhou 2008) could provide valuable tools for analyzing e interpreting the relationship between the SPI and different meteorological signals in terms of their multifractal spectra. Regarding drought indicators, Zhang et al (2010) applied MF-DFA for investigating flood/drought grade series in the Yangtze Delta, China, at a millennium time scale. Those authors called the attention on the probability of predicting future climatic changes based on historical records in the Yangtze Delta region, China. Luo et al. (2014) found stationarity and long-range correlation of drought time series in China using MF-DFA. Hou et al. (2018) investigated, through MF-DFA, both persistence and anti-persistence of SPI signals in different regions of China. Those authors used 52 years (624 months) data points. Adarsh et al. (2019) also used MF-DFA for investigating SPI time series at 3-months, 6-months and 12-months timescales. Those authors partitioned each SPI signal into four time periods (1871-1905, 1906-1940, 1941-1975 and 1976-2016) . They found that the degree of multifractality increased with the SPI timescale especially after 1976. Ogunjo (2021) applied MF-DFA to a 30 years sample (360 months) for investigating multifractal properties of meteorological drought over different regions across Nigeria. That author found that the multifractal spectrum strength increased with the timescale of SPI computation. Even though MF-DCCA is a relatively recent promising tool for investigating meteorological extremes, its application with SPI or SPEI time series is still scarce.

The literature reports other studies using MF-DFA and MF-DCCA with some meteorological signals. For example Baranowski et al. (2015) used MF-DFA with several meteorological variables from Finland, Germany, Poland and Spain. From the multifractal analysis, those authors concluded that precipitation was the most sensitive variable to climate dynamics. Philippopoulos et al. (2019) applied MF-DFA for investigating scaling and multifractal characteristics of mean daily temperature time series of the ERA-Interim reanalysis data for Greece. Those authors found positive long-term correlation with multifractal structure insensitive to larger local fluctuation. Tropical regions have been also investigated with MF-DFA. Agbazo et al. (2019) also applied MF-DFA to daily temperature time series collected over Benin (West Africa). One of their main conclusions was that geographical position influences the multifractality of temperature. MF-DCCA (Zhou 2008) represented an improvement to multifractal statistics as it allows one to investigate the mutual correlation between two time series based on their multifractal parameters. Recently, Plocoste and Pavón-Domínguez (2020) applied MF-DCCA for investigating the influence of solar radiation intensity on wind speed. Those authors found that the degree of multifractal cross-correlation was identical for each site. Sankaran et al. (2020) found that MF-DCCA was sensitive to the scale dependent relationship of reference evapotranspiration, wind speed, incoming solar radiation, air temperature, air pressure and relative humidity and their associations from weekly to inter-annual time scales. In addition, meteorological variables as related to human health have been also investigated with MF-DCCA. For example, Wang et al. (2020) studied the multifractal cross-correlations between meteorological factors and diseases of bacterial origin. That investigation included 113 samples collected between the years 2010 and 2019.

Thus, we were motivated by two basic questions. How do class A pan evaporation data and calculated virtual temperature influence the information derived from the classical SPI signal? Second, are SPI and/or SPEIm sensitive to different atmospheric variables in terms of their corresponding local regularities? The objectives of the present work were (i) to introduce Class A pan evaporation data and virtual temperature-based potential evapotranspiration (PET) into the Standardized Precipitation-Evapotranspiration Index (SPEIm) computation and (ii) to describe small and large fluctuations of SPI and its modified counterpart through multifractal detrended fluctuation (MF-DFA) and multifractal detrended cross-correlation (MF-DCCA) analyses.

Data sets corresponding to the present study were collected from Puyo meteorological station (01° 30' South Latitude 77° 57' West Longitude, 950 m above sea level) in Pastaza Province, Ecuador. Figure 1 depicts the topographical map of the studied area separating Highlands (central area), Coast region (towards the Pacific Ocean), and Amazon region to the right. We collected 40 years data of monthly rainfall, mean, maximum and minimum air temperature, Class A pan evaporation, RH and RSB (percentage of time of sun brightness duration/12). Consequently the study used 480 months with no missing data points within the investigated time period. The climate type of the region is classified as Af according to the Köppen-Geiger climate classification (Peel et al. 2007 ). The study site limits with two contrasting environments. The warm Amazonian tropical forest located at the Northeast, East and Southeast and a volcanic belt situated at the Northwest, West and Southwest parts of the study region. In particular, there are several snowe-covered volcanoes around the study site all over 5000 m above sea level which generate a cold flux air from West to East (Millán et al. 2009 ).

Virtual temperature (T v ) is understood as the temperature of an air dry parcel with the same pressure and density as an equivalent moist air parcel (Jacobson 2005) . We calculated the virtual temperature from air temperature (T), atmospheric pressure (P), relative humidity (RH) and actual vapour pressure (e a ) (Allen et al. 1998 ): where T v and T are expressed in Kelvin (K) while e a and P are given in kPa.

Actual vapour pressure (e a ) was computed from monthly air temperature (T) and monthly relative humidity (RH) (Allen et al. 1998) :

where e 0 (T) is the saturation vapour pressure (kPa) at air temperature ( 0 C).

In this case the saturation vap**our pressure (e 0 (T)) was computed as (Allen et al. 1998 ) follows:

The atmospheric pressure (P) at the height above sea level (z-z 0 ) was calculated according to Burman et al. (1987) as follows:

where P 0 = 101.325 kPa (atmospheric pressure at the sea level), z-z 0 = 950 m is the height above sea level of the studied site (z 0 = 0 at the sea level), g = 9.81 ms −2 is the gravitational acceleration, R = 287 J kg −1 K −1 is the specific gas constant, δ = 6.5 × 10 -3 K m −1 is the constant lapse rate moist air and T i [K] is the reference air temperature corresponding to the i-month.

The mathematical background of the SPI is well described in McKee et al. (1993) . The basic point is to fit the monthly precipitation time series to the following Gamma function:

where Γ( ) is the Gamma function, y is the random variable and > 0 is a shape parameter. The next step is to obtain the cumulative probability distribution of precipitation and transform it to the standard normal random Z-variable with zero mean and unit variance. The cumulative probability of precipitation, F(x), is defined as follows:

(1)

(2) e a = RH 100 e 0 (T),

In this case β > 0 is a scale parameter and x > 0 represents the amount of precipitation. If x = 0 the precipitation distribution contains zeros.

As the gamma function is undefined for x = 0, the cumulative probability can be approximated as follows:

where p is the probability of a zero value. Finally, H(x) is transformed to the Z-variable, which is the SPI value.

Additionally, we modified the classical SPEI methodology but now incorporating Class A pan evaporation (E) and potential evapotranspiration (PET). In this case where D j represents the difference in Eq. (8) at the j-month, P j is the precipitation, E j is the Class A pan evaporation and PET vj is the virtual temperature-based potential evapotranspiration. Thus D j represents an approximate estimation of the available water in the soil-plant system.

PET v values were calculated using the Thornthwaite (1948) 

where L is the day length in hours, N is the number of days in the month, T v is the mean monthly virtual temperature, I is a heat index which is the sum of 12 monthly temperature indexes, i: and the m exponent is a function of I index:

As the Gamma distribution does not admit negative values, the application of Eq. (8) requires that D j > 0. Otherwise another distribution function such as the loglogistic function can be selected. In spite of that limitation, condition D j > 0 was met within the context of the present study. From our point of view Eq. (8) allows a more realistic analysis of water deficit. The time series derived from Eq. (8) was fitted to the Gamma distribution function instead of the log-logistic function which is used in the SPEI methodology. Tarquis et al. (2003) pointed out that practical applications of multifractal analysis involve the definition of a multifractal measure. We considered two multifractal measures. First the local temporal distribution of low and high fluctuations of the detrended SPI (1,6)/SPEI(1,6) m, RH and RSB signals (MF-DFA) and second the local covariance fluctuation between SPI (1,6)/SPEI(1,6)m and different meteorological variables (RH and RSB) through MF-DCCA. MF-DCCA (Zhou 2008) is an extension to two time series of the MF-DFA devised by Kantelhardt et al. (2002) for one time series. On the other hand, MF-DFA is a generalization of the Detrended Fluctuation Analysis (DFA) formulated by Peng et al. (1994) . Thus, we introduce here only the basic background on the MF-DCCA. Let us consider two time series x j and y j (j = 1,2, 3,..,S). We compute the deviation time series as follows:

Here x and y represent the mean of the signals x j and y j , respectively. Each deviation time series (Eqs. 12, 13) is divided into S π = int(S/π) non-overlapping windows of identical time scale π. In order to take into account the whole sequence, the same division procedure is carried out from the end of the time series. This generates 2S π non-overlapping sub-windows. Local trends are attained after least square fitting of a kth-order polynomial to each λ sub-window.

The detrended cross-covariance C 2 XY (π,λ) (λ = 1,2,…, Sπ) is obtained after subtracting the fitted polynomial from the original window as follows:

The qth-order detrended covariance function, C q (π) is computed as follows:

The particular case q = 0 implies

According to the fractal (multifractal) theory, C q (π) is a power-law function of the window size π. That is where the cross-correlation exponent h xy (q) is the generalized Hurst exponent. The particular cases h(q) and h xy (q) for q = 2 quantify the persistence (h or h xy > 0.5), randomness (h = 0.5 or h xy = 0.5) or antipersistence (h or h xy < 0.5).

The mass exponent spectrum, τ xy (q), in the MF-DCCA is defined as follows:

The local Hölder exponent (or singularity strength), α xy (q), is computed as follows:

The multifractal spectrum, f xy [α xy (q)], is identified as follows:

Equations (19-20) represent the Legendre transformation from the independent variables τ xy and q to the conjugate independent variables f xy and α xy (Arnold 1978; Baveye et al. 2008) .

First of all, we carried out a monofractal detrended fluctuation analysis (DFA) on the original time series before running MF-DFA and MF-DCCA. According to Eke et al. (2002) the signal is noise-like when the Hurst exponent (H r ) is within the range 0.2 ≤ H r ≤ 0.8. In such a case MF-DFA or MF-DCCA can be used without any data transformation. The time series represents a random walk-like structure when H r is within the range 1.2 ≤ H r ≤ 1.8. In that case the signal should be transformed. We selected the following multifractal parameters for analyzing both, individual and cross-correlation between different multifractal spectra:

(a) first of all the width of the singularity strength Δα which describes the intensity of time series fluctuations;

(16)

The height of the multifractal spectrum:

(c) The spectrum Asymmetry Index (AI hereafter) (Drożdż and Oswiecimka 2015) :

where α 0 is the Hölder exponent of the multifractal support (α for q = 0). where C 2 XY ( ) is the detrended cross-covariance function between X and Y time series. V X ( ) and V Y ( ) are the detrended variance functions corresponding to X and Y time series.

In the present study MF-DFA and MF-DCCA were performed using MFDFA.m and MFDXA.m MATLAB™ codes (Ihlen 2012) . We followed some Ihlen (2012) recommendations for selecting the appropriate parameters. One of them is to set the minimum segment size equal or larger than 10 samples. Other Ihlen (2012) suggestion is to select the maximum segment size below 1/10 of the sample size. Thus, we used the following input parameters with MFDXA.m code: k = 2 as the polynomial order for detrending, π min = 10 as the lower bound of scale, π max = 20 as the upper bound of the window size and a total of 10 segments between the above limits. We selected q = ± 10 as the upper and lower bounds of the q-order, respectively, with q res = 21 as the total number of q moments. The cross-correlation coefficient (ρ DCCA ) was computed using the DCCA Package developed by Prass and Pumi (2021) for R Language version 4.0.5 (R Core Team 2021). We used the following input parameters with the rhodcca function: the window size ranged from λ = 10 to λ = 50 with step = 1, k = 2 as (21) Δ = max − min (One multifractal spectrum)

Δ xy = xy,max − xy,min (Two cross−correlated multif ractal spectra);

the degree of the polynomial fit and overlap = FALSE as the choice for non-overlapping windows. Podobnik et al (2011) stressed the importance of computing the ρ DCCA range within which cross-correlations are significant at a given statistical level. We followed the same statistical approach used by those authors. That is, the time series correspond to independent and identically distributed (i.i.d) random variables. Thus, the time series could be seem as representing a Gaussian process with Hurst exponent H = 0.5. We generated 2000 i.i.d time series pairs from a Gaussian distribution each with length equal to the investigated time series (N = 480). The probability distribution function (PDF) of the ρ DCCA cross-correlation coefficient was computed for window sizes ranging from λ = 10 to λ = 50. Finally, we selected k = 2 and non-overlapping windows for the simulation process. Our choice differs from that used by Podobnik et al (2011) which considered k = 1 and overlapping windows. The detrending variance DFA(λ), detrending covariance DCCA(λ) and ρ (DCCA) were also computed using the DCCA Package developed by Prass and Pumi (2021) for R Language version 4.0.5 (R Core Team 2021). We selected 5% and 95% percentiles of the distribution as the lower and upper limits of the confidence interval, respectively. Therefore, posi-

The causes of multifractal structures were explored by randomly shuffling the original time series. We also generated ten surrogate time series using the non-linear approach developed by Schreiber and Schmitz (1996) . We included an additional analysis for estimating the index of stability of data distributions. The index of stability (Lévy index) was approximated using the Meerschaert and Scheffler (1998) equation: where is the Lévy index, γ0.557 is the Euler constant, x j (j = 1,2,…,N) are the time series data points and x is the time series mean. When 0 < < 2 the time series presents a heavy-tail distribution (e.g. Lévy-stable distribution) typical of sub-Gaussian behaviour. Even though Lévy-stable distributions do not have closed form formulations, Eq. (27) could be a useful approximation for practical proposals. The case = 2 corresponds to the normal distribution. This additional test would support any information derived from surrogates. Figure 2a shows the cumulative Gamma distribution of rainfall (P) and D = P-(E + PET) (Fig. 2b) , respectively. The continuous curve represents the theoretical Gamma function while filled squares correspond to data points. The

goodness-of-fit was significant at p < 0.05 in both cases. Therefore, we considered the Gamma distribution as an appropriate statistical model for computing SPI and SPEIm. Table 1 depicts the descriptive statistics of the original time series. The investigated time series showed long tail to the right based on 3rd order moment statistics (skewness) (apart from Class A pan evaporation). Class A pan evaporation (E) presented long tail to the left (S 3 = − 0.12). According to 4th order statistics rainfall (P), mean temperature and relative humidity (RH) were leptokurtic while pan evaporation and RSB presented platykurtosis. Based on that statistics, the distributions could be considered as non-Gaussian (K 4 = 0 for the normal distribution). Table 2 shows the skewness (S 3 ) and kurtosis (K 4 ) of residuals after time series detrending. The distribution of residuals preserved the non-Gaussian behaviour observed from the untransformed time series. Figure 3 illustrates the evolution of P and E + PET time series (Fig. 3a) , mean air temperature (T) and virtual (Fig. 3b) , RH (Fig. 3c) and RSB (Fig. 3d) . The spectral analysis (not shown) revealed that semi-annual and annual cycles were the main periodic components of each time series. Based on the trend lines one can observe that E + PET (Fig. 3a) , T and T v (Fig. 3b) increased in the studied time period . Precipitation (Fig. 3a) and RSB (Fig. 3d ) did not show a significant positive or negative trend. The difference ΔT = T v -T ranged from ΔT = 2.34 °C to ΔT = 3.68 °C within a moist environment with RH ranging between 81 and 99% (see Table 1 ). Only RH (Fig. 3c) showed a slight decrease in the studied time period. Figure 4a -b shows 1-month SPI and SPEIm (SPI1 and SPEI1m hereafter) (Fig. 4a) , 6 months SPI and SPEIm (SPI6 and SPEI6m hereafter) (Fig. 4b ) and the anomalies (SPEIm minus SPI) (Fig. 4c) . Values of SPI1 < − 1 are usually related to meteorological droughts while SPI6 < − 1 could be associated to agricultural droughts in the studied region (World Meteorological Organization 2012). In any case drought duration and severity depend on the consecutive number of months with SPI values within the range − 2 < SPI < − 1 or less. Fig. 4a and Table 3 that, according to the SPEI1m, available soil moisture could be slightly higher than that described by the conventional SPI1. This could be particularly true for March (1998) . For that particular month the standard SPI computed a value SPI1 = − 1.32 (moderately dry) while the SPEIm found SPEIm = 1.72 (very wet). Figure 5a shows the individual MF-DFA spectra for SPI1 and SPEI1m. Both spectra show left-hand deviations which suggest that local high fluctuations are more dominant than smaller fluctuations of the multifractal measure.

Note, however, that SPEI1m spectrum is more asymmetric than SPI1 spectrum. On the other hand, Fig. 5b presents the MF-DFA spectra corresponding to SPI6 and SPEI6m. Even though both multifractal spectra also show left hand deviations, the SPI6 spectrum is truncated to the right for negative q orders of the spectrum. Table 4 shows the more relevant multifractal parameters corresponding to the MF-DFA. These are the spectrum width (Δα), its height (Δf), the generalized Hurst exponent [h(2)], the Hurst exponent of the monofractal time series (H r ) and Asymmetry Index (AI). One can note that multifractal spectrum width, Δα, increases with the time scale. Furthermore, the spectrum width corresponding to SPEI1m and SPEI6m were broader than those corresponding to SPI1 and SPI6, respectively. Figure 6 depicts the SPI(1,6) and SPEI(1,6)m time series, local Hurst exponent (Ht) fluctuations computed at scale = 10 and scale = 20 and the probability distribution (Ph) of local Hurst exponents. The panels correspond to SPI1 (Fig. 6a) , SPEI1m (Fig. 6b) , SPI6 (Fig. 6c) and SPEI6m (Fig. 6d) , respectively. In all cases the charts of Ht versus time at scales 10 and 20 represent the local Hurst exponent fluctuation. The straight line at scales 10 and 20 is the modal Ht value. That value is also represented in the probability distribution of Ph versus Ht (lower chart of each panel). The maximum Ht value (max Ht) identifies the period with local fluctuation of smaller magnitude whereas the minimum Ht value (min Ht) is related to the period with local fluctuation of larger magnitude. Figure 7 represents the MF-DCCA spectra as related to the relationship between SPI1/SPEI1m versus RH (Fig. 7a) and SPI6/SPEI6m versus the same variable (Fig. 7b) . Contrary to the individual SPI(1,6) and SPEI(1,6) m, the spectra present right-hand deviations (AI < 0). In particular, the SPEI1m/RH multifractal spectrum presents a left-hand truncation for positive q order values. On the other hand, Fig. 8 depicts the MF-DCCA spectra as related to the relationship between SPI1/SPEI1m versus RSB (Fig. 8a ) and SPI6/SPEI6m versus the same variable (Fig. 8b) . In this case the main characteristic of both spectra is their broad width α max -α min . Consequently, one could consider that RSB increased the degree of multifractality of the SPI and SPEIm signals. Figure 9 shows the multifractal spectra corresponding to RH (MF-DFA), RSB (MF-DFA) and RH/ RSB (MF-DCCA). In the last case the spectrum is shifted in the direction of negative Hölder exponents (α XY ). RH and RSB spectra show right-hand deviations while the interaction between RH and RSB produced a left hand deviation spectrum. The multifractal parameters correspond to the individual SPI1/SPEI1m and SPI6/SPEI6m spectra derived from MF-DFA. Multifractal spectra width corresponding to SPEI1m and SPEI6m are wider than those computed from SPI1 and SPI6. Table 5 depicts the most important multifractal parameters derived from the MF-DCCA spectra. Random shuffling and surrogates are common methods for analyzing the sources of multifractality. Multifractal structures can be caused by long-term correlation and/or fat tail distributions. In this sense, Table 6 presents two basic multifractal parameters [h XY (2) and Δα XY ] computed after shuffling the original time series (SPI(1,6)/SPEI(1,6) m, RH and RSB) and generation of surrogate time series. One can note that, in general, h XY (2) for the original time series is larger than that corresponding to surrogates and random shuffling. Figure 10a -c depicts the fluctuations of the ρ(DCCA) cross-correlation coefficient as a function of different non-overlapping window size ranging from λ = 10 to λ = 50 months. Figure 10a corresponds to SPI1/SPEI1m versus RH and RSB, Fig. 10b represents SPI6/SPEI6m versus RH and RSB while Fig. 10c depicts the case RH versus RSB. Additionally, we also included the ρ c (DCCA) critical thresholds separating significant from no significant crosscorrelations. Figure 10d shows an example of the Probability Distribution Function of ρ(DCCA) corresponding to N = 480 and a window size λ = 50 months. It is almost symmetrical and does not differ from the normal distribution. We found a significant linear relationship between the critical ρ c (DCCA) cross-correlation coefficient and the window size λ for the 95% percentile.

The linear relationship between ρ c (DCCA) cross-correlation coefficient and the window size (λ) for the lower bound of the distribution (5% percentile) was negatively significant (R 2 = 0.983). Figure 11a -b illustrates the effect of random shuffling and surrogate data on the nonlinear structure of SPEI6m/ RH pair. Figure 11a presents the relationship between the generalized Hurst exponent h XY (q) and q for SPEI6m/RH original, SPEI6m/RH based on surrogates and SPEI6m/RH shuffled series. Note that h XY (q) versus q for the original series is larger than those corresponding to surrogates and random shuffling for all q orders. Figure 11b shows the multifractal spectrum f XY (α) versus α XY for the SPEI6m/RH pair corresponding to the original time series, surrogates and random shuffled series. It is noticeable that the original time series presents right-hand deviation (AI < 0) while randomly shuffled and surrogate time series are deviated to the left (AI > 0). Finally, Table 7 presents (28) c (DCCA) = 0.063(0.001) + 0.0019(0.00) (n = 41, R 2 = 0.997) 

At a first sight SPI1 and SPEI1m followed approximately the same regularity as shown in Fig. 4 . Neither SPI1 nor SPEI1m showed significant positive or negative trends (Fig. 4a) . On the contrary both SPI6 and SPEI6m presented a consistent negative trend towards dry episodes in the future (SPI6/SPEI6m < 0) (Fig. 4b) . The anomalies (differences) are better visualized in Fig. 4c mainly for SPI1 minus SPEI1m. The incorporation of pan evaporation and virtual temperature-based PET data into the SPEI computation provides a better estimation of the available moisture in the soil-plant-atmosphere system. That is possibly a reason of why virtual temperature is a physical parameter used in the computation of crop coefficients. It could be due to the inclusion of RH, actual vapour pressure and atmospheric pressure, P(z), into the SPI formulation. This could be of interest for assessing meteorological and agricultural droughts. Table 3 ). This particular finding agreed well with a severe drought in the Amazon basin in 2010. Thus, from our point of view it is difficult to identify drought events based only on precipitation data. In general, the SPI6 and SPEI6m anomalies (differences) were irrelevant. Nevertheless, this finding is consistent with recent studies conducted by Pei et al. (2020) . Those authors found that the differences between SPI and SPEI decreased with increasing timescale. This way, for short-term SPI timescales (e.g. 1 month scale) the traditional SPI did not match its modified counterpart for several time periods. In summary our approach shows some differences as compared to that developed by Li et al. (2017) . While those authors used the SPEI methodology with precipitation and pan evaporation, we added virtual temperature-based PET to the standard SPEI computation.

We are well aware that our approach could fail for other climate types were D j < 0. In such a case the log-logistic or Pearson Type III distribution could be a suitable choice.

The multifractal spectrum width (Δα) increased after including pan evaporation and virtual temperature-based PET for SPEIm computations as shown in Table 4 and Fig. 5a ,b. Furthermore, Δα also increases when the SPEI or SPEIm increased from 1 to 6 months timescales. The larger influence of the modification was exposed by the SPEI6m (Δα = 0.943). This indicates, to some extent, that the SPEI6m can reveal a higher degree of multifractality than the classical SPI6 model. That is, the SPEI6m resolved a wider range of the probability distribution of wet/dry events. Therefore agricultural drought complexity could be influenced by precipitation deficit, virtual temperature and real evaporation. Ogunjo (2021) also observed an increase of Δα with the timescale in tropical locations across Nigeria. For example that author found spectra width ranging from Δα = 0.158 for 1-month SPI timescale to Δα = 0.397 for 24 months scale, respectively. Furthermore, Adarsh et al. (2019) found multifractal spectra width of Δα = 0.27 for 3 months, Δα = 0.30 for 6 months and Δα = 0.39 for 12 month SPI timescales. The increase of the multifractality degree (Δα) with the aggregation timescale seems to be a characteristic feature of drought index scaling. The Δf value corresponding to the SPEI1m signal was approximately the same as that corresponding to the SPI1. Thus, both signals spanned approximately the same Hausdorff dimension range.

However, some f(α) values were negative for SPEI1m, SPI6 and SPEI6m (− 0.411 ≤ f(α) ≤ − 0.115) [left-hand tail of Fig. 5a,b ]. This occurred mainly within a narrow range from α = 0.201 to α = 0.242 and α = 0.74 to α = 0.80, respectively.

While those negative f(α) values are usually discarded, from a physical stand point f(α) < 0 for some α values were well described theoretically by Mandelbrot (1991) . That author suggested that negative f(α) could be due to fluctuations associated to the finite sample size and/or sampling variability which fits well within the context of hydrometeorology. On the other hand, Mandelbrot (1990) also stressed that negative f(α) could be essential for understanding turbulence and Diffusion Limited-Aggregation (DLA). Both physical phenomena are also part of the atmospheric dynamics. Nevertheless, it is needed a larger sample which he called supersamples. The 2nd-order Hurst exponent, h(2), was larger than 0.5 which indicates that long-term correlation is an important factor governing the SPI dynamics (Table 4) . Eke et al. (2002) criterion confirmed that each time series, modeled as a monofractal signal, represents a like-noise process with long-term correlation (H r > 0.5) ( Table 4 ). The Asymmetry Index (AI) was positive for SPI(1,6) and SPEI(1,6)m which is consistent with left-hand deviation of spectra. This indicates that higher SPI(1,6) and SPEI(1,6)m covariance fluctuations at local time windows are more significant than smaller fluctuations at the same time scales. That is, dry or wet events at short-time scales were more frequent within the investigated time period. This agrees well with investigations conducted by Bunde et al. (2005) . Those authors have point out that long-term correlation (H > 0.5) is a natural mechanism for the clustering of extreme events. Note, however, that the individual multifractal spectrum corresponding to RH and RSB was deviated to the right hand (AI < 0) (Table 4 ). This particular finding suggests that RH and RSB produce low-amplitude oscillations of the covariance function. Thus the multifractal pattern of SPI(1,6) and SPEI(1,6)m could be influenced by such low oscillations in their interactions with RH and RSB.

The largest local Hurst exponent, H t = 0.963, at the scale = 10 was attained at normal conditions (SPI1 = 0.37) around the year 1981 (Fig. 6a) . However, the lowest local Hurst exponent, H t = 0.357, at the scale = 10 was obtained under extremely dry conditions (SPI1 = − 2.55). That value corresponded approximately to September 2011. It was a time period of local low fluctuations as observed from the upper chart in Fig. 6a . Such a local anti-persistency can indicate an inverse pattern in the next years as observed for the years 2013 and 2014. The minimum local Hurst exponent (H t = 0.357) at the scale = 20 was observed from an inverse peak of low amplitude at moderately dry conditions (SPI1 = − 1.31) around 2011. The maximum local Hurst exponent, H t = 0.807, at the scale = 20, was reached under a normal condition state (SPI1 = 0.93) around the year 1974 (top chart in Fig. 6a . In the case of SPEI1m, the maximum and minimum H t values matched similar time periods as SPI1 and under similar wet/dry conditions (Fig. 6b) . Consequently, one can depict two scenarios. First of all the largest H t value (H t > 0.5, long-term correlation) was reached under normal conditions. Second, the minimum Hurst exponent (H t < 0.5, anti-persistence) was attained under dry conditions. This occurred for scales 10 and 20 and SPI1/SPEI1m. Thus, the dynamics of Hurst exponents was almost the same for SPI1 and SPEI1m. The modal H t value was approximately the same for SPI1 and SPEI1m (Fig. 6a,b) .

The dynamics of Hurst exponents for SPI6 and SPEI6m seems to be different when compared to SPI1 and SPEI1m as shown in Fig. 6c -d. The maximum local Hurst exponent (H t = 1.302) at the scale = 10 was attained at normal conditions (SPI6 = 0.23) around the year 2003 (Fig. 6c) . On the other hand, the lowest local Hurst exponent (H t = 0.897) at the scale = 10 was also obtained under normal conditions (SPI6 = 0.56) in 1984. The smallest local Hurst exponent (H t = 0.885) at the scale = 20 was observed at normal conditions (SPI6 = 0.13) around 2012. The largest Hurst exponent (H t = 1.587) at the scale = 20, was also observed under normal condition (SPI1 = 0.64) around the year 2008 (top chart in Fig. 6c ). Thus, maximal and minimal Hurst exponents were obtained under normal hydrometeorological conditions. The SPEI6m shows some differences as compared to SPI6 (Fig. 6d) . The maximum Hurst exponent (H t = 1.644) at the scale = 10 was also attained at normal conditions (SPEI6m = 0.43) around the year 2004 (Fig. 6d) . At the same time the lower local Hurst exponent (H t = 0.875) was also computed under normal conditions (SPEI6m = 0.59) in 1984. The largest local Hurst exponent (H t = 1.602) at the scale = 20 was reached under severely wet conditions (SPEI6m = 1.98) around the year 1974. However the smallest H t (H t = 0.887) was also computed at normal conditions (SPEI6m = 0.56) but in the year 1984 (Fig. 6d) . Accordingly the nonlinear dynamics of Hurst exponent fluctuation corresponding to SPI6/SPEI6m is different from that observed for SPI1/SPEI1m. Figure 7a shows the MF-DCCA spectrum for SPI1 and SPEI1m while Fig. 7b refers to the SPI6 and SPEI6m versus RH in both cases. In terms of Δα XY and Δf XY , the SPEI1m and SPEI6m versus RH revealed approximately the same degree of multifractality as the conventional SPI1 and SPI6 (Table 5 ). The Asymmetry Index (AI XY ) was negative in all cases which suggests right-hand deviation of the spectra (Fig. 7a,b) . Therefore the low amplitude fluctuation of the covariance function is due mainly to the influence of RH. In general, the influence of RH on low fluctuations at local time scales was the same for both SPI timescales. It is interesting that the left-hand multifractal spectra for SPI1/ RH and SPEI1m/RH are located within a range of negative Hölder exponents (− 0.36 ≤ α XY ≤ − 0.175) (Fig. 5a) . In general negative Hölder exponents have been theoretically defined within the p-exponent and p-leader formalism (Jaffard et al. 2016) . In this case the p-exponent meets the condition 0 < p ≤ − 1/α (Leonarduzzi et al. 2014) . Negative regularity has been recently computed from real-world data. Even though it was associated to an impulse-like singularity of wind speed, it is necessary a p-leader wavelet-based analysis. From our point of view the only useful information is that α XY < 0 could be associated to infrequently occurring hydrometeorological events. In this case the local crosscovariance function drops within some scales. RH showed a larger influence on SPI1 as compared to SPEI1m. In that case the multifractal spectrum resolved a wider range of the probability distribution for q > 0. Figure 7b depicts the multifractal interaction between SPI6 and SPEI6m versus RH. The whole curve is positioned within the range of positive Hölder exponents (0.15 ≤ α XY ≤ 0.78). Both multifractal spectra presented right hand deviation (AI XY < 0, q < 0). This suggests that smaller fluctuations of the covariance function were significant at the 6 months time scale. It is consistent with the fact that extreme hydrological events are less predictable for higher q-order. Ramanathan and Satyanarayana (2019) have pointed out that predictability drops for larger q-order moments of the covariance function as higher q-order moments take into account rarer events. Figure 8a illustrates the multifractal spectrum for SPI1/ RSB and SPEI1m/RSB whereas Fig. 8b refers to the SPI6/ SPEI6m versus RSB, respectively. Regarding Fig. 8a , approximately the left hand part of both multifractal spectra is located in the region of the so-called p-exponents (negative Hölder exponents). On the other hand some f XY (α) values were negative at both tails of the spectra (Fig. 8a) . Furthermore the AI for each spectrum was AI XY < 0 (Table 5) . Thus, the RSB can influence the SPI1 and SPEI1m in three ways. First, long-range small fluctuations of the SPI1/RSB covariance function were dominant (AI XY < 0). Second, those low-amplitude fluctuations could suggest clustering of dry/wet states (α XY > 0, q < 0) and impulse-like singularities could be also likely (α XY < 0, q > 0). Third, some low fluctuations could be due to the finite sample size (f XY < 0) and/or the effect of atmospheric turbulent flux. Figure 8 (b) suggests different information as compared to Fig. 8a . In this case Δα > 0 within the whole range − 10 ≤ q ≤ 10. On the other hand, the f XY (α) < 0 values spanned a smaller range of the right hand tail of the spectra than those shown in Fig. 8a . The AI for each spectrum was AI XY < 0 (Table 5) . Thus, both multifractal spectra also suggest clustering (Δα > 0) and low fluctuations of the covariance function (AI XY < 0). In general the multifractal width of SPI(1,6)/SPEI(1,6)m versus RH ranged from Δα = 0.455 to Δα = 0.574. On the other hand the multifractal strength of SPI(1,6)/SPEI(1,6)m versus RSB varied between Δα = 0.692 to Δα = 0.735. This signifies that SPI(1,6) and SPEI(1,6)m were more sensitive to RSB than RH with higher long-term correlations. One can note that Δα increased after performing the SPEI modification. According to Agbazo et al. (2019) the larger Δα is, the more heterogeneous the probability distribution is. It is likely to advance a mixture of short and long drought periods. The left-hand tail of the RH/RSB singularity spectrum is shifted toward p-exponents (negative α XY ) ( Fig. 9) with high fluctuations of the covariance function (AI XY > 0). Thus, intermittence and/or fast change in the atmospheric pattern can be also associated with the combined effect of RH and RSB. This effect drops as time scale (window size) increases.

The ρ DCCA cross-correlation coefficient is an interesting avenue for exploring the evolution in time of a drought indicator and its interaction with different meteorological variables. For example, in Fig. 10a one can observe significant long-term cross-correlations between SPI1/SPEI1m versus RH at all considered time scales. All the ρ DCCA values are larger than the critical value represented by the 95% percentile of the distribution function. In particular, the maximum ρ DCCA value was computed for the 38-month time scale (approximately 3 years) (ρ DCCA (SPI1) = 0.643 and (ρ DCCA (SPEI1m) = 0.656, respectively). On the contrary, SPI1/SPEI1m versus RSB showed negative crosscorrelations between 5 and 25 months' time scales. After 25 months' window scale the cross-correlation coefficient between SPI1/SPEI1m versus RSB could be considered as significant for time scales λ = 35, λ = 40 and λ = 45 months (Fig. 10a) . The dotted lines represent an approximation to the critical ρ c value (ρ c = 0.256) suggested by Table 2 in Podobnik et al. (2011) . Thus, according to our study RH was positively cross-correlated with the occurrence of drought events (SPI1/SPEI1m) at all window scales. This could be due to the relationship between RH and rainfall occurrence. On the contrary, RSB was negatively cross-correlated with SPI1/SPEI1m only at specific short-term time scales. Figure 10b shows a different picture as compared to Fig. 10a . The ρ DCCA cross-correlation coefficient between SPI6/ SPEI6m versus RH and RSB could be considered as not significant. However, we found a large positive peak around the 35 months' time scale (approximately 3 years) for SPEI6m versus RH (ρ DCCA = 0.561). That positive maximum coincided with an extreme negative minimum ρ DCCA value for SPEI6m versus RSB (ρ DCCA = − 0.439) (Fig. 10b) . We found a negative RhoDCCA cross-correlation coefficient for RH/ RSB within each time lag (Fig. 10c ). An extreme negative value was detected around 12 months' time scale (1 year) (ρ DCCA = − 0.862). The critical ρ c (DCCA) value corresponding to the lower limit of the PDF (5% percentile) ranged from ρ DCCA = − 0.087 (λ = 10 months) to ρ DCCA = − 0.159 (λ = 50 months). Note that all negative ρ DCCA cross-correlation values depicted in Fig. 10c fall below that range of the distribution function. Thus, negative ρ DCCA cross-correlations between RH/RSB pair were all significant. As a particular example, the critical points of the PDF for a window size λ = 50 ranged from ρ DCCA = − 0.1591 (5% percentile) to ρ DCCA = 0.1568 (95% percentile) (Fig. 10d) . From our point of view the ρ DCCA cross-correlation coefficient allows one to identify those window scales within the SPI or SPEIm timescale where the influence of any meteorological variable is stronger. The shuffling process removed long-term correlations from SPI1/RH, SPI6/RH, SPEI6m/RH and SPI1/RSB spectra (Table 6 ). Then for those four cases long-term correlation of small and large fluctuations of the covariance function could be the dominant cause of multifractality. Figure 11a ,b illustrates the case of SPEI6m/RH spectrum. The h(q) versus q relationship was the smallest for the randomly shuffled data while h(q) versus q for surrogates was also smaller than the original (Fig. 11a) . Thus a long-tail distribution also contributes, in part, to the multifractal structure. The Δα XY difference between the original signal and surrogates was larger than the difference between the original and random shuffled data for SPEI1m, SPEI1m/RSB, SPEI6m/RSB and RH/RSB spectra. This also indicates that long-tail distribution (e.g. Levy distribution) could be a cause of the multifractal structure in the previous cases. The random shuffling and surrogates transformed the symmetry of the original multifractal spectrum. Note that the original SPEI6m versus RH spectrum presented a negative AI XY (AI XY = − 0.436). On the contrary, randomly shuffled data and surrogates showed positive AI XY values (AI XY = 0.32 and AI XY = 0.53, respectively) as depicted in Fig. 11(b) . Furthermore, the singularity spectra corresponding to random shuffling and surrogates have right truncation and long left tails. The application of Meerschaert and Scheffler (1998) approach allowed us to identify the meteorological variables determining the degree of multifractality and long-memory effect. According to Table 7 precipitation (κ = 0.824), RSB (κ = 1.008), pan evaporation (κ = 1.078) and RH (κ = 1.603) are the main sources of multifractal behaviour. This allows us to address the following three basic issues: First, the probability of occurrence of extreme events is higher than that predicted by the normal distribution. Second, the duration and intensity of the extreme events (e.g. droughts) are inversely proportional to the frequency of occurrence. Third, the Meerschaert and Scheffler (1998) Eq. (27) could be used in similar studies for conclusive evaluation of heavy tail distribution as a source of multifractality.

The present study compared the classical SPI, which is based only on precipitation data, with a modification to the SPEI (SPEIm) introducing Class A pan evaporation and virtual temperature-based potential evapotranspiration (PET) into the model distribution function. The study was conducted using 1 months and 6 months SPI and SPEI timescales. MF-DFA and MF-DCCA were performed on the traditional SPI and its SPEIm counterpart. Some concluding remarks are obtained as follows:

(a) The integration of pan evaporation and virtual temperature-based PET into the SPEI computation can improve our perception of the available moisture in the soil-plant-atmosphere system. However, our approach seems to be more effective for humid climates (e.g. tropical environments) where precipitation is larger than E + PET. When D i < 0 other statistical distribution (e.g. log-logistic) function should be used instead of the Gamma distribution. (b) The SPEIm did not match the traditional SPI for several periods at the 1-month SPI time scale. In particular the SPEI1m did not agree with the SPI classification for some time periods. (c) The SPEIm allowed us to detect a more heterogeneous probability distribution than that derived from the traditional SPI. Short-and long-time scale fluctuations of the SPEIm were significantly cross-correlated with RH and RSB. (d) The nonlinear dynamics of local-scale invariance of SPI(1,6)/SPEI(1,6)m time series experienced structural changes approximately at the same time periods independent of the window scale. (e) Finally, the ρ DCCA cross-correlation coefficient offers an interesting avenue for exploring extreme events occurrence (e.g. droughts) and the meteorological variables connected to such events. That constitutes a topic difficult to investigate with more traditional statistical methods.

IM and NV computed the SPI and SPEIm; HM, IM and NV reviewed and edited and HM wrote the manuscript. All the authors and academic authorities revised and approved the present manuscript.

Funding No funding was received for conducting this investigation.

Data availability All the data used in the present study are available upon request.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Improved ground hydrology calculations for global climate models (GCMs): soil water movement and evapotranspiration

Multifractal characterization of meteorological drought in India using detrended fluctuation analysis

Multifractal behaviours of daily temperature time series observed over Benin synoptic stations (West Africa)

Crop evapotranspiration: Guidelines for computing crop water requirements. FAO Irrigation and drainage paper

Multifractal analysis of meteorological time series to assess climate impacts

Biophysical chemistry of fractal structures and processes in environmental systems, IUPAC series on analytical and physical chemistry of environmental systems

Standardized precipitation evapotranspiration index (SPEI) revisited: parameter fitting, evapotranspiration models, tools, datasets and drought monitoring

Lon-term memory: a natural mechanism for the clustering of extreme events and anomalous residual times in climate records

Thermodynamic factors in evapotranspiration

Future meteorological droughts in Ecuador: decreasing trends and associated spatio-temporal features derived from CMIP5 models

Social protection and agriculture: introduction to the special issue

The effect of neglecting the virtual temperature correction on CAPE calculations

Detecting and interpreting distortions in hierarchical organization of complex time series

Fractal characterization of complexity in temporal physiological signals

A century of nonlinearity in the geosciences

The physics of climate variability and climate change

The Lincoln declaration on drought indices: universal meteorological drought index recommended

Multifractal analysis of the drought area in seven large regions of China from 1961 to 2012

Introduction to multifractal detrended fluctuation analysis in Matlab

Fundamentals of atmospheric physics

) p-exponent and p-leaders, part I: negative pointwise regularity

Multifractal detrended fluctuation analysis of nonstationary time series

Extreme events in surface winds: predicting turbulent gusts

Extending multifractal analysis to negative regularity: p-exponents and p-leaders

A revised drought index based on precipitation and pan evaporation

Multiscaling properties of drought conditions in China

Modification of the standardized precipitation evapotranspiration index for drought evaluation

Negative fractal dimensions and multifractals

Random multifractals: negative dimensions and the resulting limitations of the thermodynamic formalism

The relationship of drought frequency and duration to time scales

A simple robust estimator for the thickness of heavy tails

Meteorological complexity in the Amazonian area of Ecuador: an approach based on dynamical system theory

A review of drought concepts

On the dynamic interpretation of the virtual temperature

Exploring the multi-fractal nature of the air flow and pollutant dispersion in a turbulent urban atmosphere and its implications for long range pollutant transport

Multifractal properties of meteorological drought at different time scales in a tropical region

Updated world map of the Köpen-Geiger climate classification

Comparative analysis of drought indicated by the SPI and SPEI at various timescales in inner Mongolia

Mosaic organization of DNA nucleotides

Multifractal Detrended fluctuation analysis of temperature reanalysis data over Greece

Multifractal detrended crosscorrelation analysis of wind speed and solar radiation

Statistical tests for power-law cross-correlated processes

On the behavior of the DFA and DCCA in trend-stationary processes

R: A Language and Environment for Statistical Computing, version 4.0.5. R Foundation for Statistical Computing

Higher-order statistics based multifractal predictability measures for anisotropic turbulence and the theoretical limits of aviation weather forecasting

Heat and drought 2003 in Europe: a climate synthesis

Multifractal cross correlation analysis of agro-meteorological datasets (including reference evapotranspiration) of California

Improved surrogate data for nonlinearity tests

Co-benefits and synergies between urban climate change mitigation and adaptation measures: a literature review

Scaling and multiscaling of soil pore systems determined by image analysis

An approach toward a rational classification of climate

Spatiotemporal drought variability on the Mongolian Plateau from 1980-2014 based on the SPEI-PM, intensity analysis and Hurst exponent

A multiscalar drought index sensitive to global warming: the standardized precipitation evapotranspiration index-SPEI

Cross-correlations between bacterial food borne diseases and meteorological factors based on MF-DCCA: a case in South Korea

Agrometeorology related to extreme events

Standardized precipitation index user guide. WMO-No 1090

The global climate in 2015-2019

DCCA cross-correlation coefficient: quantifying level of cross-correlation

Multifractal analysis of measure representation of flood/drought grade series in the Yangtze Delta, China, during the past millennium and their fractal model simulation

Predicting climate change mitigation and adaptation behaviors in agricultural production: a comparison of the theory of planned behavior and the Value-Belief-Norm Theory

Multifractal detrended cross-correlation analysis for two nonstationary signals

The authors are grateful to INAMHI (Instituto Nacional de Meteorologia e Hidrología, Ecuador) for providing free access to the meteorological data set and its quality control process. Dr. Humberto Millán wants to express gratitude to the Amazonic State University (Puyo, Pastaza) and several undergraduate students for the valuable support during his stage as an invited professor of Physics. We also want to dedicate this modest contribution to the memory of Sr. John T. Houghton (an outstanding Oxford Atmospheric Physicist and IPCC leader) who died as a consequence of COVID-19 pandemic.Author contributions HM and IM designed research; IM, HM and NV contacted INAMHI, acquired data sets and performed research; HM, and IM analysed data; HM conducted the MF-DFA and MF-DCCA;