Multifractal Cross Correlation Analysis of Agro-Meteorological Datasets (Including Reference Evapotranspiration) of California, United States

: The multifractal properties of six acknowledged agro-meteorological parameters, such as reference evapotranspiration (ET0), wind speed (U), incoming solar radiation (SR), air temperature (T), air pressure (P), and relative air humidity (RH) of ﬁve stations in California, USA were examined. The investigation of multifractality of datasets from stations with di ﬀ ering terrain conditions using the Multifractal Detrended Fluctuation Analysis (MFDFA) showed the existence of a long-term persistence and multifractality irrespective of the location. The scaling exponents of SR and T time series are found to be higher for stations with higher altitudes. Subsequently, this study proposed using the novel multifractal cross correlation (MFCCA) method to examine the multiscale-multifractal correlations properties between ET0 and other investigated variables. The MFCCA could successfully capture the scale dependent association of di ﬀ erent variables and the dynamics in the nature of their associations from weekly to inter-annual time scales. The multifractal exponents of P and U are consistently lower than the exponents of ET0, irrespective of station location. This study found that joint scaling exponent was nearly the average of scaling exponents of individual series in di ﬀ erent pairs of variables. Additionally, the α -values of joint multifractal spectrum were lower than the α values of both of the individual spectra, validating two universal properties in the MFCCA studies for agro-meteorological time series. The temporal evolution of cross-correlation determined by the MFCCA successfully captured the dynamics in the nature of associations in the P-ET0 link.


Introduction
Characterization of agro-meteorological time series is a major concern for irrigation and water resource planners. Agro-meteorological time series often possess non-linear, non-stationary and complex temporal scaling characteristics [1], which make the accurate prediction of their course difficult. Any improved understanding on the fluctuations and scaling characteristics of such time series may enhance the modeling efforts aimed at a fuller understanding of their structures. Many of agro-meteorological time series possesses multiple time scales ranging from sub-daily to yearly, ET0. In this context, the determination of correlation between different causal meteorological variables with ET0 is of great significance in agro-meteorological modeling. However, such associations were never investigated in the past in a multiscale-multifractal perspective.
Therefore, the specific objectives of this work include: (i) to investigate the multifractal properties of ET0 and meteorological variables like wind speed, air pressure, solar radiation, air temperature, relative air humidity at five stations located at distinctly different terrain conditions in California, USA; (ii) to investigate the power law correlations of different meteorological variables with ET0 using the novel MFCCA approach. The next section presents the theoretical details of MFDFA and MFCCA methods. The details of data used in the study are presented in the section thereafter. Subsequently, results of MF analysis and MFCC are presented along with relevant discussions. The major conclusions drawn from the study are presented in the last section.

Materials and Methods
In this study, firstly MFDFA method is applied to examine the multifractality of different agro-meteorological time series of California. Subsequently, MFCCA method is used to investigate the mutual association of different contributory meteorological variables with ET0. As the nonlinear properties (i.e., the long-term volatility correlations and the width of the multifractal spectrum) are strongly affected by the periodicities in agro-meteorological time series, the periodic seasonal trends in all the studied time series have to be filtered out. To remove them, the Seasonal and Trend decomposition using Loess (STL) method developed by Cleveland et al. [70] was used before performing the MFDFA and MFCCA analyzes.

MFDFA Method
The different steps of MFDFA are: (1) Determine the accumulated deviation of the time series signal X (x1, x2 . . . xL) from its mean (so called 'profile') as: where i = 1, 2,..., L, with mean x and length L. (2) Identify different independent sub-series of different lengths (called scale s) both in forward and backward directions. (3) Fit polynomial to each sub-series and find the difference between original sub-series and its fitted polynomial as: and: Here y υ (i) is the polynomial used for fitting in fragment υ. The non-linear behavior of the plot between h(q) and q indicates the multifractality of the time series. For monofractal time series h(q) is independent of q. Multifractal behavior (a significant dependence of h(q) on q) will manifest itself only if small and large fluctuations scale differently in the time series. For positive values of q, h(q) describes the scaling behavior of the segments with large fluctuations, whereas for negative values of q, the scaling behavior of the segments with small fluctuations is described. While the precision of h(q) estimation depends mostly on the length of the time series, for large negative values of q (very small fluctuations) the precision of measurement is also starting to play a role. A detailed analysis on the precision of h(q) estimation for the MFDFA method was performed by López and Contreras [71]. The GHE value for q = 2 is considered to be equivalent to the classical Hurst exponent (H) [27,28]. From h(q) other useful exponents like mass exponent τ(q) or singularity exponent f (α) can be derived as follows: The plot of f (α) vs. α is called as singularity spectrum, which will be inverted parabolic in shape with apex at unity for multifractal series. The base width of the spectrum (w = α max − α min ) is an indicative of the strength of the multifractality in such a way that larger width indicates higher multifractal strength and vice versa. The symmetry of the spectrum (symmetric, left or right-skewed) indicates the dominancy for high (or low) fluctuations. It is quantified by a parameter called asymmetry index (AI), computed based on the width of right and left wings of the parabola as: with ∆α left = α 0 − α min and ∆α right = α max − α 0 , respectively. If AI is zero, the spectrum is symmetric, positive values indicates a right skewed and negative values indicates a left-skewed spectrum. Negative values of AI suggest the existence of complex structures at the level of large fluctuation amplitudes in the respective series; they also suggest that the time series are characterized by a multifractal structure, which is insensitive to local fluctuations with small amplitudes. This means that a low weight can be attributed to low fractal exponents, which in turn suggests the frequent occurrence of extreme events [72]. The value of α at q = 0 (indicated as α 0 and called as Holder exponent) infers the complexity of the time series. The mathematical expressions for the computation of multifractal properties and the guidelines for selection of algorithm specific control parameters such as polynomial order, range of moment order, scaling range (s min − s max ) etc. can be can be found in literature [73].

MFCCA Method
The different steps of MFCCA computational procedure can be described as follows: (1) For two time series x i and y i (I = 1, 2, . . . , N), determine the profiles as: and: Atmosphere 2020, 11, 1116 6 of 24 where j = 1, 2, . . . , N and x and y are the mean of the two time series x i and y i . (2) Divide each profile P x (j) and P y (j) into N s non-overlapping segments both in progressive and retrograde directions, to avoid any omission of time series data from the head or tail end of the series. (3) For each 2N s segments, compute local trend of both "profiles" P x (j) and P y (j) by fitting polynomials p X,υ m (j) and p Y,υ m (j) of appropriate order m. The subtraction of the fitted polynomial from the original segment gives the covariance as given below: where s is the length of υ segment. (4) Calculate detrended covariance by summation over all segments: (5) Check if F q XY (s) behaves as a power-law function of s (the scaling behavior), where s is the segmental sample size: The scaling exponent λ XY (q) is similar to the generalized Hurst exponent h(q) in MFDFA and it can be obtained by observing the slope of log-log plot of F(s) versus s by ordinary least squares. The ratio between the detrended covariance function F q XY (s) obtained from MFCCA and the detrended variance functions F q X (s) and F q Y (s) obtained from MFDFA defines the detrended cross correlation ρ XY (s) as: Theoretically, the value of ρ XY (s) ranges between −1 ≤ ρ XY (s) ≤ 1. The MFCCA analysis facilitate the estimation of scale dependent correlation between two candidate time series, which can provide better insight into the physical association between the variables.

Study Area and Data
The multifractal characteristics of time series changes with latitude, altitude and distance from the coast [74]. Therefore, it is better to consider datasets from different terrain conditions and elevations. Considering this, the daily data of five automated weather stations of California, USA, namely Dagget (latitude 34 [75], with hot, dry summers, and cooler wetter winters with lesser precipitation than typical Mediterranean climate. All the stations considered are located at radial distances within 180-250 km from the city of Los Angeles. Dagget city is located at north-east of Los Angeles by 200 km at the interior part of Southern California with no oceanic proximity. According to the KCC, Daggett has a semi-arid climate. There is very little rainfall year-round and diurnal temperature variation is high. In the winter, Atmosphere 2020, 11, 1116 7 of 24 nights are often very cold while days are mildly cool. Meanwhile, summer nights are mildly warm while days can be extremely hot. The major land resource areas of the Daggett are Sonoran Basin and Range (Arizona, California, and Nevada, western range and irrigated region). The evaporation losses in Daggett site are moderately high due to the high temperature. The Bakersfield station is 180 km to the north of Los Angeles, situated in the center of Kern County, located in the extreme southern end of the Great San Joaquin Valley, and is partially surrounded by a horseshoe-shaped rim of mountains. The climate in Bakersfield has a hot desert type climate with very hot, dry summers, and winters that consist of mild days with chilly/cold nights. San Diego, with approximate height of 4 m above m.s.l., is located about 200 km to the south of Los Angeles on the coast of Pacific ocean. The basic climate is typical Mediterranean as per KCC with hot, sunny, and dry summers, and cooler, wetter winters. However, the city is much more arid, and winters are still dry on comparing with most other zones with this type of climate. The wet winter season of San Diego is influenced by the El Nino Southern Oscillation (ENSO) and it receives more winter storms with warmer and more humid conditions. During the La Nina phase, San Diego becomes drier with cooler and less humid conditions. Santa Maria is located in the central coast of California in the northern Santa Barbara county. It is approximately over 240 km north west of Los Angeles. The region experiences a cool Mediterranean climate (KCC typical of coastal areas of California north of Point Conception). The climate is mostly sunny, refreshed by the ocean breeze. The attitudes of these stations are 586, 151, 77, 30 and 4 m from the mean sea level, respectively. This helped to capture the differences in multifractal behavior with the altitudinal differences. Long daily records of air temperature (T), solar radiation (SR), wind speed (U), air pressure (P), relative air humidity (RH), and reference evapotranspiration (ET0) of 1961-1990 period are used for multifractal analysis. FAO−56 Penman-Monteith equation was used by United States Environmental Protection Agency for the estimation of ET0. The map showing locations of different stations is provided in Figure 1.  Figure 1. The measured daily climatic data for these stations were downloaded from the US EPA web server [76]. The respective time series are plotted in Figure 2. The quality control is a major prerequisite for using meteorological information, so it is highly necessary to ensure that the record structure of the data is correct and complete before proceeding with modeling of any time series and there are definite procedures to perform the quality control of the datasets [77][78][79]. The measured daily climatic data for these stations were downloaded from the US EPA web server [76]. The respective time series are plotted in Figure 2. The quality control is a major prerequisite for using meteorological information, so it is highly necessary to ensure that the record structure of the Atmosphere 2020, 11, 1116 8 of 24 data is correct and complete before proceeding with modeling of any time series and there are definite procedures to perform the quality control of the datasets [77][78][79]. However, in this case, the network is highly sophisticated and data was published by USEPA after rigorous quality examinations and used widely by other researchers [48,80]. However, along with the visual examination the data was checked for their continuity and also the chances of outliers was examined. No such anomalous observations and their replacement were needed and separate homogeneity analysis was not warranted in this case. The basic statistical parameters of the different datasets are provided in Table 1. However, in this case, the network is highly sophisticated and data was published by USEPA after rigorous quality examinations and used widely by other researchers [48,80]. However, along with the visual examination the data was checked for their continuity and also the chances of outliers was examined. No such anomalous observations and their replacement were needed and separate homogeneity analysis was not warranted in this case. The basic statistical parameters of the different datasets are provided in Table 1.

MFDFA of Agro-Meteorological Time Series
The multifractality of the five contributing time series (T, U, SR, RH, P) and ET0 time series of all the five stations were examined using MFDFA method. The minimum and maximum scale ranges can be fixed appropriately by following the guidelines stated in literature [72,73]. The minimum scale s min greater than (m + 2), where 'm' is the polynomial order, is a widely accepted thumb rule. The maximum scale s max is suggested to be chosen as up to L/10 [27]. Here, a rigorous trial and error exercise is made in the selection of the scale range considering the physics of the problem. The fluctuation functions are developed by considering the statistical moment orders from −4 to 4 are fixed and the first order polynomial (m = 1) was chosen for invoking MFDFA [41]. The fitting is considered to be acceptable on getting an R 2 ≥ 0.98 [81]. Thus a scaling range of [8:600] days was found to be appropriate to capture the scaling behavior, which corresponds from weekly to inter-annual period.
Different graphical representations such as GHE plot, mass exponent plot, multifractal spectrum were developed for all the six time series of each station and they are presented Figure 3. From these plots, the important multifractal parameters (H, W, AI, ∆h(q), ∆f(α) and α 0 ) were estimated and the values were provided in Table 2. The value of Hurst exponent is well above 0.5 for all the time series considered in this study, which indicates strong long-term persistence of the datasets. Also, all the 30 time series considered in the study exhibited multifractal behavior. The temperature time series of all the stations consistently displayed high value of H (>0.75) irrespective of the location and altitude, whereas the value of H of air pressure is relatively less (0.58 to 0.64) and that of wind speed are moderate (varying in the ranges 0.61-0.73 and 0.62-0.7 respectively). Higher persistence is noted in the ET0 and T series indicating higher predictability in all the stations, which is quite anticipated. From the time series of U of different stations the complexity it is well evident indicating lesser predictability of the wind speed and air pressure series. The persistence and hence the predictability of RH is less (H ranges between 0.68 and 0.71) at the three low altitude stations. The positive value of asymmetry parameter of different time series indicates fairly strong weighted fractal exponents that are typical in fine structured series. It is also noted that in general, the persistence of SR time series increases with increase in altitude of the stations 0.6 for San Diego to 0.70 for Bakersfield station, which is quite anticipated because of high reception of SR at higher altitudes. The highly complex wind speed series of SM possesses also the high degree of multifractality. The ∆h(q) values show similar trend as that of spectral width indicating the extend of strength of multifractality. The difference between maximum and minimum values of singularity ∆f(α) provides an estimate of the spread in changes in fractal patterns, which denotes the frequency ratio of the largest to the smallest fluctuations. A positive ∆f(α) means that the largest fluctuations are more frequent than smallest fluctuations and vice versa. The value of α0 is indicating the complexity of the time series, and delivers valuable information about the structure of the studied process, with a high value indicating that it is less correlated and possesses fine structure. From   The value of Hurst exponent is well above 0.5 for all the time series considered in this study, which indicates strong long-term persistence of the datasets. Also, all the 30 time series considered in the study exhibited multifractal behavior. The temperature time series of all the stations consistently displayed high value of H (>0.75) irrespective of the location and altitude, whereas the value of H of air pressure is relatively less (0.58 to 0.64) and that of wind speed are moderate (varying in the ranges 0.61-0.73 and 0.62-0.7 respectively). Higher persistence is noted in the ET0 and T series indicating higher predictability in all the stations, which is quite anticipated. From the time series of U of different stations the complexity it is well evident indicating lesser predictability of the wind speed and air pressure series. The persistence and hence the predictability of RH is less (H ranges between 0.68 and 0.71) at the three low altitude stations. The positive value of asymmetry parameter of different time series indicates fairly strong weighted fractal exponents that are typical in fine structured series. It is also noted that in general, the persistence of SR time series increases with increase in altitude of the stations 0.6 for San Diego to 0.70 for Bakersfield station, which is quite anticipated because of high reception of SR at higher altitudes. The highly complex wind speed series of SM possesses also the high degree of multifractality. The ∆h(q) values show similar trend as that of spectral width indicating the extend of strength of multifractality. The difference between maximum and minimum values of singularity ∆f(α) provides an estimate of the spread in changes in fractal patterns, which denotes the frequency ratio of the largest to the smallest fluctuations. A positive ∆f(α) means that the largest fluctuations are more frequent than smallest fluctuations and vice versa. The value of α 0 is indicating the complexity of the time series, and delivers valuable information about the structure of the studied process, with a high value indicating that it is less correlated and possesses fine structure. From Table 2, it is noted that there exists a strong association between α 0 and value of H. On considering the values of these exponents of different stations, the correlation values are found to be 0.966, 0.912, 0.986, 0.99, 0.966 and 0.992 respectively for ET0, T, U, SR, RH and P time series, indicating a strong association between the two exponents in all the time series. This is strongly in agreement with the observation reported by Burgeano et al. [34] for daily temperature series of Catalonia, Spain.
From the non-linear behavior of GHE plots of different time series of different variables (Figure 4), the multifractality of the variables can be confirmed. The scaling exponents of the two high altitude stations such as Dagget and Bakersfield are consistently higher than that of other stations for all moment orders for ET0 and SR time series, whereas a contrasting behavior is noted for the air pressure series. For RH series also, this behavior is noted for most of the moment orders, while there is no consistent pattern for the U time series. In general, the GHE plots of San Diego, Los Angeles and Santa Maria are quite similar, while that of Dagget and Bakersfield are alike. The wind speed time series are highly irregular and accordingly, the multifractal exponents display distinctly different pattern when compared with the rest of the variables. The GHE plots of pressure of time series from all stations except Dagget are similar in magnitude and its multifractal spectra are more symmetrical for all cases.
Atmosphere 2020, 11, x FOR PEER REVIEW 12 of 24  The behavior in GHE is also reflected in the spectra of different time series, while there is no major differences exist in mass exponents of this variable between time series from different stations. Further, the Fourier spectrum of different agro-meteorological time series are developed and the results of Dagget station are presented in Figure 5. The spectral density plots for different variables for remaining stations are provided in the Supplementary Materials, Figures S1-S4.
The behavior in GHE is also reflected in the spectra of different time series, while there is no major differences exist in mass exponents of this variable between time series from different stations. Further, the Fourier spectrum of different agro-meteorological time series are developed and the results of Dagget station are presented in Figure 5. The spectral density plots for different variables for remaining stations are provided in the Supplementary Materials, Figures S1-S4. Fourier spectral analysis estimates clearly detected the scaling from few days to inter-seasonal time scales of up to five months in different time series. The scale break was noted only in the air pressure time series of different time series. One of the scale range is only within a short period of few days within a month, such time scales are not directly captured in MFDFA, as smin is fixed at weekly scale of 8 days considering the stability. The Hurst exponent estimate from weekly to inter seasonal range showed a larger persistence in temperature and relative humidity; medium persistence (0.7-0.8) in ET0 and air pressure, and low persistence in wind speed time series.

Origin of Multifractality
Investigating the origin for multifractality is an important step in multifractal analysis of timeseries. In general, multifractality may result from: (i) long-term temporal correlations and (ii) the broadness of probability density function (PDF). The 'shuffling' procedure and the use of surrogate data are two commonly followed methods for the detection of the source of multifractal behaviour of time-series [82]. The shuffling operation removes any temporal correlations in the time-series, by preserving the probability distribution. The use of a surrogate time-series (generated from the original one) may help to estimate the effect of broadness of PDF. In the first step of shuffling procedure, firstly (m, n) pairs of random integer numbers which satisfies (m, n) < = N, were generated, in which N is the length of the time-series, subsequently the order of entries were interchanged. The process was repeated for large number of times (20N) to ensure that the ordering of entries in the time-series was fully shuffled such that the long or short-term memories were completely destroyed [83]. The shuffling was repeated with ten different initial random realizations to avoid systematic errors in the random number generators. The randomization of the phases of the data in the Fourier space may help in generating the surrogate series [84]. In this study, the most popular Iterative Amplitude Adjusted Fourier Transform (IAAFT) propounded by Schreiber and Schmitz [85] was used for generating surrogate series. The IAAFT method involves the following steps (modified from [86]): (i) arrange the time-series in ranked order and determining the amplitudes; (ii) determine the Fourier transform of initial surrogate (which is a random shuffle of the time-series) and replacing the  Fourier spectral analysis estimates clearly detected the scaling from few days to inter-seasonal time scales of up to five months in different time series. The scale break was noted only in the air pressure time series of different time series. One of the scale range is only within a short period of few days within a month, such time scales are not directly captured in MFDFA, as s min is fixed at weekly scale of 8 days considering the stability. The Hurst exponent estimate from weekly to inter seasonal range showed a larger persistence in temperature and relative humidity; medium persistence (0.7-0.8) in ET0 and air pressure, and low persistence in wind speed time series.

Origin of Multifractality
Investigating the origin for multifractality is an important step in multifractal analysis of time-series. In general, multifractality may result from: (i) long-term temporal correlations and (ii) the broadness of probability density function (PDF). The 'shuffling' procedure and the use of surrogate data are two commonly followed methods for the detection of the source of multifractal behaviour of time-series [82]. The shuffling operation removes any temporal correlations in the time-series, by preserving the probability distribution. The use of a surrogate time-series (generated from the original one) may help to estimate the effect of broadness of PDF. In the first step of shuffling procedure, firstly (m, n) pairs of random integer numbers which satisfies (m, n) < = N, were generated, in which N is the length of the time-series, subsequently the order of entries were interchanged. The process was repeated for large number of times (20N) to ensure that the ordering of entries in the time-series was fully shuffled such that the long or short-term memories were completely destroyed [83]. The shuffling was repeated with ten different initial random realizations to avoid systematic errors in the random number generators. The randomization of the phases of the data in the Fourier space may help in generating the surrogate series [84]. In this study, the most popular Iterative Amplitude Adjusted Fourier Transform (IAAFT) propounded by Schreiber and Schmitz [85] was used for generating surrogate series. The IAAFT method involves the following steps (modified from [86]): (i) arrange the time-series in ranked order and determining the amplitudes; (ii) determine the Fourier transform of initial surrogate (which is a random shuffle of the time-series) and replacing the resulting amplitudes by retaining the phases; (iii) find the inverse Fourier transform and rank them; (iv) replace the amplitudes with the value from the original data with the same position in the ranked list. Hence, the IAAFT algorithm has the advantage over standard phase randomization that the values for the original data are also preserved. The method was also implemented with different initial random seeds. If the multifractality results from temporal correlations, the GHE values of shuffled data for different moment orders are expected to get converged to 0.5. If broadness of PDF is the reason for multifractality, GHE values of surrogate series will be independent of q [82]. If both factors are responsible for multifractal behavior, both the shuffled and surrogated series will show lower multifractal properties than the original series. The GHE plots of original, shuffled and surrogate series can be compared for getting preliminary information on the cause of multifractality.
The GHE plots were prepared for original, shuffled and surrogate series of all the six variables of all the five stations are estimated. The typical plots of Dagget station are given in Figure 6, whereas plots for remaining stations are included in the Supplementary Materials, Figures S5-S8. From Figure 6 it is noticed that, in general, for all variables the GHE plots of shuffled and surrogate series fall below the GHE plot of original series. This indicates the chances that multifractality can be due to both broadness of PDF and temporal correlations. This is clearly evident in ET0 series for all the moment orders, which was also noted the high altitude Bakersfield station ( Figure S5  resulting amplitudes by retaining the phases; (iii) find the inverse Fourier transform and rank them; (iv) replace the amplitudes with the value from the original data with the same position in the ranked list. Hence, the IAAFT algorithm has the advantage over standard phase randomization that the values for the original data are also preserved. The method was also implemented with different initial random seeds. If the multifractality results from temporal correlations, the GHE values of shuffled data for different moment orders are expected to get converged to 0.5. If broadness of PDF is the reason for multifractality, GHE values of surrogate series will be independent of q [82]. If both factors are responsible for multifractal behavior, both the shuffled and surrogated series will show lower multifractal properties than the original series. The GHE plots of original, shuffled and surrogate series can be compared for getting preliminary information on the cause of multifractality. The GHE plots were prepared for original, shuffled and surrogate series of all the six variables of all the five stations are estimated. The typical plots of Dagget station are given in Figure 6, whereas plots for remaining stations are included in the Supplementary Materials, Figures S5-S8. From Figure  6 it is noticed that, in general, for all variables the GHE plots of shuffled and surrogate series fall below the GHE plot of original series. This indicates the chances that multifractality can be due to both broadness of PDF and temporal correlations. This is clearly evident in ET0 series for all the moment orders, which was also noted the high altitude Bakersfield station ( Figure S5   In none of the cases a perfect q-independency of surrogates are noted, but in all cases the GHE values of different moment order dropped down centering about 0.5. Hence it can be concluded that multifractality of different agro-metrological time series are dominantly due to temporal correlations. In none of the cases a perfect q-independency of surrogates are noted, but in all cases the GHE values of different moment order dropped down centering about 0.5. Hence it can be concluded that multifractality of different agro-metrological time series are dominantly due to temporal correlations.

MFCCA of Meteorological Variables with ET0
Multifractal cross correlation analysis was used to investigate the multiscale cross correlation of different meteorological variables. The scaling exponents were computed for individual time series (with themselves) and the joint scaling exponents were determined for each pair-wise combination with ET0. It is to be noted that in this study results of MFCCA are retrieved for the moment order q = 2 for interpretations and the cross-correlations were determined for different time scales to investigate the change in magnitude and sign [65,66]. Here also a polynomial order m = 1 and a scaling range of [8:600] are chosen for the analysis. Different analysis are performed for q = −4 to 4 to avoid possible elevation bias in the results [60,61]. As the seasonal (90 days) and annual (365 days) time scales are more important for practical applications, these two correlations (seasonal ρ XY (90) and annual ρ XY (365)) are calculated separately to investigate the change in magnitude and sign along with the overall Pearson correlation (ρ o ). The graphical representations such as GHE plot, multifractal spectra and cross correlation plots for each link of different stations were prepared. We presented the plots for Dagget and San Diego (highest and lowest altitude stations) as representatives in the Figures 7 and 8, whereas MFCCA plots for other stations are presented in the Supplementary Materials, Figures S9-S11. Detailed examination of different graphical representations of multifractal exponents is made to get an insight into the spatial diversity associations of different variables with ET0 in a multifractal perspective.
The different GHE plots show that the multifractal exponents of air temperature are higher and that of air pressure are lower than the exponents of ET0 irrespective of the station location. Also, the exponents of RH and SR maintain a consistent pattern with respect to the exponents of ET0, for most of the moment orders. In the case of RH time series, the multifractal exponents are less than that of ET0 for all the moment orders more than unity for the three low altitude stations. The multifracral exponents of wind speed are quite similar to that of ET0 for the dataset of LA station, whereas the exponents of U are consistently lower than that of ET0 for the datasets of remaining stations. The similarity in the multifractal exponents of individual time series is reflected in the multifractal spectra also, i.e., if the GHE values, are closer the spectra will be closer and gets dispersed if the GHE plots are distinctly different. But as a universal property, it is noted that the joint multifractal spectra comprise least α values in comparison with the α values of individual spectra. The pattern of correlation plots is strikingly similar in the SR-ET0 link irrespective of the latitude, altitude and location of the station for a time scale up to~4 months.
In general, it can be noted that the pattern of correlations in different linkages maintain a stable pattern up to 4 months approximately and beyond annul scales, more irregularities are introduced in the correlation behavior in the associations. The pattern of T-ET0 link is similar for the high altitude stations, while a different pattern is noticed for low altitude stations. The correlation behavior mains a definite and consistent pattern in the RH-ET0 link and SR-ET0 link at all the stations, at least up to intra-seasonal scales. The behavior of U-ET0 link is also found to be similar at all the stations except SM, where the association cannot be generalized because of highly complex behavior of the wind speed time series ( Figure S10   Additionally, it was observed that air temperature and solar radiation are positively cross-correlated with ET0, regardless of time scale and location, whereas relative humidity is negatively (anti) cross-correlated with ET0. For air pressure situation is more complicated, as the character of cross correlation with ET0 depends both on time scale and location. For the stations with the highest altitudes (Dagget and Bakersfield) P-ET0 is negatively cross-correlated at most of the time scales, while for stations closer to the coastline the sign of cross-correlation changes from positive to negative at larger time scales beyond the intra-seasonal range. There exists a weak positive power law correlation in the U-ET0 links of all the stations, except the two altitude stations of Los Angeles and San Diego. In the U-ET0 linkages of these two stations, a transition in cross-correlation behavior from positive to negative at larger time scales is noticed.
A more important point is that MFCCA could successfully capture such transitions in the nature of power law cross correlation in different agro-meteorological linkages. Such scale specific transitions in the nature of associations could also be explored using other potential methods like wavelet coherence or Time Dependent Intrinsic Correlation [87,88]. One cannot ignore the influence of local climatological drivers and interplay between different meteorological variables resulted in such dynamics in the nature of association of different variables with ET0. The annual and seasonal power law correlations along with the overall correlations in different linkages are presented in Table 3. The overall correlations are consistently positive in the T-ET0 and SR-ET0 links, while it is strongly negative in RH-ET0 link, irrespective of the location of the station. There exists a moderate positive correlation in the SR-ET0 link at all stations (0.48-0.74 at different locations), while in the U-ET0 link such correlation is noted only at the two high altitude stations of Dagget and Bakersfield. In the P-ET0 and U-ET0 the overall correlation is very weak, which could also be attributed to the alterations in the nature of correlations with the time scale. i.e., in positive correlations in some of the time scales may get nullified by the negative correlation in some other scale and vice versa. Table 3 shows that strength of the association also could also be time scale dependent. In majority of the linkages, the annual correlations are greater than seasonal correlations. But it is to be noted that, against the general notion of strong associations at annual time scale, a stronger correlations can also be deciphered at seasonal time scales in some cases (see, in the linkages of T and SR with ET0 of stations of Bakersfield and Santa Maria stations). The spectral properties such as asymmetry and spectral width are also calculated for each pair wise analysis of different stations, and all the results are summarized in Table 3. From Table 3, it is noticed that the joint scaling exponent is nearly half of the scaling exponent of the individual series. This is universal property in multifractal cross correlation studies [89], which is found to be valid in different pairs of linkages in different stations of California. Apart from this, the joint multifractal spectral width is found to be lower that of the width of individual spectra in all cases. This can be considered as another universal property in multifractal cross correlation studies. Also, it is noted that despite the similar character (positive or negative) of AI of individual spectra, the joint spectra may be of opposing character. This preliminary study is performed to investigate the differences in the multifractal properties and multifractal correlations considering stations primarily based on altitude of the stations. The location of the different stations shows that the three low altitude stations are located at the coastal belts of California. The oceanic proximity plays a dominant role in the climatology of the country in general, and the stations in particular. The large scale atmospheric circulations from oceans line El Nino and La Nina play a major role on the climatology in the stations like San Diego in particular. Hence the role of these circulations on the multifractal properties of the time series also cannot be ignored. Apart from the global parameters like terrestrial radiations and temperature, the local factors like the latitude, topography and local processes, moisture, vegetation etc. along with oceanic and atmospheric circulations may influence the multifractality of the rainfall series. However it will be hard to find a universal pattern such as how the exponents change with latitude, altitude or distance from the coast and attributing the physical reason for multifractality to single indicator may not be justifiable. The multiscale power law correlations using MFCCA is first of this kind in agro-meteorological studies. The analysis shows that MFCCA is a robust alternative to investigate the pair wise associations of agro-meteorological datasets in different time scale. The study could successfully capture the dynamics (i.e., the transitions from positive to negative or vice versa) of associations over the time scale. Even though a preliminary investigation is made on differences in the multifractal properties based on altitude and location of the stations, an in-depth investigation on local meteorological processes will be complementary to this understanding in physical context. Also the uncertainty analysis in the quantification of scaling exponents is yet another open research problem in the domain. Moreover, stemming from this pre-requite analysis and parameters derived from it, the multifractal modeling of different variables need to be performed and its potential applicability against the classical time series or and data driven approaches are to be analyzed.

Conclusions
This study investigated the multifractality of six agro-meteorological variables from California, USA using MFDFA. The analysis showed that all agro-meteorological time series possess strong long term persistence and multifractality. The persistence of T, ET0 and SR are high at all stations, indicating better predictability of these variables. The Hurst exponent (H) of air pressure and wind speed are moderate (varying between 0.61 and 0.75), irrespective of the location or altitude, whereas the H-value of SR and ET0 increases with increase in altitude of the station. MFCCA is found to be a robust method for investigating the power law correlations in agro-meteorological datasets. MFCCA analysis showed that for all the time series considered, the joint scaling exponent is roughly the average of individual scaling exponents, validating the universal property of multifractal cross correlation studies for agro-meteorological datasets. The multifractal exponents of air pressure and wind speed consistently lower than the exponents of ET0, irrespective of station location. The α-values of joint multifractal spectrum are lower than the α values of both of the individual spectra, indicating a universality in the multifractal cross correlation studies. Additionally, MFCCA method allows assessing the impact of time scale on cross-correlation strength. It successfully captured the dynamics in the nature of associations of P and U time series with ET0 over the time scale. The pattern of correlation plots is strikingly similar in the SR-ET0 link irrespective of the latitude, altitude and location of the station. It was found that in general, the temporal evolution of cross-correlation display diversity in different pair-wise associations involving ET0 for low altitude stations and high altitude stations. Multifractal modeling of agro-meteorological datasets, uncertainty analysis of derived multifractal parameters and in-depth investigation on local meteorological processes in conjunction with the MFCCA are some possible avenues for future research.