Complexity of Forces Driving Trend of Reference Evapotranspiration and Signals of Climate Change

: Understanding the trends of reference evapotranspiration ( ET o ) and its inﬂuential meteorological variables due to climate change is required for studying the hydrological cycle, vegetation restoration, and regional agricultural production. Although several studies have evaluated these trends, they su ﬀ er from a number of drawbacks: (1) they used data series of less than 50 years; (2) they evaluated the individual impact of a few climatic variables on ET o , and thus could not represent the interactive e ﬀ ects of all forces driving trends of ET o ; (3) they mostly studied trends of ET o and meteorological variables in similar climate regions; (4) they often did not eliminate the impact of serial correlations on the trends of ET o and meteorological variables; and ﬁnally (5) they did not study the extremum values of meteorological variables and ET o . This study overcame the abovementioned shortcomings by (1) analyzing the 50-year (1961–2010) annual trends of ET o and 12 meteorological variables from 18 study sites in contrasting climate types in Iran, (2) removing the e ﬀ ect of serial correlations on the trends analysis via the trend-free pre-whitening approach, (3) determining the most important meteorological variables that control the variations of ET o , and (4) evaluating the coincidence of annual extremum values of meteorological variables and ET o . The results showed that ET o and several meteorological variables (namely wind speed, vapor pressure deﬁcit, cloudy days, minimum relative humidity, and mean, maximum and minimum air temperature) had signiﬁcant trends at the conﬁdence level of 95% in more than 50% of the study sites. These signiﬁcant trends were indicative of climate change in many regions of Iran. It was also found that the wind speed ( WS ) had the most signiﬁcant inﬂuence on the trend of ET o in most of the study sites, especially in the years with extremum values of ET o . In 83.3% of the study sites (i.e., all arid, Mediterranean and humid regions and 66.7% of semiarid regions), both ET o and WS reached their extremum values in the same year. The signiﬁcant changes in ET o due to WS and other meteorological variables have made it necessary to optimize cropping patterns in Iran.


Introduction
Assessment of changes in reference evapotranspiration (ET o ) is required in climate change studies, agricultural and forest meteorology, irrigation scheduling, surface water balance, drought analysis, long-term decision-making in food and water security policies, and optimum allocation of water resources [1][2][3][4][5][6][7]. Evaluating the trends of meteorological variables may help determine the effects of major factors on ET o and climate change [8,9].
Dadaser-Celik et al. [10] evaluated the 32-year trend of ET o in Turkey. Analysis of climatic data showed an upward trend in air temperature, and downward trends in wind speed and relative humidity in Turkey. Changes in these three variables explained the majority of variations in ET o . Song et al. [11] assessed the 46-year trend of ET o in the North China Plain. Their results indicated that the downward trends of net radiation and wind speed had a bigger impact on ET o compared to the upward trends of maximum and minimum air temperature. Wang et al. [12] characterized the 31-year trend of ET o in the western Heihe River Basin in China. They found that wind speed and sunshine duration were the two key meteorological variables that decreased ET o .
Darshana et al. [13] investigated the 30-year trend of ET o in the Tons River Basin in central India. Their outcomes showed that maximum air temperature and net radiation had a stronger impact on ET o compared to minimum air temperature and relative humidity. Zongxing et al. [14] studied the 49-year trend of ET o in the south-west of China. The decrease in wind speed was the main driving force for the reduction of ET o . This happened because the lower wind speed raised the vapor pressure, which ultimately reduced the evaporative demand of atmosphere. Li et al. [15] examined the 46-year trend of ET o in the Upper Mekong River Basin. They showed that sunshine duration had a more significant (at the confidence level 95%) effect on ET o compared to air temperature, relative humidity, and wind speed. Gao et al. [16] evaluated the 45-year trend of ET o in China, and found that sunshine duration, wind speed, and relative humidity had a more important influence on ET o than air temperature. Zhang et al. [17] assessed the changes in ET o and its controlling factors in China. They indicated that maximum air temperature, relative humidity, and wind speed affected ET o . A number of studies (e.g., [8,9,[18][19][20]) showed that the increasing and decreasing trends of ET o were mainly due to the increase in air temperature and decrease in wind speed, respectively.
Tabari et al. [21,22] evaluated the 40-year (1966-2005) trend of ET o in the west and southwest of Iran with arid and semiarid climates. They studied the effect of air temperature, relative humidity, vapor pressure, wind speed, and rainfall on ET o , and found that wind speed had the most influence on ET o . Unfortunately, they did not consider the impact of serial correlations (autocorrelation) on the trends of ET o and meteorological variables. Kousari and Ahani [23] assessed trends of ET o and six meteorological variables (mean, minimum and maximum air temperature, relative humidity, wind speed and sunshine hours) in three climate regions (arid, semiarid and humid) of Iran during 1975-2005. They did not take into account the influence of serial correlations on the trends. Shadmani et al. [24] evaluated the 41-year (1965-2005) trend of ET o from 11 weather stations in arid regions of Iran but did not specify which meteorological variables controlled trend of ET o .
The existing studies (1) evaluated the impact of various climatic variables on ET o and (2) showed the signals of climate change in Iran based on upward/downward trends of ET o and meteorological variables [8,9,[20][21][22][23]. However, they suffer from the following shortcomings: (1) they used data series of less than 50 years which are suitable only for evaluating climate variability and not climatic change. At least a 50-year period is necessary to study climate change [25][26][27][28][29][30][31][32]. Borman [26] used a 50-year period to analyze the sensitivity of ET o to climate change in Germany. Kingston et al. [27] found a considerable uncertainty in evaluating the changes of ET o due to climate change using 30-year data. Several studies assessed the response of ET o to climate change in China using 50-year data [28][29][30][31]. Hence, there is a consensus in the literature to use at least a 50-year period to evaluate the trend of ET o due to climate change. ( of serial correlations (autocorrelation) on the trends of ET o and meteorological variables. Finally (5) they did not assess if the extremum values of meteorological variables and ET o occur simultaneously. This study overcame the abovementioned drawbacks by analyzing the 50-year (1961-2010) trends of ET o estimates and 12 meteorological variables in Iran using the trend-free pre-whitening Mann-Kendall (TFPW-MK) and Spearman's Rho tests. These sites were chosen to cover four different climates, namely arid, semiarid, Mediterranean and humid. Moreover, variations of 12 meteorological variables (mean, maximum, and minimum air temperature, difference between the maximum and minimum air temperature, rainfall, wind speed and direction, mean and minimum relative humidity, vapor pressure deficit, number of cloudy days, and sunshine hours) were investigated to characterize forces that drive trends of ET o . Furthermore, the effect of serial correlations on the trends was eliminated by the TFPW approach [24]. Finally, during the period 1961-2010, the coincidence of annual extremum values of ET o and meteorological variables was studied to identify the main meteorological variables that affect variations of ET o .

Data
Daily meteorological data in the 18 study sites were downloaded from the Iran Meteorological Organization (IMO) archive. These data cover a period of 50 years (1961-2010), and include mean, minimum, and maximum daily air temperature ( • C), vapor pressure deficit (kPa), mean and minimum relative humidity (%), wind speed (m/s) at the screen-height of 2 m, rainfall (mm/day), number of cloudy days, and number of sunshine hours per day (hr/day).
Characteristics of the study sites are indicated in Table 1. Figure 1 shows the location of these sites in Iran. They are chosen to cover various climate regions (arid, semiarid, humid, and Mediterranean) across Iran.

FPM Equation
ETo is the rate of evapotranspiration from a uniform height, actively growing, well-watered, and completely shaded hypothetical crop [33]. The hypothetical crop has a height 0.12 (m), a fixed surface resistance of 70 (s/m), and an albedo of 0.23, closely resembling the evapotranspiration from an extensive surface of green grass [26]. In this study, daily meteorological variables from the 18 study sites were used in the Agricultural Organization of the United Nation (FAO)-Penman-Monteith (FPM) equation to estimate ETo on a daily basis. ETo estimates from the FPM equation were validated against the lysimeter measurements in the 18 study sites [34][35][36][37].
The FPM equation is given by [33]: where ETo is the daily reference evapotranspiration (mm/day), γ is the psychrometric constant (kPa/℃), Δ is the slope of the saturation vapor pressure-temperature curve (kPa/℃), T is the mean daily air temperature (℃), and u is the mean daily wind speed at the screen-height of 2 m (m/s). G is the ground heat flux (MJ/m 2 /day) and is negligible on daily timescales [21,33]. es and ea are, respectively, the saturation and actual vapor pressures (kPa), which are given by [33]: where ( ) and ( ) are the saturation vapor pressures (kPa) at daily maximum and minimum air temperatures, respectively. RHmax and RHmin are the daily maximum and minimum relative humidity (%), respectively.

FPM Equation
ET o is the rate of evapotranspiration from a uniform height, actively growing, well-watered, and completely shaded hypothetical crop [33]. The hypothetical crop has a height 0.12 (m), a fixed surface resistance of 70 (s/m), and an albedo of 0.23, closely resembling the evapotranspiration from an extensive surface of green grass [26]. In this study, daily meteorological variables from the 18 study sites were used in the Agricultural Organization of the United Nation (FAO)-Penman-Monteith (FPM) equation to estimate ET o on a daily basis. ET o estimates from the FPM equation were validated against the lysimeter measurements in the 18 study sites [34][35][36][37].
The FPM equation is given by [33]: where ET o is the daily reference evapotranspiration (mm/day), γ is the psychrometric constant (kPa/ • C), ∆ is the slope of the saturation vapor pressure-temperature curve (kPa/ • C), T is the mean daily air temperature ( • C), and u is the mean daily wind speed at the screen-height of 2 m (m/s). G is the ground heat flux (MJ/m 2 /day) and is negligible on daily timescales [21,33]. es and ea are, respectively, the saturation and actual vapor pressures (kPa), which are given by [33]: where e (Tmax) and e (Tmin) are the saturation vapor pressures (kPa) at daily maximum and minimum air temperatures, respectively. RHmax and RHmin are the daily maximum and minimum relative humidity (%), respectively. R n is the net radiation (MJ/m 2 /day), and is estimated by [21,33]: where n is the number of sunshine hours per day (hr/day), n max is the maximum possible duration of sunshine per day (hr/day), and R a is the extraterrestrial radiation (MJ/m 2 /day). R a and n max depend on the latitude and julian day [33]. a s and b s are empirical coefficients, which are obtained from [21] for each study site. In this study, the annual averages of daily FPM ET o estimates and meteorological variables were used to analyze the annual trends.

Mann-Kendall (MK) Test
One of the most well-known methods for detecting trends in a time series is the non-parametric Mann-Kendall (MK) test [38][39][40]. Unlike the parametric statistical approaches, the non-parametric statistical tests are more suitable for non-Gaussian distributed data, which are frequently observed in hydrologic time series [34]. The MK test is defined as follows [38][39][40]: If S > 0 then Z = (S − 1)/VAR (S) 0.5 (8a) where N is the number of data, S is the summation of signs, VAR is the variance, and x j and x i are the data values in years j and i, respectively (with j > i). t p and g indicate ties and the number of ties, respectively. The MK test determines whether to reject the null hypothesis (H0) or accept the alternative hypothesis (Ha), where H0: no monotonic trend is present and Ha: a monotonic trend is present. If 1.65 < Z ≤ 1.96, 1.96 < Z ≤ 2.58, Z > 2.58, the trend is significant at the confidence levels of 90%, 95%, and 99%, respectively [38][39][40]. Local (at-site) significance levels for each trend test can be obtained from Atmosphere 2020, 11, 1081 6 of 26 where | | denotes the absolute value, and Φ ( ) is defined as, where θ is the random variable, and t is the variable for which the cumulative distribution function should be calculated. If p ≤ α, the existing trend is statistically significant at the significance level of α [41,42]. Sometimes the MK test detects trends because of the serial correlations of time series data, which leads to an increased rejection rate of the null hypothesis [41,43].
Hydrologic time series often have significant serial correlations that can undermine the ability of the MK test to correctly assess the significance of trend [41,43].
This method is defined by the following steps: 1. Calculate the slope of trend using the Sen's slope estimator [43,44]: where Q stands for the slope of trend. Q i is the Sen's slope estimator for each value of i, x j and x k are the numerical values at times j and k (j > k), respectively. 2. There is no trend if Q is equal to zero. Otherwise, it is assumed that the existing trend is monotonic, and the time series data are de-trended as follows: where X t is the lag-1 autocorrelation of the de-trended time series and X t is the autocorrelation of time series. 3. Using the rank correlation coefficient estimator, the lag-1 autocorrelation of the de-trended time series is estimated by replacing the sample data by their ranks as follows [45]: where r j is the lag-j autocorrelation coefficient, and X is the average of the data. Then, the estimated lag-1 autocorrelation is removed from the time series as follows: 4. Adding the removed trend in step 2:

Spearman's Rho Test
The Spearman's Rho test was applied to investigate correlation among all the meteorological variables and ET o estimates from the FPM equation. The Spearman's Rho test is defined as [18]: where ρ is the correlation coefficient of linear regression between series i (the order of the data in the original series) and x (data values) [18]. If |Z| > Z α at a significance level of α, then the null hypothesis of no trend is rejected [18]. Following the literature , the confidence levels of 90%, 95%, and 99% were used in this study to evaluate the trends of ET o and meteorological variables.
If there is no mention of confidence level, it meant a confidence level of 95% was adopted. Otherwise, we explicitly mention the confidence levels of 90% and 99% wherever they were used. Table 2 compares p-values from the MK and TFPW-MK tests for ET o and the 12 meteorological variables in the 18 study sites.

Comparison of p-Values from the MK and TFPW-MK Tests
The p-value is a random variable derived from the distribution of the test statistic to analyze a data set and to test a null hypothesis [41,42]. If p-values from the abovementioned tests are different, but they result in an identical trend for a particular time series (i.e., an insignificant trend or a significant trend at the same confidence level), the cells in Table 2 are highlighted in yellow color. If the two tests yield different trends for a specific time series (i.e., one leads to a significant trend, while the other one results in an insignificant trend, or they both lead to a significant trend but with different confidence levels), the cells are highlighted in red color.
As shown, in each study site, p-values from the MK and TFPW-MK tests are different at least for three meteorological variables (yellow and red cells). Overall, 38% of the obtained p-values from the two tests are different due to the serial correlations of time series data. In each study site, the two tests also lead to different trends for at least one variable (red cells). As shown in Table 2, 15% of the trends from two tests are different. These results indicate that the serial correlation undermine the ability of the MK test to correctly determine the trends and their confidence level [41,43]. Hence, the TFPW method should be applied to eliminate the impact of the serial correlations on the MK test.
As can be seen in Table 2, some meteorological variables are subject to different p-values from the MK and TFPW-MK test more often than others. For example, wind direction (WD), wind speed (WS), and rainfall (P) were flagged for more study sites compared to mean air temperature (Tmean) and minimum air temperature (Tmin). This happens because WD, WS, and P time series have more serial correlations compared to Tmean and Tmin time series. Similar findings were reported in other studies [46][47][48][49].

Trends of ET o and Meteorological Variables in the Study Sites
Variations of the meteorological variables were assessed in all of the study sites to identify the ones that control the trend of annual mean of ET o . However, herein, the trends are presented only in two sites, Kerman and Qazvin (Figures 2 and 3).
Atmosphere 2020, 11, x FOR PEER REVIEW 9 of 27 the highest Tmax and n values (during the 50-year period) were observed in 2010, which may an indication of drought in this region [50,[52][53][54][55][56][57]. These results show the complexity of forces driving the trend of ETo and may highlight importance of considering different ETo models and variables in various climates for local and global studies.  in the FPM ETo estimates. It is worth mentioning that these significant trends reduced ETo and may cause alarming climate change in Qazvin [57][58][59]. During the 50-year study period , Tmax reached its highest value in 2010. It should be noted that the minimum es-ea and ETo occurred in 1991 (see the arrows in Figure 3). The highest decreasing rate was for WS (−6.8%/decade). RH showed an upward trend as a result of decreasing WS and es-ea (Figure 3). The warm dry south-eastern winds (called Raz or Shareh) originated from the arid areas of central Iran (Yazd). They were prevailing from 1961 to 2010 and decreased RHmin in Qazvin [60][61][62][63].   Figure 2 shows the 50-year trends in the annual averages of daily FPM ET o estimates, and meteorological variables including mean (Tmean), minimum (Tmin), and maximum (Tmax) air temperature, the difference between maximum and minimum air temperature (Tmax-Tmin), vapor pressure deficit (es-ea), relative humidity (RH), minimum relative humidity (RHmin), rainfall (P), wind speed and direction (WS and WD), number of cloudy days (CD), and sunshine hours (n) in Kerman.

Kerman
As can be seen, the significant upward trends in Tmean, Tmin, Tmax, and n (at confidence level 95%) as well as the substantial downward trends in CD, Tmax-Tmin, P, and WS (at confidence level 95%) led to a negative trend (at confidence level 95%) in the FPM ET o estimates.
Although the upward trends in Tmax, Tmin, Tmean, and n could potentially increase ET o , a downward trend in WS and (es-ea) led to a significant decreasing trend in the FPM ET o values. It is worth mentioning that the minimum of both WS and ET o occurred in the same years (i.e., 1982 and 1984) (see the arrows in Figure 2). Similarly, the course of (es-ea) shows lower values in 1982-1986. Thus, WS and to a lesser extent (es-ea) may be the driving force for the variations of ET o in Kerman.
Although CD showed the strongest downward trend among all the meteorological variables in Kerman, it could not capture the variations of ET o as good as WS. Both WS and ET o showed significant trends at a confidence level of 95%, while CD showed a significant trend at a confidence level of 99%.
During 1960-1980s, the north-western wind was prevailing, which blows from the Zangi Abad region with an arid climate into Kerman [50,51]. However, from 1980 until 2010, the northern wind was dominant, which comes from Chatroud region with dry temperate climate [50,51]. Hence, the change in wind direction (WD) reduced the impact of the arid climate in Kerman, and consequently decreased ET o .
The decreasing rate of Tmax-Tmin (−0.3%/decade) was due to a more rapid increase in Tmin (+3.8%/decade) compared to Tmax (+0.8%/decade). Tmin, P, WS, and CD had the highest changes (larger than ±2% per decade), implying climate change in Kerman [50,51]. The lowest P and CD, and the highest Tmax and n values (during the 50-year period) were observed in 2010, which may an indication of drought in this region [50,[52][53][54][55][56][57]. These results show the complexity of forces driving the trend of ET o and may highlight importance of considering different ET o models and variables in various climates for local and global studies. Figure 3 shows the 50-year trends of annually averaged ET o and meteorological variables in Qazvin. As can be seen in Figure 3, the significant upward trend in RH as well as substantial downward trends in es-ea, WS, CD, and RHmin (at the confidence level 95%) led to a negative trend in the FPM ET o estimates. It is worth mentioning that these significant trends reduced ET o and may cause alarming climate change in Qazvin [57][58][59]. During the 50-year study period (1961-2010), Tmax reached its highest value in 2010. It should be noted that the minimum es-ea and ET o occurred in 1991 (see the arrows in Figure 3). The highest decreasing rate was for WS (−6.8%/decade). RH showed an upward trend as a result of decreasing WS and es-ea ( Figure 3). The warm dry south-eastern winds (called Raz or Shareh) originated from the arid areas of central Iran (Yazd). They were prevailing from 1961 to 2010 and decreased RHmin in Qazvin [60-63]. Table 3 and Figure 4 summarize trends of FPM ET o estimates and meteorological variables in the 18 study sites. The FPM ET o estimates showed significant upward and downward trends, respectively, in 33.3% (28.6% of arid regions vs. 44.4% of semiarid regions) and 22.2% (28.6% of arid regions vs. 22.2% of semiarid regions) of the study sites (Table 3 and Figure 4). Hence, the FPM ET o retrievals had significant variations in most parts of Iran (57.2% of arid regions vs. 66.6% of semiarid regions), which may indicate climatic change signals.  *, **, *** are significant trends at a confidence level of 90%, 95%, and 99%, respectively.

Trends of ET o and Meteorological Variables in Iran
Atmosphere 2020, 11, 1081

of 26
Atmosphere 2020, 11, x FOR PEER REVIEW 13 of 27  Results showed that ET o had a significant positive (negative) trend in the west and east (center) of Iran. WS was the only meteorological variable with the same rising/falling trend. Table 3, in eight of the study sites (i.e., 44.4% of sites), there were significant trends in both ET o and Tmean/RHmin/WS. Moreover, both ET o and Tmin/Tmax/es-ea showed significant trends in 33.3% of the study sites. Hence, these meteorological variables could be considered as driving forces to control variations of ET o in Iran.

As shown in
WS had significant upward and downward trends in 22.2% (14.3% of arid regions vs. 33.3% of semiarid regions) and 33.3% (57.1% of arid regions vs. 22.2% of semiarid regions) of the study sites, respectively. Therefore, significant variations in WS were observed in 55.5% (71.4% of arid regions vs. 55.5% of semiarid regions) of the sites. Table 3 and Figures 2-4 indicate that WS is the driving force for the trends of ET o in 55.5% of the study sites. Thus, WS may be the most important variable that controls the variations of ET o , particularly in the arid regions. These findings are further validated by our outcomes in Table 5 (Section 4.6) that lists the three meteorological variables with the highest correlation with ET o . As shown, WS has the highest significant correlation with ET o in all the study sites expect Rasht (RA). Studying the trends of ET o and meteorological variables over longer time periods and more sites across Iran leads to more reliable results. Table 3 and Figure 4 also showed that there were significant upward trends for Tmean and Tmax in 77.8% (85.7% of arid regions vs. 55.6% of semiarid regions) and 66.7% (100.0% of arid regions vs. 33.3% of semiarid regions) of the study sites, respectively. In Jiroft, Mashhad, and Urmia (3 out of 18 study sites, i.e., 16.7% study sites), there were upward significant trends between ET o and Tmean/Tmax (Table 3). There was a significant downward trend for RHmin in 66.7% (71.4% of arid regions vs. 66.7% of semiarid regions) of the study sites ( Figure 4).
Tmax, Tmean, Tmin, es-ea, and n had significant increasing trends, respectively, in 100.0%, 85.7%, 85.7%, 57.1%, and 57.1% of the arid sites. On the other hand, RHmin, Tmax-Tmin, and WS had significant decreasing trends in 71.4%, 57.1%, and 57.1% of the arid sites, respectively. These results indicated that the arid sites were affected more than the semiarid ones by climate change. In most of the semiarid sites, Tmin and Tmean showed increasing, and RHmin and CD indicated decreasing trends.
As mentioned above, variations of ET o in each study site are controlled by meteorological variables such as air temperature and wind speed. The Arak and Esfahan sites show contrasting trends in ET o although they are relatively close to each other. This happens because wind speed (which is the dominant controlling variable in both sites) has an increasing (decreasing) trend in Arak (Esfahan). Furthermore, Arak and Esfahan have different climate types. Arak is located in a cold mountainous area, while Esfahan is placed in a hot flat desert.
The focus of this study was to evaluate the annual trends of ET o and meteorological variables over a 50-year period (1961-2010). However, according to Table 3 and Figure 4, we can conclude that their seasonality is important. For example, the annual minimum air temperature, Tmin (that occurs in cold seasons, i.e., fall-winter) had significant increasing trends in 83.4% of sites. While, the annual maximum air temperature, Tmax (which happens in warm seasons, i.e., spring-summer) showed significant rising trends in 66.7% of sites. In addition, the significant downward trends of Tmax-Tmin were observed in 41.5% of sites, whereas its upward trends were seen only in 9.3% of sites. Its significant downward trends were due to the higher increasing rate of Tmin in cold seasons than that of Tmax in warm seasons. These results imply that the cold seasons have more influence on the significant trends than warm seasons. It should be noted that these findings are primitive, and future studies should be directed towards evaluating the seasonal variations of ET o and other climatic variables over long periods.  Figure 5 shows the box plots for the 50-year variations of ETo and meteorological variables in each study site.  The second highest variations of ET o was seen in Ahvaz (AH). AH had also the highest variations of Tmean compared to the other study sites. Hence, the variations of ET o in Ahvaz can be related to that of Tmean.

Range of Variations of ETo and Meteorological Variables
The lowest variations in Tmax-Tmin, es-ea, and n values were observed in Rasht (RA), which led to the smallest 50-year change in ET o [9,10]. Most of the outliers occurred in 2010, which signaled a drought condition in different regions of Iran. The FPM ET o values ranged from 2 to 11 mm day −1 in various climates of Iran.
As shown in Equation (1), ET o is affected by the atmospheric vapor pressure deficit (es-ea). As expected, a high positive correlation exists between ET o and es-ea in Figure 5. For example, high values of ET o and es-ea are seen in Zabol (ZA). Moreover, the lowest ET o and es-ea values are observed in Rasht (RA). This happens because the atmospheric demand to water vapor increases as the vapor pressure deficit (es-ea) rises [33]. Positive correlations are also observed between ET o and Tmax, and ET o and n. A larger Tmax leads to more potential for converting soil moisture to water vapor and causes plants to open up their stomata and release more water vapor [33]. Furthermore, larger values of n yield higher R n (Equation (4)) and ET o (Equation (1)). On the other hand, a negative correlation is found between ET o and RH. The highest ET o and lowest RH values are seen in RA. This is due to the fact that the atmospheric demand for water vapor decreases as RH increases. Table 4 shows the years of occurrence of maximum and minimum values of ET o and meteorological variables in each study site during the 50-year study period (i.e., 1961-2010). If the annual extremum values of ET o and a particular meteorological variable coincide, they are highlighted by green color (Table 4).

Annual Extremum Values of ET o and Meteorological Variables during the Study Period (1961-2010)
As can be seen, the maximum of both ET o and WS in Ahvaz, Bushehr, Hamedan, Jiroft, Moghan, Rasht, Sanandaj, Shahrekord, Shiraz, Tabriz, Yazd, and Zabol happened in 1964, 1967, 1967, 2007, 1989, 1975, 1971, 2008, 1981, 1961, 1971, and 1984, respectively. Similarly, the minimum of ET o and WS in Arak, Esfahan, Jiroft, Kerman, Moghan, Sanandaj, Shahrekord, Yazd, and Zabol occurred in 1966, 1996, 1996, 1984, 1994, 1986, 1968, 1995, and 1968, respectively. Overall, in 83.3% (all arid, Mediterranean and humid regions and 66.7% of the semiarid regions) of the study sites, annual extremum values of ET o and WS occurred in the same year. These results imply that WS is the primary driver of the variations in ET o in Iran. In addition, the annual extremum values of ET o and RH coincided in 33.3% of the study sites, namely Ahvaz, Arak, Hamedan, Mashhad, Moghan, and Qazvin. Hence, RH is one of the main driving forces for variations of ET o in these study sites.  Figure 7 shows the correlation coefficient maps among the meteorological variables and reference evapotranspiration for the study sites in arid, Mediterranean, and humid regions. Similarly, Figure 8 indicates the correlation coefficient maps for the sites in semiarid region.

Correlation between ETo and Meteorological Variables
As indicated, WS was the only meteorological variable that had high correlation (ρ ≥ 0.8) with ETo in Ahvaz, Esfahan, Jiroft, Yazd, Zabol, Sanandaj, Arak, Mashhad, Moghan, Qazvin, Shahrekord, Tabriz, and Urmia. These results imply that WS was the driving force for the variations of ETo in most (72.2%) of the study sites in Iran, including 87.5% of arid regions and 77.8% of semiarid regions. This   Figure 7 shows the correlation coefficient maps among the meteorological variables and reference evapotranspiration for the study sites in arid, Mediterranean, and humid regions. Similarly, Figure 8 indicates the correlation coefficient maps for the sites in semiarid region.   Spearman's Rho correlation coefficient (ρ) maps for weather stations located in semiarid region. Colored boxes show a significant correlation at the confidence level of 95%. Dark, medium, and light purples denote 0.8 ≤ ρ ≤ 1, 0.5 ≤ ρ < 0.8, and ρ < 0.5, respectively.

Correlation between ET o and Meteorological Variables
As indicated, WS was the only meteorological variable that had high correlation (ρ ≥ 0.8) with ET o in Ahvaz, Esfahan, Jiroft, Yazd, Zabol, Sanandaj, Arak, Mashhad, Moghan, Qazvin, Shahrekord, Tabriz, and Urmia. These results imply that WS was the driving force for the variations of ET o in most (72.2%) of the study sites in Iran, including 87.5% of arid regions and 77.8% of semiarid regions. This is in agreement with the results obtained by the TFPW-MK test (Table 3 and Figure 4).
WS and WD had the lowest correlation with the other meteorological variables. Hence, these variables may be controlled mainly by human activities such as desertification, land use change, and urban development [67][68][69][70][71][72][73][74][75][76][77][78][79]. The highest correlations were observed among the air temperature-related variables (i.e., Tmean, Tmax, Tmin, Tmax-Tmin) and humidity. In Mashhad and Urmia, the FPM ET o estimates had a good correlation with all of the meteorological variables. The upward trends of Tmin in all the study sites in Iran (except Shahrekord due to its high elevation) indicated a signal of climate change and signaled the need for the optimization of cropping pattern [80,81]. Table 5 shows the total number of meteorological variables in each study site with significant correlations with ET o . It also lists the three meteorological variables that have the highest significant correlation with ET o . As can be seen, ET o had the highest significant correlation with WS in all the sites except Rasht. ET o indicated the largest (the second largest) correlation with Tmax in Rasht (Mahshad, Tabriz, Urmia, and Zabol). According to Table 5, in 66.7% of the study sites, more than five meteorological variables had significant correlations with ET o . This highlights the complexity of the forces driving variations of ET o . In this study, the TFPW-MK and Spearman's Rho tests were used to identify significant (at the confidence levels 90%, 95% and 99%) trends of ET o and meteorological variables across Iran. These significant trends are indicative of climate change signals . For instance, in the Mashhad site, significant upward (downward) trends were seen for ET o , Tmean, Tmin, Tmax, and es-ea (Tmax -Tmin, WD, RH and RHmin), suggesting climate change alarms. These findings are consistent with those of other studies [82][83][84][85][86][87]. In the Shahrekord site, Tmax -Tmin (Tmean, Tmin, RHmin, and CD) showed significant increasing (decreasing) trends, which can be considered as climate change signals. These findings are also in agreement with those of [8,9,[20][21][22][23]. In the Shiraz site, significant upward (downward) trends of ET o , Tmean, Tmin, Tmax, and es-ea (Tmax -Tmin, WS, RH and RHmin) represented signals of climate change [88]. In the Urmia site, growing (reducing) trends of ET o , Tmax, WS, n, Tmean and Tmax (RHmin and CD) signified climate change [89][90][91][92]. Finally, in the Zabol site, upward (downward) trends of ET o , Tmean, Tmin, and WS, and Tmax (RHmin) may be considered as alarms of climate change [84,93].
The FPM ET o retrievals mainly depend on climate variables such as air temperature and humidity, wind speed, and solar radiation. This equation assumes a stomatal resistance of 70 s/m and an albedo of 0.23 for a standard 12 cm grass, and thus it ignores changes in stomatal resistance response to the elevated CO 2 and the surface albedo, ultimately causing uncertainty in the ET o estimates. Surface albedo varies with changes in vegetation and soil moisture [94].
Elevated CO 2 emission influences plant physiology by diminishing stomatal conductance [95][96][97]. Ignoring changes in the surface albedo and vegetation response to the elevated CO 2 emission lead to uncertainty in the FPM ET o estimates and results of this study. A similar uncertainty in the FPM ET o estimates was reported by other studies. For example, Milly and Dunne [98] showed that ignoring the stomatal resistance response to the increased CO 2 emission in the FPM equation led to discrepancy between the ET o estimates from climate models and the FPM equation. In a similar effort, Li et al. [99] indicated that the variations in the FPM ET o estimates are only less than 10% for 40% changes in stomatal resistance. However, there are still high uncertainties in the amount of stomatal response to the raised CO 2 emission. For example, Domec et al. [100] reported that vegetation responses of pine to long term elevated CO 2 were manifested only when soil moisture was at a high level. Future studies should be directed toward evaluating the uncertainties in the FPM ET o estimates with respect to climate change and drought condition [101,102].

Conclusions
This study assessed the trends of meteorological variables (namely, mean, minimum, and maximum air temperature, difference between the maximum and minimum air temperature, mean and minimum relative humidity, wind speed and direction, vapor pressure deficit, rainfall, and sunshine hours) and reference evapotranspiration (ET o ) in the 18 study sites in different climate regions of Iran.
The effect of meteorological variables on ET o was also evaluated. Using the trend-free pre-whitening Mann-Kendall (TFPW-MK) and Spearman's Rho tests, wind speed (WS) was found to be the most important variable that controls the trend of ET o in most regions of Iran. Moreover, the increasing/decreasing trends of other meteorological variables (air temperature and humidity, rainfall, wind speed and direction, and sunshine hours) were indicative of climate change in many regions of Iran. The upward trend of minimum air temperature (Tmin) in all of the study sites in Iran (except Shahrekord due to its high elevation) indicated a signal of climate change and signaled the need for the optimization of cropping system. The significant increasing rate of Tmin may lead to a reduction in the growing degree day (GDD), ultimately decreasing the production of strategic and vital crops in future.
WS had significant upward and downward trends in 22.2% and 33.3% of the study sites, respectively. Therefore, significant variations in WS were observed in 55.5% of the investigated regions. The results also indicated that there were significant upward trends for mean air temperature (Tmean) and maximum air temperature (Tmax) in 77.8% and 66.7% of the study regions, respectively. In most of semiarid climates, Tmin and Tmean showed increasing, and minimum relative humidity (RHmin) and cloudy days (CD) indicated decreasing trends. This indicates that arid regions were affected more than semiarid areas by climate change. Precipitation (P), WS, and CD showed significant trends in 42.9%, 28.6%, and 28.6% of arid regions, respectively. However, WS, Tmin, RHmin, P, and CD had significant trends in 44.4%, 22.2%, 11.1%, 11.1%, and 11.1% of semiarid regions, respectively.