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

**Mohammad Valipour 1,2, \*, Sayed M. Bateni 1 , Mohammad Ali Gholami Sefidkouhi 3 , Mahmoud Raeini-Sarjaz <sup>3</sup> and Vijay P. Singh 4**


Received: 1 September 2020; Accepted: 5 October 2020; Published: 10 October 2020

**Abstract:** Understanding the trends of reference evapotranspiration (*ETo*) and its influential 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 suffer 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 *ETo*, and thus could not represent the interactive effects of all forces driving trends of *ETo*; (3) they mostly studied trends of *ET<sup>o</sup>* and meteorological variables in similar climate regions; (4) they often did not eliminate the impact of serial correlations on the trends of *ETo* and meteorological variables; and finally (5) they did not study the extremum values of meteorological variables and *ETo*. This study overcame the abovementioned shortcomings by (1) analyzing the 50-year (1961–2010) annual trends of *ET<sup>o</sup>* and 12 meteorological variables from 18 study sites in contrasting climate types in Iran, (2) removing the effect 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 *ETo*, and (4) evaluating the coincidence of annual extremum values of meteorological variables and *ETo*. The results showed that *ET<sup>o</sup>* and several meteorological variables (namely wind speed, vapor pressure deficit, cloudy days, minimum relative humidity, and mean, maximum and minimum air temperature) had significant trends at the confidence level of 95% in more than 50% of the study sites. These significant trends were indicative of climate change in many regions of Iran. It was also found that the wind speed (*WS*) had the most significant influence on the trend of *ET<sup>o</sup>* in most of the study sites, especially in the years with extremum values of *ETo*. In 83.3% of the study sites (i.e., all arid, Mediterranean and humid regions and 66.7% of semiarid regions), both *ET<sup>o</sup>* and *WS* reached their extremum values in the same year. The significant changes in *ET<sup>o</sup>* due to *WS* and other meteorological variables have made it necessary to optimize cropping patterns in Iran.

**Keywords:** reference evapotranspiration; climate change; drought; meteorological extremes; climatic variables; wind speed

#### **1. Introduction**

Assessment of changes in reference evapotranspiration (*ETo*) 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–7]. Evaluating the trends of meteorological variables may help determine the effects of major factors on *ETo* and climate change [8,9].

Dadaser-Celik et al. [10] evaluated the 32-year trend of *ET<sup>o</sup>* 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 *ETo*. Song et al. [11] assessed the 46-year trend of *ETo* in the North China Plain. Their results indicated that the downward trends of net radiation and wind speed had a bigger impact on *ET<sup>o</sup>* compared to the upward trends of maximum and minimum air temperature. Wang et al. [12] characterized the 31-year trend of *ET<sup>o</sup>* in the western Heihe River Basin in China. They found that wind speed and sunshine duration were the two key meteorological variables that decreased *ETo*.

Darshana et al. [13] investigated the 30-year trend of *ETo* in the Tons River Basin in central India. Their outcomes showed that maximum air temperature and net radiation had a stronger impact on *ET<sup>o</sup>* compared to minimum air temperature and relative humidity. Zongxing et al. [14] studied the 49-year trend of *ET<sup>o</sup>* in the south-west of China. The decrease in wind speed was the main driving force for the reduction of *ETo*. 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 *ETo* in the Upper Mekong River Basin. They showed that sunshine duration had a more significant (at the confidence level 95%) effect on *ET<sup>o</sup>* compared to air temperature, relative humidity, and wind speed. Gao et al. [16] evaluated the 45-year trend of *ETo* in China, and found that sunshine duration, wind speed, and relative humidity had a more important influence on *ETo* than air temperature. Zhang et al. [17] assessed the changes in *ET<sup>o</sup>* and its controlling factors in China. They indicated that maximum air temperature, relative humidity, and wind speed affected *ETo*. A number of studies (e.g., [8,9,18–20]) showed that the increasing and decreasing trends of *ETo* 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<sup>o</sup>* 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 *ETo*, and found that wind speed had the most influence on *ETo*. Unfortunately, they did not consider the impact of serial correlations (autocorrelation) on the trends of *ET<sup>o</sup>* and meteorological variables. Kousari and Ahani [23] assessed trends of *ET<sup>o</sup>* 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 *ETo* from 11 weather stations in arid regions of Iran but did not specify which meteorological variables controlled trend of *ETo*.

The existing studies (1) evaluated the impact of various climatic variables on *ETo* and (2) showed the signals of climate change in Iran based on upward/downward trends of *ET<sup>o</sup>* and meteorological variables [8,9,20–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–32]. Borman [26] used a 50-year period to analyze the sensitivity of *ETo* to climate change in Germany. Kingston et al. [27] found a considerable uncertainty in evaluating the changes of *ETo* due to climate change using 30-year data. Several studies assessed the response of *ET<sup>o</sup>* to climate change in China using 50-year data [28–31]. Hence, there is a consensus in the literature to use at least a 50-year period to evaluate the trend of *ET<sup>o</sup>* due to climate change. (2) They evaluated the individual impact of a few climatic variables on *ETo*, and thus could not represent the interactive effects of all forces driving trends of *ETo*. (3) They mostly studied trends of *ETo* and meteorological variables in similar climate regions. (4) They often did not eliminate the impact of serial correlations (autocorrelation) on the trends of *ET<sup>o</sup>* and meteorological variables. Finally (5) they did not assess if the extremum values of meteorological variables and *ETo* occur simultaneously. This study overcame the abovementioned drawbacks by analyzing the 50-year (1961–2010) trends of *ETo* 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 *ETo.* 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<sup>o</sup>* and meteorological variables was studied to identify the main meteorological variables that affect variations of *ETo*.

#### **2. 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.


**Table 1.** Location, altitude, and climate of the 18 study sites.

ICAO: International Civil Aviation Organization; masl: meter above sea level.

**Figure 1.** Location of the 18 study sites in Iran.

#### **3. Models and Methods**

#### *3.1. FPM Equation*

– – *ET<sup>o</sup>* 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–37].

The FPM equation is given by [33]:

() ()

=

$$ET\_o = \frac{0.408(R\_\text{\textmu} - G) + \gamma \frac{900}{T + 273} \mu (es - ea)}{\Delta + \gamma (1 + 0.34 \mu)} \tag{1}$$

*γ* ℃ Δ – ℃ ℃ where *ET<sup>o</sup>* 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<sup>2</sup> /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]:

$$\text{es} = \frac{e\_{(T\text{max})} + e\_{(T\text{min})}}{2} \tag{2}$$

100

2

100

$$ea = \frac{e\_{(Tmin)}\frac{RHmax}{100} + e\_{(Tmax)}\frac{RHmin}{100}}{2} \tag{3}$$

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<sup>n</sup>* is the net radiation (MJ/m<sup>2</sup> /day), and is estimated by [21,33]:

$$R\_{\rm nl} = \left(a\_{\rm s} + b\_{\rm s} \frac{n}{n\_{\rm max}}\right) \mathbf{R}\_{\rm a} \tag{4}$$

where *n* is the number of sunshine hours per day (hr/day), *nmax* is the maximum possible duration of sunshine per day (hr/day), and *R<sup>a</sup>* is the extraterrestrial radiation (MJ/m<sup>2</sup> /day). *R<sup>a</sup>* and *nmax* depend on the latitude and julian day [33]. *a<sup>s</sup>* and *b<sup>s</sup>* are empirical coefficients, which are obtained from [21] for each study site.

In this study, the annual averages of daily FPM *ETo* estimates and meteorological variables were used to analyze the annual trends.

#### *3.2. 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–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–40]:

$$S = \sum\_{i=1}^{N-1} \sum\_{j=i+1}^{N} \text{sign}(\mathbf{x}\_{j} - \mathbf{x}\_{i}) \tag{5}$$

$$sign(\mathbf{x}\_{j} - \mathbf{x}\_{i}) = \begin{cases} 1 \text{ if } \left(\mathbf{x}\_{j} - \mathbf{x}\_{i}\right) > 0 \\ 0 \text{ if } \left(\mathbf{x}\_{j} - \mathbf{x}\_{i}\right) = 0 \\ -1 \text{ if } \left(\mathbf{x}\_{j} - \mathbf{x}\_{i}\right) < 0 \end{cases} \tag{6}$$

$$VAR(S) = \frac{1}{18} \left[ N(N-1)(2N+5) - \sum\_{p=1}^{g} t\_p \binom{t\_p - 1}{p} (2t\_p + 5) \right] \tag{7}$$

$$\text{If } \mathbf{S} > 0 \text{ then } Z = (\mathbf{S} - 1) \text{\textquotedbl{}VAR (S)\textquotedbl{}}{}^{0.5} \text{\textquotedbl{}}\tag{8a}$$

$$\text{If } S = 0 \text{ then } Z = 0 \tag{8b}$$

$$\text{If } S < 0 \text{ then } Z = (S+1) \text{\%} \\ \text{VAR}(S) \stackrel{0.5}{\longrightarrow} \tag{8c}$$

where *N* is the number of data, *S* is the summation of signs, *VAR* is the variance, and *x<sup>j</sup>* and *x<sup>i</sup>* are the data values in years *j* and *i*, respectively (with *j* > *i*). *t<sup>p</sup>* 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–40].

Local (at-site) significance levels for each trend test can be obtained from

$$p = 2[1 - \Phi(|Z|)]\tag{9}$$

where | | denotes the absolute value, and Φ ( ) is defined as,

$$\Phi(|Z|) = \frac{1}{\sqrt{2\pi}} \int\_0^{|Z|} \exp\left(\frac{t^2 \theta}{2}\right) dt \tag{10}$$

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].

In this study, trend-free pre-whitening (TFPW) was used to effectively eliminate the impact of serial correlations on the MK test [41,43–45].

This method is defined by the following steps:

1. Calculate the slope of trend using the Sen's slope estimator [43,44]:

$$Q\_i = \frac{\mathbf{x}\_j - \mathbf{x}\_k}{j - k} \tag{11}$$

$$Q = \begin{cases} \begin{array}{c} Q\_{\frac{N+1}{2}} \text{ N is odd} \\ \frac{1}{2} \left( Q\_{\frac{N}{2}} + Q\_{\frac{N+2}{2}} \right) N \text{ is even} \end{array} \tag{12}$$

where *Q* stands for the slope of trend. *Q<sup>i</sup>* is the Sen's slope estimator for each value of *i*, *x<sup>j</sup>* and *x<sup>k</sup>* 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:

$$\mathbf{X}\_t^{'} = \mathbf{X}\_t - \mathbf{Q}t \tag{13}$$

where *X<sup>t</sup>* ′ is the lag-1 autocorrelation of the de-trended time series and *X<sup>t</sup>* 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]:

$$r\_j = \frac{\frac{1}{n-j} \sum\_{i=1}^{N} \left( \mathbf{X}\_i - \overline{\mathbf{X}} \right) \left( \mathbf{X}\_{i+j} - \overline{\mathbf{X}} \right)}{\frac{1}{n} \sum\_{i=1}^{N} \left( \mathbf{X}\_{i+j} - \overline{\mathbf{X}} \right)^2} \tag{14}$$

where *r<sup>j</sup>* 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:

$$X\_{t}^{\prime\prime} = X\_{t}^{\prime} - r\_{j} X\_{t-1}^{\prime} \tag{15}$$

4. Adding the removed trend in step 2:

$$X\_{t}^{\prime\prime\prime} = X\_{t}^{\prime\prime} + \mathcal{Q}t \tag{16}$$

#### *3.3. Spearman's Rho Test*

The Spearman's Rho test was applied to investigate correlation among all the meteorological variables and *ETo* estimates from the FPM equation.

The Spearman's Rho test is defined as [18]:

$$\rho = 1 - \frac{6\sum\_{i=1}^{N} (\mathbf{x}\_i - \mathbf{i})^2}{N(N^2 - 1)}\tag{17}$$

$$Z = \rho \sqrt{\mathbf{N} - 1} \tag{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*<sup>α</sup> at a significance level of α, then the null hypothesis of no trend is rejected [18].

Following the literature [1–24], the confidence levels of 90%, 95%, and 99% were used in this study to evaluate the trends of *ETo* 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.

#### **4. Results and Discussions**

#### *4.1. Comparison of p-Values from the MK and TFPW-MK Tests*

Table 2 compares *p*-values from the MK and TFPW-MK tests for *ET<sup>o</sup>* and the 12 meteorological variables in the 18 study sites.

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–49].


**Table 2.**Estimated*p*-values from the Mann–Kendall (MK) and trend-free pre-whitening Mann–Kendall (TFPW-MK) tests in the 18 study sites.

#### *4.2. Trends of ETo 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 *ETo*. However, herein, the trends are presented only in two sites, Kerman and Qazvin (Figures 2 and 3). –

– **Figure 2.** The 50-year (1960–2010) trends of annually averaged reference evapotranspiration and meteorological variables in Kerman. Circles indicate an insignificant trend. One, two, and three triangles indicate a significant trend at a confidence level of 90%, 95%, and 99%, respectively. The value on the left of triangular symbols shows the increasing/decreasing rate over the 50-year period.

– –

–

−

– **Figure 3.** The 50-year (1960–2010) trends of annually averaged reference evapotranspiration and meteorological variables in Qazvin. Circles indicate an insignificant trend. One, two, and three triangles indicate a significant trend at a confidence level of 90%, 95%, and 99%, respectively. The value on the left of triangular symbols shows the increasing/decreasing rate over the 50-year period.

#### 4.2.1. Kerman

Figure 2 shows the 50-year trends in the annual averages of daily FPM *ET<sup>o</sup>* 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.

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 *ETo* estimates.

Although the upward trends in *Tmax, Tmin, Tmean*, and *n* could potentially increase *ETo*, a downward trend in *WS* and (*es-ea*) led to a significant decreasing trend in the FPM *ET<sup>o</sup>* values. It is worth mentioning that the minimum of both *WS* and *ET<sup>o</sup>* 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 *ETo* in Kerman.

Although *CD* showed the strongest downward trend among all the meteorological variables in Kerman, it could not capture the variations of *ETo* as good as *WS*. Both *WS* and *ETo* 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 *ETo*.

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–57]. These results show the complexity of forces driving the trend of *ET<sup>o</sup>* and may highlight importance of considering different *ET<sup>o</sup>* models and variables in various climates for local and global studies.

#### 4.2.2. Qazvin

Figure 3 shows the 50-year trends of annually averaged *ET<sup>o</sup>* 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 *ETo* estimates. It is worth mentioning that these significant trends reduced *ET<sup>o</sup>* and may cause alarming climate change in Qazvin [57–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<sup>o</sup>* 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].

#### *4.3. Trends of ETo and Meteorological Variables in Iran*

Table 3 and Figure 4 summarize trends of FPM *ETo* estimates and meteorological variables in the 18 study sites. The FPM *ETo* 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<sup>o</sup>* 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.


**Table 3.** Trends of Agricultural Organization of the United Nation (FAO)-Penman–Monteith (FPM) evapotranspiration (*ETo*) estimates and 12 meteorological variables in the 18 study sites. U and D represent the significant upward and downward trends, respectively. NS denotes no significant trend.

\*, \*\*, \*\*\* are significant trends at a confidence level of 90%, 95%, and 99%, respectively.

**Figure 4.** The average of variations in reference evapotranspiration and meteorological variables in Iran from 1961 to 2010. Big, medium, and small triangles imply significant values at a confidence level of 99%, 95%, and 90%, respectively. Blue, red, yellow, and green colors indicate Mediterranean, arid, semiarid and humid climates, respectively.

Results showed that *ETo* 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.

As shown in Table 3, in eight of the study sites (i.e., 44.4% of sites), there were significant trends in both *ET<sup>o</sup>* and *Tmean*/*RHmin*/*WS*. Moreover, both *ET<sup>o</sup>* 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 *ETo* in Iran.

*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<sup>o</sup>* in 55.5% of the study sites. Thus, *WS* may be the most important variable that controls the variations of *ETo*, 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 *ETo*. As shown, *WS* has the highest significant correlation with *ET<sup>o</sup>* in all the study sites expect Rasht (RA). Studying the trends of *ET<sup>o</sup>* 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<sup>o</sup>* 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.

According to the results, *P*, *WS*, and *CD* showed significant trends in 42.9%, 28.6%, and 28.6% of the arid sites, respectively. *WS, Tmin, RHmin, P*, and *CD* had substantial trends in 44.4%, 22.2%, 11.1%, 11.1%, and 11.1% of the semiarid sites, respectively. Furthermore, *WS, Tmin,* and *CD* had a rate of more than ±2%/decade in 66.7% (85.7% of arid regions vs. 66.7% of semiarid regions), 66.7% (42.9% of arid regions vs. 77.8% of semiarid regions), and 55.6% (57.1% of arid regions vs. 55.6% of semiarid regions) of the study sites, respectively.

As mentioned above, variations of *ETo* 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<sup>o</sup>* 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<sup>o</sup>* 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 *ETo* and other climatic variables over long periods.

#### *4.4. Range of Variations of ETo and Meteorological Variables*

Figure 5 shows the box plots for the 50-year variations of *ET<sup>o</sup>* and meteorological variables in each study site.

**Figure 5.** Statistical boxplots of meteorological variables and reference evapotranspiration in each study site. For each distribution, the horizontal line within the box indicates the median (50% percentile). The upper and lower edges of the box represent the 75% and 25% percentile, respectively. The upper and lower ends of the whiskers represent the maximum and minimum values. Outliers are observations beyond the end of whiskers.

As can be seen, the largest 50-year variations in both *ET<sup>o</sup>* and *WS* happened in Zabol (ZA), which highlighted the significant influence of *WS* on *ETo*.

The second highest variations of *ETo* was seen in Ahvaz (AH). AH had also the highest variations of *Tmean* compared to the other study sites. Hence, the variations of *ET<sup>o</sup>* in Ahvaz can be related to that of *Tmean*.

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<sup>o</sup>* [9,10]. Most of the outliers occurred in 2010, which signaled a drought condition in different regions of Iran. The FPM *ET<sup>o</sup>* values ranged from 2 to 11 mm day−<sup>1</sup> in various climates of Iran.

As shown in Equation (1), *ET<sup>o</sup>* is affected by the atmospheric vapor pressure deficit (*es-ea*). As expected, a high positive correlation exists between *ET<sup>o</sup>* and *es-ea* in Figure 5. For example, high values of *ETo* and *es-ea* are seen in Zabol (ZA). Moreover, the lowest *ETo* 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 *ETo* and *Tmax*, and *ETo* 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<sup>n</sup>* (Equation (4)) and *ET<sup>o</sup>* (Equation (1)). On the other hand, a negative correlation is found between *ET<sup>o</sup>* and *RH*. The highest *ET<sup>o</sup>* 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.

#### *4.5. Annual Extremum Values of ETo and Meteorological Variables during the Study Period (1961–2010)*

Table 4 shows the years of occurrence of maximum and minimum values of *ETo* and meteorological variables in each study site during the 50-year study period (i.e., 1961–2010).


**Table 4.** The years of occurrence of maximum and minimum values of *ETo* and meteorological variables in each study site during the 50-year study period. If the annual extremum values of *ETo* and a particular meteorological variable coincide, they are highlighted in green color.

If the annual extremum values of *ETo* and a particular meteorological variable coincide, they are highlighted by green color (Table 4).

As can be seen, the maximum of both *ETo* 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<sup>o</sup>* 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<sup>o</sup>* and *WS* occurred in the same year. These results imply that *WS* is the primary driver of the variations in *ETo* in Iran. In addition, the annual extremum values of *ETo* 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<sup>o</sup>* in these study sites.

Figure 6 illustrates the frequency of maximum and minimum values of *ET<sup>o</sup>* and meteorological variables (namely *Tmean, Tmax, Tmin, Tmax-Tmin, es-ea, P, RH, RHmin, CD, WS, WD,* and *n*) in all the study sites for each year. Higher frequencies are observed in 1972, 1982, and 1992, indicating a large number of study sites reached their maximum and minimum values of *ETo* and meteorological variables in these years. This happened because 1972, 1982, and 1992 are El Niño-Southern Oscillation (ENSO) years [64–66]. The highest frequency during the 50-year period is observed in 2010. This year was considered as a dry year in Iran and most regions of Iran experienced drought alarms [50,52–57].

**Figure 6.** Frequency of maximum and minimum values of *ETo* and 12 meteorological variables in 18 study sites for each year.

–

*ρ ≥*

#### *4.6. Correlation between ETo and Meteorological Variables*

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 (ρ) maps for ≤ ρ ≤1, 0.5 ≤ ρ 0.8, and ρ **Figure 7.** Spearman's Rho correlation coefficient (ρ) maps for weather stations located in arid, Mediterranean, and humid regions. 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.

Spearman's correlation coefficient (ρ) maps for ≤ ρ ≤ ≤ ρ 0.8, and ρ **Figure 8.** 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.

As indicated, *WS* was the only meteorological variable that had high correlation (ρ ≥ *0.8*) with *ET<sup>o</sup>* 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 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–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<sup>o</sup>* 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 *ETo*. It also lists the three meteorological variables that have the highest significant correlation with *ETo*. As can be seen, *ET<sup>o</sup>* had the highest significant correlation with *WS* in all the sites except Rasht. *ET<sup>o</sup>* 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 *ETo*. This highlights the complexity of the forces driving variations of *ETo*.


**Table 5.** The total number of meteorological variables in each site with significant correlations with *ETo*, and the three meteorological variables with the highest significant correlation with *ETo*.

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<sup>o</sup>* and meteorological variables across Iran. These significant trends are indicative of climate change signals [1–24]. For instance, in the Mashhad site, significant upward (downward) trends were seen for *ETo*, *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–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–23]. In the Shiraz site, significant upward (downward) trends of *ETo*, *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 *ETo*, *Tmax*, *WS*, *n*, *Tmean* and *Tmax* (*RHmin* and *CD*) signified climate change [89–92]. Finally, in the Zabol site, upward (downward) trends of *ETo*, *Tmean*, *Tmin*, and *WS*, and *Tmax* (*RHmin*) may be considered as alarms of climate change [84,93].

The FPM *ETo* 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<sup>2</sup> and the surface albedo, ultimately causing uncertainty in the *ET<sup>o</sup>* estimates. Surface albedo varies with changes in vegetation and soil moisture [94].

Elevated CO<sup>2</sup> emission influences plant physiology by diminishing stomatal conductance [95–97]. Ignoring changes in the surface albedo and vegetation response to the elevated CO<sup>2</sup> emission lead to uncertainty in the FPM *ET<sup>o</sup>* estimates and results of this study. A similar uncertainty in the FPM *ET<sup>o</sup>* estimates was reported by other studies. For example, Milly and Dunne [98] showed that ignoring the stomatal resistance response to the increased CO<sup>2</sup> emission in the FPM equation led to discrepancy between the *ET<sup>o</sup>* estimates from climate models and the FPM equation. In a similar effort, Li et al. [99]

indicated that the variations in the FPM *ET<sup>o</sup>* 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<sup>2</sup> emission. For example, Domec et al. [100] reported that vegetation responses of pine to long term elevated CO<sup>2</sup> were manifested only when soil moisture was at a high level. Future studies should be directed toward evaluating the uncertainties in the FPM *ETo* estimates with respect to climate change and drought condition [101,102].

#### **5. 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 (*ETo*) in the 18 study sites in different climate regions of Iran.

The effect of meteorological variables on *ET<sup>o</sup>* 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<sup>o</sup>* 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.

In 83.3% (all arid, Mediterranean and humid regions and 66.7% of semiarid regions) of the study sites, *ET<sup>o</sup>* and *WS* reached their maximum and/or minimum in the same year(s). Furthermore, *ET<sup>o</sup>* had the highest significant correlation with *WS* in all the sites except Rasht. This indicates that *WS* is the primary driver of the variations in *ET<sup>o</sup>* in most weather stations of Iran, especially in the years with extremum values of *ETo*. However, in a number of weather stations, the significant correlation between *ET<sup>o</sup>* and other meteorological variables such as *Tmax, Tmax-Tmin, es-ea, P, RH, RHmin,* and *n* made it difficult to characterize driving forces for the trend of *ET<sup>o</sup>* and signals of climate change. The results of this study highlighted the complexity of forces that drive variations of *ETo*.

**Author Contributions:** M.V. conceived and designed the study, performed the calculations, and prepared the original draft. S.M.B. supervised the study and revised the manuscript. M.A.G.S. and M.R.-S. supervised the study. V.P.S. reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

1. McNally, A.; Arsenault, K.; Kumar, S.; Shukla, S.; Peterson, P.; Wang, S.; Funk, C.; Peters-Lidard, C.D.; Verdin, J.P. A land data assimilation system for sub-Saharan Africa food and water security applications. *Sci. Data* **2017**, *4*, 170012. [CrossRef] [PubMed]


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
