Multiscale coherence analysis of reference evapotranspiration of north-western Iran using wavelet transform

This study applies different wavelet coherence formulations for investigating the multiscale associations of reference Evapotranspiration (ET0) of Tabriz and Urmia stations in North West Iran with five climatic variables, mean temperature (T), pressure (P), relative humidity (RH), wind speed (U) and Solar Radiation (SR). The relationships between different variables are quantified using the Average Wavelet Coherence (AWC) and the Percentage of Significant Coherence (PoSC). The Bivariate Wavelet Coherence (BWC) analysis showed that mean temperature (AWC1⁄4 0.73, PoSC1⁄4 59.18%) and wind speed (AWC1⁄4 0.63, PoSC1⁄4 49.55%) are the dominant predictors at Tabriz and Urmia stations. On considering the Multiple Wavelet Coherence (MWC) analysis, it is noticed that among the two-factor combinations, the T-P and P-RH combinations resulted in the highest coherence values for Tabriz and Urmia stations. T-U-SR combination produced the highest multiple wavelet coherence values among the three-factor cases for both the stations. The Partial Wavelet Coherence (PWC) analysis indicated a drastic reduction in coherence from the values of respective BWC analysis, indicating a strong interrelationship between different variables and ET0. The interrelationship between meteorological variables and ET0 is more apparent at Tabriz, while it is controlled more by the local-scale meteorology at Urmia.


INTRODUCTION
Understanding the association of meteorological variables with reference evapotranspiration (ET 0 ) in diverse process scales is of immense importance in the agricultural water management (Ding et al. 2013). Evapotranspiration is one of the most important components of the hydrologic cycle. ET 0 is the estimation of the evapotranspiration from a 'reference surface'. The reference surface is a hypothetical grass reference crop with an assumed crop height of 0.12 m, a fixed surface resistance of 70 s/m and an albedo of 0.23 (Allen et al. 1998). ET 0 is influenced by the coupled effects of multiple meteorological variables like temperature, pressure, humidity, solar radiation and wind speed. Numerous models ranging from physical to statistical and hybrid artificial intelligence (AI) models were proposed in the past for the accurate prediction of ET 0 (Tabari et al. 2013;Nourani et al. 2019aNourani et al. , 2020. Improved understanding of the relationships among the predictor variables is essential for accurate prediction of ET 0 . The role of different variables in the process of ET 0 and mutual association between the causal variables is of prime importance in developing AI models for the prediction of ET 0 (Nourani et al. 2019). This will help in the selection of the most appropriate set of predictor variables for the estimation of ET 0 . The processes are well understood to be multiscaling in character and an advanced spectral analysis method only can give insight into the processes in multiple time scales (Ding et al. 2013). The influence of different predictor variables on ET 0 may vary at different process scales. Apart from the notion of strong link of different meteorological variables with ET 0 , the interrelationships among these variables, the scale specific association and local-scale processes need to be accounted in modelling. Such observations may help in developing hybrid decomposition modelling of ET 0 with better insights in selection of input variables, accounting for their role in the process of evapotranspiration Araghi et al. 2020;Chia et al. 2020). Also, the identification of predictor contributing the most to the process is very important for its estimation using empirical equations under the datascarce environment with the reasonable degree of accuracy (da Valle Júnior et al. 2021). ET 0 is one of the most important hydrometeorological variables which control the hydrological balance in irrigation water management. Any changes in ET 0 in the agricultural-based region can have severe implications in irrigation demand, land suitability and crop productivity ( Jerin et al. 2021).
Continuous wavelet transform (CWT) and bivariate wavelet coherence (BWC) analysis are popularly used for capturing the dominant periodicity and multiscale association between of time series (Torrence & Compo 1998;Grinsted et al. 2014;Joshi et al. 2016;Rehmati et al. 2020;Jerin et al. 2021). BWC analysis is widely used to study the multiscale associations between two variables. But it is to be noted that most of the hydrological processes are influenced by multiple contributory meteorological variables and there will be strong interactions between them, which makes the governing mechanism too complex. Therefore to unravel such processes completely, the concurrent influence of the different contributory variables in different time scales need to be investigated. To facilitate this, the BWC has been extended to study the relations among multiple variables using the multiple wavelet coherence (MWC) method propounded by Hu & Si (2016). This method has been used in many recent works to study the relationship between variables in a diverse set of problems ranging from geophysical studies to the recent pandemic (Nalley et al. 2019;Su et al. 2019;Das et al. 2020;Song et al. 2020;Islam et al. 2021;Rahman et al. 2021). The variation in the response variable cannot be predicted with an increase in the number of predictor variables because of the cross-correlation among them (Hu & Si 2016). Because of the existing interplay among the contributory variables, a specific bivariate relationship can only be explained clearly by untangling the role of other contributory variables. This can be facilitated by performing partial wavelet coherence (PWC) analysis by excluding the effect of other predictor variables. The PWC method was suggested by Mihanovic et al. (2009) and has been used in several studies (Rathinasamy et al. 2019;Song et al. 2020;Islam et al. 2021;Rahman et al. 2021). Most of the studies in the past have considered one excluding variable at a time while studying the effect of a response variable with a predictor variable (Ng & Chan 2012). In short, the BWC analysis alone cannot capture the coherence relationships between response and contributory variables in hydrological processes, but instead multiple and their partial effects should also be investigated. Hence, using the MWC and PWC analyses along with simple BWC may help in the comprehensive multiscale coherence study of hydrological variables.
The bivariate coherence relationships of hydrological variables of precipitation in Iran have been studied extensively in the past. Araghi et al. (2017) investigated the influence of three large-scale climate oscillations, namely Arctic Oscillation (AO), North Atlantic Oscillation (NAO) and Southern Oscillation Index (SOI), on precipitation of 1960-2014 periods at 30 locations in Iran. Those researchers have reported SOI as the most influential climatic oscillation on the Iranian precipitation. The major effective period of NAO was found to be greater than 64months while that of AO was found to be less than 32 months. The period of SOI was found to be less than 64months in most parts of the country, except in the north-western part. Ebrahimi et al. (2021) examined the relationships between winter precipitation of Iran with sea surface temperature (SST) anomaly of South Pacific, North Atlantic and Indian oceans using the CWT, cross-wavelet transform (XWT) and the BWC analysis. The monthly gridded oceanic SST with a 2.5°Â2.5°of 1984-2019 was used for the analysis. The results of XWT showed that the 8-16-month period is the predominant period in most of the precipitation zones. There was no significant causal relationship between the Indian Ocean to the precipitation zones. Nourani et al. (2019b) used BWC analysis to investigate the impacts of hydroclimatological variables on the water level fluctuations of Urmia lake with rainfall, runoff, temperature, relative humidity and evaporation. They reported that runoff has the highest coherence (0.9-1) with the lake water level. Nourani et al. (2021) used BWC analysis to investigate the associations between water loss (WL), precipitation (P), discharge (Q), temperature (T) and SST along with the satellite images of the Normalized Difference Vegetation Index (NDVI) from 1990 to 2019. The study reported that the impact of T and SST on water level fluctuations was not perceptible. The impact of discharge was found to be more intense in revealing the dominant role of anthropogenic impacts. Rezaei & Gurdak (2020) investigated interannual to multidecadal climate variability on a diverse set of variables of Lake Urmia water level using singular spectrum analysis (SSA) and BWC analysis. They have considered the impact of El Niño Southern Oscillation (ENSO), North Atlantic Oscillation (NAO), Pacific Decadal Oscillation (PDO) and Atlantic Multidecadal Oscillation (AMO) on precipitation, temperature, lake level, groundwater fluctuations, soil moisture, vegetation coverage and insolation clearness index in the Lake Urmia watershed. Overall, ENSO and PDO have a stronger influence than NAO and AMO on the variability in the water level of Lake Urmia as well as on other hydroclimate variables, except for temperature.
From the detailed review of literature, it is noticed that (i) many studies applied wavelet transform and wavelet coherence for the analysis of hydrometeorological data and most of them investigated the teleconnections of hydrological variables with the large-scale climatic oscillations. None of them investigated the coherence relationships between local meteorological factors with the target variable ET 0 ; (ii) many of the studies investigated the multiscale teleconnections of hydrological variables using simple bivariate coherence alone. Even though many studies considered the effect of partial and concurrent effect of causal variables in the recent past (Rahman et al. 2021), such studies on the ET 0 processes have never been explored. Moreover, the declination of the water level of Lake Urmia is a matter of major concern in the region in north-western Iran (Schulz et al. 2020;Fathian & Vaheddoost 2021). Hence, the studies associated with Lake Urmia are of great significance in the water management of NW Iran. Addressing the above research gaps, the present study investigates the coherence relationships of monthly ET 0 of Tabriz and Urmia stations of North-West Iran, with five predictor variables, namely mean temperature (T), surface pressure (P), relative humidity (RH), wind speed (U) and Solar Radiation (SR) using different variants of BWC methods. The specific objectives of the study are framed as follows: (i) to assess the individual influence of various predictor variables on ET 0 of Tabriz and Urmia stations and identify the dominant predictor variable; (ii) to quantify the joint influence of various predictor variables on ET 0 and identify the most influential combination of climatic variables on the process and (iii) to estimate the effect of an excluding variable in predicting the values of ET 0 of the two locations of NW Iran. The study area, data details and methodology are presented in the next section. In the subsequent section, the results of different types of wavelet coherence analysis are presented along with relevant discussion. The final section presents the important conclusions from the study. The studies on the coherence relationships ET 0 may helpful for the improved forecasting and eventually for planning and management of irrigation and water resources and ecological monitoring. The outcomes of this study will eventually help for determining crop water requirements and proper management of water resources for irrigation, irrigation scheduling, designing irrigation systems, mitigating the climate change impacts, etc.

DATA AND METHODS
In this section, the details of study area, data used and the methodology followed are presented.

Study area
In this study, the data of two stations namely Tabriz and Urmia in the north-west Iran, located to the banks of Lake Urmia are considered. Lake Urmia is the second largest salty lake of the world facing dramatic problem of to be dried and the researchers in the region are investigating hydroclimatic and anthropogenic factors leading to this issue (Schulz et al. 2020;Nourani et al. 2021). One task in this regard is which province (left or right) has more contribution in the problem (in both natural and anthropogenic points of view). So any additional insights gained from the studies related to the region are believed to be beneficial in the improved modelling of ET 0 and subsequent irrigation water management. The locations of the meteorological stations are marked in Figure 1. Tabriz city having an average altitude of 1,350 m is located in north-western Iran (Lat. 38°8 0 N -Long. 46°15 0 E), after the confluence of river Qurichai and river Ajichai. The annual precipitation of Tabriz is nearly 380 mm, which experiences dry summer. A semi hot climate in the summer season and mild climate in the spring season prevails in the area. The mean annual temperature is ∼13°C and annual potential evapotranspiration is 1,733 mm. Urmia (Lat. 37°34 0 N -Long. 44°58 0 E and altitude ∼1,300 m) is a city located in north-western Iran. The city experiences very little rainfall in the summer season with an average annual rainfall of 300 mm (Nourani et al. 2019a). In most of the years, freezing happens for about 120 days. High altitude locations receive solar irradiance at a faster rate, which results in the increase in air temperature and finally results in the increased rate of evapotranspiration. This is reflected in the maximum temperature and daily pan evaporation of Tabriz station (35.2°C and 15.3 mm) against the values of 33.9°C and 11 mm for Urmia station. However, the Urmia Lake increases the evaporation from the water bodies and as a result, Urmia receives more maximum and mean precipitation at daily scale (3.73 and 0.66 mm) than Tabriz (2.9 and 0.59 mm).

Data
In this study, the meteorological variables such as temperature (T) (°C), surface pressure (P) (hPa), relative humidity (RH) (%), wind speed (U) at 2 m height (ms À1 ), solar radiation (SR) (J cm À2 day À1 ) and reference evapotranspiration (ET 0 ) (mm) at daily time scale, collected from meteorological stations of Tabriz and Urmia cities. Time-series plots of data (1992-2016 and 1992-2018 periods) used in the study are provided in Figures 2 and 3, respectively, for Tabriz and Urmia stations. Statistical properties of time series of hydrometeorological variables of Tabriz and Urmia stations are provided in Table 1.
The collected meteorological data should be checked for its correctness and continuity before any analysis and standard procedures must be followed to ensure the data quality (Feng et al. 2004;Islam et al. 2020). The internal consistency test, step test, persistence test and range (fixed or dynamic) test are commonly followed for quality control of the datasets (Estévez et al. 2011(Estévez et al. , 2016Herrera-Grimaldi et al. 2019). Apart from the visual examination to detection of unacceptable values like negativity or outliers, the value of RH exceeding 100, etc., the range test procedure was applied prior to the analysis, because of its ability to detect erroneous data (data outside acceptable fixed range). More details on this test and ranges are available in Nourani et al. (2019a).

Estimation of ET 0
The FAO Penman-Monteith (PM) method is used to determine ET 0 . The PM equation derived by combining the energy balance and mass transfer methods to compute the evaporation from an open water surface using standard meteorological parameters. The factors like sunshine, temperature, humidity and wind speed were extended to cropped surfaces by incorporating aerodynamic and bulk resistance factors (Allen et al. 1998).
ET 0 in mm day À1 is computed as where R n is the net radiation at the crop surface (MJ m À2 day À1 ), G is the soil heat flux density (MJ m À2 day À1 ), T is the mean daily air temperature at 2 m height (°C), u 2 is wind speed at 2 m height (m s À1 ), e s is the saturation vapour pressure (kPa), e a is the actual vapour pressure (kPa), (e se a ) is the saturation vapour pressure deficit (kPa), Δ is the slope of vapour pressure curve (kPa°C À1 ) and γ is the psychrometric constant (kPa°C À1 ). Subsequently, the monthly mean values are considered for wavelet coherence analysis, as the ET 0 calculated with mean monthly weather data is indeed very similar to the average of the daily ET 0 values calculated with daily average weather data for that month (http://www.fao.org/3/X0490E/x0490e08.htm). Multiscale spectral analysis tools of Continuous Wavelet Transform (CWT), Wavelet Coherence (WC), Multiple Wavelet Coherence (MWC) and Partial Wavelet Coherence (PWC) are used to study the influence of various climatic variables on reference evapotranspiration. The theory of wavelet transforms is well debated in literature (Torrence & Compo 1998;Grinsted et al. 2004). The extension of CWT namely wavelet coherence and its variants can be used for investigating the relationship between different time series. The analytical framework of the procedures is given in Figure 4.

Continuous wavelet transform
The Morlet wavelet is commonly used in the field of hydroclimatology and it is of the mathematical form: where c o defines the Morlet wavelet with ω 0 and η representing dimensionless frequency and time, respectively. Morlet wavelet (with ω 0 ¼ 6) provides a good balance between time and frequency localization. It allows separating the phase and amplitude of any time series and is localized in scale, which helps in achieving high resolution in frequency. The wavelet power spectrum is averaged over time and scale, to make it convenient for interpretations. The time average of power at a given frequency known as global wavelet spectrum (GWS) is an 'efficient' estimator of the 'true' power  Journal of Water and Climate Change Vol 00 No 0, 6 Uncorrected Proof spectrum. GWS was used to test for the presence of significant periodic components in the time series. The peak of the GWS indicates the dominant periodicity.

Cross-wavelet transform
The cross-wavelet spectrum exposes regions with high common power and reveals information about phase relationship. Considering two signals, X (n ¼ 1,…, N) and Y (n ¼ 1,…, N), the XWT can be used to find the cross-power as where W X (s, t) is the continuous wavelet transform of X series and W Y* (s, t) is their complex conjugation. The cross-wavelet spectrum exposes regions with high common power and reveals information about phase relationship.

Bivariate wavelet coherence
The cross-spectrum is normalized with respect to the wavelet energy of X and Y in the case of wavelet coherence. The wavelet coherence evaluates the linkage between two time series within the time-frequency space by measuring the correlation between the time series, whose values vary between 0 and 1. The wavelet coherence between two time series X and Y is given by: where R(X,Y ) is the measure of the wavelet coherence between Y and X; R 2 (X, Y ) is the measure of the squared wavelet coherence between Y and X; W(X, Y ) denotes the corresponding cross-wavelet transform; W(·) denotes the wavelet transform; and j denotes a smoothing operator that can be used to balance between the desired time-frequency resolution and statistical significance. The phase in the wavelet analysis is usually measured as an angle in degrees over a specific point in the waveform cycle and represents the relative displacement among or between two signals (Araghi et al. 2017). When the two signals have no phase difference, then they are called in-phase (anti-phase).

Multiple wavelet coherence
Multiple Wavelet Coherence (MWC) is based on auto-and cross-wavelet power spectra among the analyzed variables. In this study, the response variable (Y ) is ET 0 . The MWC s 2 m (s, t) at scale s and time t is expressed as (Hu & Si 2016): is the matrix of the smoothed cross-wavelet power between the response and predictor variables; W t X,X (s, t) is the matrix of the smoothed auto-and cross-wavelet power among the multiple predictor variables and W t Y,Y (s, t) is the matrix of the smoothed wavelet power of the response variable. The W Y,XÃ t (s, t) represents the complex conjugate of W t Y,X (s, t). The 5% significance level and the confidence intervals of MWC can be calculated using a Monte Carlo simulation (1,000 repetitions) based on the first-order autocorrelation coefficient (Hu & Si 2016).

Partial wavelet coherence
The Partial Wavelet Coherence (PWC) is appropriate for finding the partial correlation between response variable Y and predictor variable X after eliminating the influence of the predictor variables Z. It can be defined as per Ng & Chan (2012) as: This formulation is similar to the simple BWC, for which the coherence value ranges from 0 to 1. In this case, a low PWC squared shown at where a high wavelet coherence squared was found implies that the time series X1 does not have a significant influence on the time series Y at that particular time-frequency space, and time series X2 dominates the effect on the variance of Y, and vice versa for the opposite case. If both σ p 2 (Y, X1, X2) and σ p 2 (Y, X2, X1) still have significant bands, both X1 and X2 have a significant influence on Y. Very recently, Hu & Si (2021) extended this approach by simultaneous accounting for the removal of multiple variables simultaneously and phase differences. But the proposers reported some inconsistency of results in some cases, emphasizing the necessity of more investigations on their proposed formulation.

Significance of WC
The Average Wavelet Coherence (AWC) and the Percentage of Significant Coherence (PoSC) can be used as measurements to examine the relative dominance of individual or combinations of predictor variables (Nalley et al. 2019;Nalley 2020;Song et al. 2020). AWC can be calculated by averaging the wavelet coherence produced overall scales according to the coherence values computed. The PoSC can be obtained by calculating the ratio of the number of significant values of power over the total number of values of the power produced in the WC computation. Significant power occurs when the ratio of the power over the significance level is greater than 1. Higher overall AWC and PoSC values indicate more dominance (Nalley 2020). This implies that the climatic variable(s) involved in MWC computations are needed to explain the variation in the respective ET 0 dataset analyzed. Although an additional variable may increase the coherence with ET 0 , the increase in PoSC by at least 5% should be observed before concluding the additional variable has a practical significance (Nalley 2020).

RESULTS
Firstly, a statistical cross-correlation between different variables and ET 0 is made and the correlation matrix is presented in Table 2. Positive correlation was obtained for temperature, wind speed and solar radiation, while negative correlation was obtained for pressure and relative humidity for both regions. The maximum correlation value was obtained for temperature (0.97 for Tabriz and 0.933 for Urmia). The correlation pattern is slightly different at two stations, even though the nature of association between the variables is identical for the two stations. The correlation values of U with other variables of Urmia are found to be very less (0.25, 0.24, À0.41, À0.36 and 0.129, respectively, with E, T, P, RH and SR) at this station. Also except P-U and RH-U, the correlations are not statistically significant at a 5% significance level. The highly irregular behaviour of U of Urmia station is one of the reasons which results in this weak linear associations with any of the time series. Moreover, ET 0 process at Urmia is modulated significantly by many local factors because of the presence of Lake Urmia. SR is expected to display high correlation with other variables of the high altitude stations. In this case, the correlation of SR with other variables is found to be high at Tabriz station, which is a high altitude station compared with Urmia. According to the sensitivity analysis performed in the preface study of this work (Nourani et al. 2019a), it was found that temperature is the most predominant physical factor in the evapotranspiration process and modelling, followed by solar radiation for both stations. However, wind speed and relative humidity are other important parameters in Tabriz and Urmia cities, respectively. However, it is well known that in a multiscale process like evapotranspiration, the simple Pearson correlation alone cannot decipher the relationships between the variables. At some of the scale/time spell the correlation may be positive, but at some other scales, the correlation may be negative, which eventually nullify the overall correlation. Therefore, in multiscale processes with non-linear relationship between the variables, broad inferences only can be made using simple correlation analysis (Adarsh & Janga Reddy 2018). This warrants a multi time scale investigation using coherence analysis, as one cannot ignore a strong correlation at some of the localized time spell over the time domain. Therefore, first the wavelet transform was applied to determine the dominant periodic scales of different time series. Apart from this the coherence relationship of an individual variable with ET 0 is studied using BWC. The concurrent and standalone roles of multiple variables on ET 0 are analysed with the help of MWC and PWC methods.

Continuous wavelet transform analysis of ET 0 and causal variables
The localized variations in the time series of predictor and response variables are studied using CWT using the Morlet wavelet. The statistically significant wavelet power at 5% significance levels in the time-frequency space is represented by thick contour lines (Grinsted et al. 2004). Figures 5 and 6 show the CWT spectra of different variables of Tabriz and Urmia stations, respectively. The wavelet transform analysis detected annual periodicity (8-16 months) for all the time series, except for the wind speed of Urmia station. Lower periodicities (2, 4 and 6 months) are observed for pressure, relative humidity and solar radiation for both the stations, but they are very localized in time domain, indicating the role of factors of local meteorology in the governing process. The wavelet coefficients obtained from CWT analysis were used to perform the BWC analysis. Here, bivariate relationships are investigated first, considering one governing factor with ET 0 at a time, so-called BWC analysis. This may help to identify the most contributing factor to ET 0 process of the two locations in a broad sense. The results of BWC analysis of different variables of both the stations are presented in Figure 7.
The coherence relationships are dominantly ranging between periodicities of 8 to 16 months in general and extend even up to 32 months in some cases like that of the temperature series of both the stations. Also, significant coherences are noticed in the relationship of different meteorological variables with ET 0 of the respective stations in the intra-annual scales (4-8 months) and seasonal scales (,4 months). However, such relationships are highly intermittent in the case of both stations except for the wind speed of Tabriz station. Here, some intra-annual scales are also dominant in most of the time spells,  in contrary to the highly intermittent behaviour of the U-ET 0 relationship of Urmia station. In general, such local periodic scales are more for the relationships of Tabriz station than that for Urmia station. In the relationships of both stations, there exists in-phase relationship in T-ET 0 , U-ET 0 and SR-ET 0 links, while anti-phase relationship prevails in the P-ET 0 and RH-ET 0 links. This reinforces the universal behavioural property of the ET 0 process governed by the controlling meteorological parameters, irrespective of time scale and time domain. Also, it is found to be in line with the linear association of these variables with ET 0 estimated using the simple Pearson correlation method (Table 2).
In a process of identification of dominant predictor variable governing evapotranspiration, it is essential to quantify the coherence relationships in a statistical sense. AWC and PoSC (at a 5% significance level) are considered as two statistical measures used for this purpose (Nalley et al. 2019;Song et al. 2020). Table 3 depicts the values of AWC and PoSC obtained for the BWC for the datasets of two stations. Table 3 shows that the highest (AWC of 0.73) and lowest coherences are noted, respectively, for ET 0 -T and ET 0 -U relationship of Tabriz station. For all cases of Urmia station, the AWC is ranging from 0.55 to 0.62. For Tabriz, temperature (AWC ¼ 0.73; PoSC ¼ 59.18) is the most important factor influencing the reference evapotranspiration, followed by Solar Radiation (AWC ¼ 0.60; PoSC ¼ 46.27) and Relative Humidity (AWC ¼ 0.60; PoSC ¼ 40.45). The results obtained from the BWC confirm the correlation results obtained for this station. Strong coherence is evident in the case of ET 0 and temperature at annual time scale and for longer time scales (64 months). The highest correlation exhibited by the temperature (R ¼ 0.97) with the ET 0 of Tabriz is well reflected in the BWC spectrum as well as in the AWC value. For Urmia station, wind speed (AWC ¼ 0.63; PoSC ¼ 49.55) is the highest contributing factor for ET 0 , followed by Pressure (AWC ¼ 0.59; PoSC ¼ 43.87) and Relative Humidity (AWC ¼ 0.58; PoSC ¼ 36.57). It is noteworthy to mention that the difference in the coherence strength is marginal (0.63, 0.59 and 0.58) on considering the dominant variables of Urmia than that of Tabriz. It means that the interplay among the variables is more influential at Urmia than that at Tabriz. The reason of dominancy of SR in Tabriz could be attributed to high altitude, windy regime of Tabriz, the low strength of interrelationships between the variables, the lack of external controlling factors (CFs) like the presence of large water body, etc. The dominance of multiple variables including wind speed and humidity in Urmia station can be attributed to the closer distance of Urmia station to the lake. It is to be noted that the overall correlation between U and different variables at Urmia station was found to be very less (Table 2). But in the multiscale and non-linear process of ET 0 , multiple factors and the interrelationships among them influences mechanism and the use of wavelet coherence is well justified. The relatively high coherence of U of Urmia station (despite the absence of annual periodic scale in the series) could strongly be attributed to the role of Lake Urmia in governing the process. The lake system is the source of humidity, contributes more moisture, influences the local meteorology, controls the interrelationship among the other variables and eventually influences the overall process of evapotranspiration of the locality. A better picture on such interrelationships could be unveiled by considering the partial or concurrent role of different variables on evapotranspiration, which in turn can be performed using MWC and PWC approaches.

Wavelet analysis using multiple wavelet coherence
The combined effect of climatic variables of Tabriz and Urmia stations contributing to ET 0 of the respective station is investigated using MWC (Hu & Si 2016). In this exercise, a total of 23 combinations are considered, considering multiple number of factors varying from 2 to 5. The AWC and PoSC values of MWC analysis of different combinations are presented in Table 4. On considering the different two-factor combinations, the combined influence of temperature and pressure resulted in highest coherence values for Tabriz (AWC ¼ 0.89, PoSC ¼ 87.5). For Urmia, the combination of pressure and relative humidity produced high AWC and PoSC values (AWC ¼ 0.83, PoSC ¼ 76.95). Temperature, wind speed and relative humidity combination produced the highest coherence value among the three-factor cases for both of the stations. On increasing the variables to four, the high coherence value is noted in the case of both the stations. The highest coherence is noted for E-T-RH-U-SR combination for both the station (AWC of 0.99 and 0.98, respectively, for Tabriz and Urmia). It may also be noted that higher wavelet coherence can be noted even with the lesser number of influential variables (for example, E-T-P-SR, ET-RH-SR, etc., of Tabriz and E-T-U-SR, E-P-RH-U, etc., of Urmia) and hence simply increasing the number of variables alone will not result in high coherence. That is, the physical role of the included variable is crucial in the governing process, as at some of the process scales the role of the variable may be to modulate the evapotranspiration, while it may dampen it at some other scales. Unlike many of the teleconnection processes of Climatic Oscillations and hydrologic variables like streamflow/rainfall (Nalley et al. 2019;Su et al. 2019), in this study, the correlation of different variables with ET 0 is also very strong. This is also reflected in the coherence relationships by showing strong coherence in the MWC relationships.

Wavelet analysis using partial wavelet coherence
The wavelet coherence between a specific response variable and ET 0 is studied by eliminating the influence of other variables using PWC analysis using the procedure proposed by Ng & Chan (2012). The PWC spectrum obtained for the stations is shown in Figures 8 and 9. A periodicity of 8-32 months is visible in the relationship between temperature and ET 0 by excluding other variables. The significant high-power region is reduced upon eliminating the effect of temperature from other variables in the case of both stations. The AWC and PoSC values are computed for the partial coherence analysis of different cases and their values are presented in Table 5.
The highest changes in AWC are noticed on the removal of SR and RH from the ET 0 -T relationship of Tabriz station (AWC becomes 0.52 and 0.53 in these cases against the value of 0.73 in BWC between T and ET 0 ). On removal of different variables from the ET 0 -T relationship of Urmia station, it is noted that the difference in PoSC value was not significant, as the changes are less than 5% (Nalley et al. 2019). This is also well reflected in the plots of PWC (Figure 9), as most of the significant power remains on comparing with BWC (Figure 7), even on considering the partial association. More local factors due to the presence of Lake Urmia may be governing the process in the ET 0 -T relationship of this station. On considering the partial relationship, only for ET 0 -T-P case (removal of P from ET 0 -T relationship), the coherence values are highest for both the stations, but all the other cases, the diverse set of variables display the highest coherence. In all the cases, the AWC (and PoSC) values of PWC are lower than the corresponding bivariate relationship, indicating that there exists an interrelationship of different meteorological variables in the process of ET 0 . This is also logical observation, as the physical link between the meteorological variables and ET 0 is more perceptible in most of the regions than hydroclimatic teleconnections. In general, it is noted that the change in PoSC and AWC are more for the datasets of Tabriz station that than of Urmia station, indicating a strong interrelationship among the variables of this station. At the location of Urmia, because of the presence of Lake Urmia, the local-scale meteorology plays an important role which scatters the interrelationship among the variables, more This resulted in a drastic drop in coherence (PoSC) to 9.52% from 46% for Tabriz station on removing the effect of T from ET 0 -SR link. Similar observations can also be made in all the remaining PWC plots. In this study, we attempted to recognize the possible influences of the local meteorological variables on ET 0 of two locations of NW Iran using three types of wavelet coherence approaches. It is worth mentioning that we have considered only the local meteorological variables for performing the multiscale coherence analysis of ET 0 . Instead, the role of large-scale climate oscillations along with local meteorological variables can be considered in analysis. The simultaneous removal of more than one variable in the PWC is a possible analysis to be explored and some of the very recent formulations by Hu & Si (2021) offer flexibility to perform such analysis. Also, improved formulations like terrestrial combined terrestrial evapotranspiration index (CTEI) (Dharpure et al. 2020;Elbeltagi et al. 2021) involving both satellite data and hydrometeorological parameters can be used instead of ET 0 . In such cases, the more number of input variables and their roles can be considered which may give an improved understanding of the processes, which in turn may help better irrigation water management of the study area. Uncorrected Proof CONCLUSION Evapotranspiration is one of the most important processes in the global water cycle. This study used CWT-based WC approaches for investigating the scale specific association between reference evapotranspiration (ET 0 ) and its CFs for the first time. The Bivariate, Multivariate and Partial variants of WC are used to investigate the standalone, concurrent and partial influence of five meteorological variables on ET 0 of Tabriz and Urmia in North-West Iran. The CWT analysis performed on the time series of each variable revealed annual periodicity for all variables, except for the wind speed of Urmia station. The mean temperature is identified as the dominant factor for Tabriz by Pearson correlation as well as BWC analysis. The presence of Lake Urmia resulted in a high coherence in the relationship of wind speed with ET 0 , despite the low linear association between the two. The influence of multiple variables studied using MWC showed that among the two-factor combinations, the combined influence of temperature and pressure resulted in the highest coherence values for Tabriz (AWC ¼ 0.89, PoSC ¼ 87.5) station. For Urmia station, the combination of pressure and relative humidity produced the highest AWC and PoSC values (AWC ¼ 0.83, PoSC ¼ 76.95). Temperature, wind speed and solar radiation combination produced the highest coherence among the different three-factor combinations. Simply increasing the number of meteorological factors alone not necessarily brings the higher coherence, instead their role in governing process decide the coherence relationship with ET 0 . The partial correlation analysis between variables performed using PWC detected a visible periodicity of 8-32 months in the relationship between temperature and ET 0 , even on excluding other variables, for the data of both the stations. In a specific CF-ET 0 PWC analysis, the significant high-power region is reduced upon excluding other variables over the BWC analysis, in the data of both the stations, indicating a high degree of mutual association among different CFs. This study found that multiple CFs have a significant influence on the ET 0 of Urmia than that of Tabriz, as the local-scale processes play a crucial role in the interrelationships because of the proximity with the Lake Urmia. Further investigative research accounting for this factor is recommended to completely unravel the evapotranspiration process of Urmia. The present study considered only the local-scale meteorological variables (MVs) for the coherence analysis of ET 0 , whereas incorporating the climatic oscillations along with MVs is a possible scope to improve this investigation. The usage of datasets of shorter duration is another shortcoming of the study. Using the datasets for a longer duration may help to capture influence of changing climate conditions on the ET 0 at any region, as the methodology presented is a generic one. The scale specific information gained from this study may be helpful in improving the modelling efforts of ET 0 and hence the irrigation water management of the selected locations of NW Iran.