Reference evapotranspiration (ET₀) refers to the regional evapotranspiration capacity under sufficient water supply conditions, defined as the evapotranspiration rate on the hypothetical underlying surface with fixed parameters (Allen et al., 1998). ET₀ is a key element of estimating actual evapotranspiration (Fan et al., 2016; Liu et al., 2017), hydrological modeling, and predicting the impact of changes in meteorological conditions on water supply (Liu et al., 2010). In the context of global climate change, the accurate estimation of ET₀ is indispensable for improving water management, strengthening crop efficient water use, and predicting crop agricultural water demand (Han et al., 2014; She et al., 2017a; Shiri, 2017). The revised Penman-Monteith equation recommended by the Food and Agriculture Organization (FAO) in 1998 is one of the most widely methods used to calculate ET₀ (Tabari et al., 2012; Kisi, 2016; Wang et al., 2017). From the perspectives of meteorology and hydrology, ET₀ is also an important topic in studies of the energy balance and water balance on land surface (Liu and Sun, 1999). Quantifying the effects of different meteorological variables on ET₀ facilitates the explanation of the impact of climate change on hydrological processes in terrestrial ecosystems (Qin et al., 2017).

As global warming is exacerbating recently, ET₀ has shown an increasing trend owing to its high dependence on temperature in some regions in the world, e.g., the Mississippi River Basin (Qian et al., 2007a), the Mediterranean region (Palumbo et al., 2012; Vicente-Serrano et al., 2014), the Taiwan Province of China (Yu et al., 2002), the Platte River Basin in central Nebraska, USA (Irmak et al., 2012), and Florida, USA (Abtew et al., 2011). Although air temperature is increasing continuously, evapotranspiration is observed to a decrease trend but not an expected increase trend in many regions (Hobbins et al., 2004; Roderick and Farquhar, 2005; Burn and Hesch, 2007; Bandyopadhyay et al., 2009; Jhajharia et al., 2012; Hosseinzadeh Talaee et al., 2014; Zheng and Wang, 2015). Such a contradiction between the expectations and the observations is called the ‘evaporation paradox’ (Roderick and Farquhar, 2002).

Furthermore, spatial-temporal changes of ET₀ are complex in China (Xing et al., 2016). Based on wide reviews of the literatures, it is suggested that: 1) the changing trends of ET₀ are different in China and other regions of the world; 2) the meteorological factors that affect ET₀ include multiple aspects, such as wind speed, precipitation, and solar radiation in addition to air temperature, and affect ET₀ to varying degrees; 3) the major meteorological factors that cause ET₀ changes in different regions of China are different (Table 1); 4) few studies on the Loess Plateau (LP) have considered the effect of vapor pressure deficit on ET₀; and 5) detecting meaningful change points is very important for the analysis of long time series in the field of hydrometeorology. Therefore, it is inevitable to analyze the trends of the ET₀ changes in a long-term period and different periods and the comprehensive effects of various meteorological variables on ET₀ on the LP.

The LP is characterized the serious soil erosion and the fragile ecological environment. To solve the deteriorating ecological problems on the LP, revegetation was carried out several decades ago, which led to a significant increase in the average vegetation coverage and a significant change in land use patterns (Lü et al., 2012). Revegetation also resulted in an ecological problem associated with water resource depletion (Feng et al., 2016). Meanwhile, vegetation greening also further affect the changes of meteorological variables and ET₀ (Zhang et al., 2011). Atmospheric demand becomes the main influencing factor in the process of the ET₀ changes instead of the air temperature alone (Qin et al., 2017). Moreover, due to the decadal and interannual changes in the climate variables and human activities, the ET₀ series may have multiple change points and different changing trends in different periods instead of following a simple monotonic linear trend, which may also lead to differences in the contributing factors of ET₀ in each period (She et al., 2017a). Previous studies on the LP focused on detecting the long-term trend of ET₀ and attributing it to independent meteorological factors, which may affect the understanding of the ET₀ changes and their contributing factors. Therefore, it is concluded in this paper that under the regional background of vegetation greening and regarding vapor pressure deficit as a key contributing factor, it is of great significance for regional water resource planning as well as ecosystem restoration to identify the changing trend of ET₀ in different periods and its major contributing factors in the process of studying the trend of ET₀ changes over a long period.

The main objectives of this study are to: 1) detect the change points of the annual ET₀ series from 1960 to 2017 by utilizing the segmented regression model; 2) analyze the spatiotemporal distribution of ET₀ on annual and seasonal scales on the LP by using the linear regression method and the Mann-Kendall test and further explain the trend of the ET₀ changes in different periods; 3) analyze the periodicity of ET₀ under climate change conditions through wavelet analysis; 4) assess the relative contribution of the various meteorological factors to the ET₀ changes and the temporal evolution of the relative contribution utilizing multiple stepwise regressive analysis and the generalized linear model; and 5) quantitatively identify the dominant meteorological variables of the ET₀ changes in different periods using the detrending method.

The Loess Plateau (LP) (33°43′N–41°16′N; 100°54′E–114°33′E) is a typical region with water and soil losses and ecological fragility in the northern China (Fig. 1). The LP has a total area of 6.4 × 10⁵ km² (Liu et al., 2016). It has a continental monsoon climate characterized by hot summers and autumns with rainfall dominated by rainstorms, and cold and dry winters and springs with winds and sandstorms. From northwest to southeast, the annual mean temperature ranges from 3.6°C to 14.3°C, and the multiyear average annual precipitation ranges from 150 mm to 800 mm and mostly occurs in summer (Wang et al., 2019). Considering the differences in ET₀ contributing factors in different climatic regions, this study divided the Loess Plateau region into three regions according to the climatic conditions: (Ⅰ) the semihumid region, (Ⅱ) the semiarid region, and (Ⅲ) the arid region (Fig. 1).

Figure 1.  Topographic map of the study areas and meteorological stations: (I) the semihumid region, (II) the semiarid region, and (III) the arid region

The daily meteorological data of 69 meteorological stations on the LP collected by the China Meteorological Administration (CMA) (http://data.cma.cn/). They were applied to calculate of ET₀ from 1960 to 2017 by the FAO56-PM model. The meteorological variables include maximum air temperature, minimum air temperature, average air temperature, wind speed, mean relative humidity, and sunshine duration (Allen et al., 1998).

Because of the difficulties in obtaining accurate field measurements, ET₀ is generally calculated from the weather data. In this study, ET₀ was calculated using the FAO Penman-Monteith model, which was widely used all over the world (Vicente-Serrano et al., 2014; Kisi, 2016; Maček et al., 2018; Zhao et al., 2020). The FAO recommended it as the standard method for the estimation of ET₀ (Allen et al., 1998). The Penman-Monteith equation can be expressed as follows:

where ET₀ is the reference evapotranspiration (mm/d); G is the soil heat flux density (MJ/(m²·d)); R_n is the net radiation at the crop surface (MJ/(m²·d)); T is the average air temperature at 2 m height (°C); U₂ is the wind speed at 2 m height (m/s); e_s is the saturation vapor pressure (kPa); e_a is the actual vapor pressure (kPa); e_s− e_a is the saturation vapor pressure deficit (kPa); Δ is the slope of the vapor pressure curve (kPa/°C); and γ is the psychrometric constant (kPa/°C) (Allen et al., 1998).

The wind speed at 2 m height can be calculated from the wind speed at 10 m by using the equation based on the wind profile relationship as follows (Allen et al., 1998):

where z is the height of the measurement above the ground surface (m) and u_z is the measured wind speed at z m above the ground surface (m/s). Here, z = 10 m.

The global solar radiation (R_s) can be estimated by the following equation:

where R_a is the extraterrestrial radiation (MJ/(m²·d)), n is the actual sunshine duration, N is the potential sunshine duration, and a and b are empirical coefficients of 0.25 and 0.5, respectively (Allen et al., 1998).

The vapor pressure deficit (VPD) can be calculated by the following equation:

In the above equation, e_s (kPa) and e_a (kPa) can be calculated as follows:

where T_max (°C) and T_min (°C) are the maximum and minimum air temperatures, respectively, and RH is the mean relative humidity (%).

The linear regression method was used to detect and analyze temporal trends (Feng et al., 2017). Generally, the least squares method was applied to estimate the linear trend as follows:

where $ {\widehat{y}}_{n} $ is the response variable; t is the year in the series; a₀ is the intercept; and a₁ is the slope of the estimated trend. If a₁ < 0, then the trend is negative; and if a₁ > 0, then the trend is positive. a₁ × 10 is called the climate tendency rate, which is used to calculate the climate change rate per decade.

The segmented regression mode, which can detect step changes and shift trends simultaneously even when the specific pattern of change points is lacking (Shao and Campbell, 2002; Shao et al., 2010; She et al., 2017b), was employed to identify the multiple potential change points of the annual ET₀ between 1960–2017 in this study. The least squares method was used to fit the model, and the modified Akaike information criterion (AICc) (Hurvich and Tsai, 1989) was utilized to find the suitable break points (Shao and Campbell, 2002; Shao et al., 2010). When the AICc has the minimum value, the model is optimal. The linear regression method was applied to identify the trends of the annual ET₀ series in different segmented periods after detecting the break points (She et al., 2017a).

To analyze change trends in a series of meteorological variables and ET₀, the Mann-Kendall test (Mann, 1945; Kendall, 1975), a typical nonparametric statistical analysis, was used here. Details of the Mann-Kendall test are provided as follows.

Given a data series composed of x₁, x₂, …, x_n that is independent and identically distributed, the statistic S of Kendall’s tau is defined as follows:

where n is the length of the dataset, x_i and x_j are the sequential data values in time series i and j, and

Mann (1945) and Kendall (1975) documented that when n ≥ 8, the statistic S is approximately normally distributed, and the mean and variance are as follows:

where t_m is the number of ties of extent m. The standardized test statistic Z_c can be computed according to the following equation:

where Z_c is the test statistic. When |Z_c| > Z_{1 − α/2}, where Z_{1 − α/2} is the standard normal deviation and α is the significance level for the test, the data series has a significant trend at the significance level of α. In this study, a significance level of 0.05 was used. The magnitude of the trend of the data can be obtained using Theil-Sen’s estimator as follows:

where β is the estimated slope of the data series (Sen, 1968). If β < 0, then the data series has a decreasing trend, and if β > 0, then the data series has an increasing trend.

Using a cluster of wavelet functions to represent or approximate a certain signal or function is the basic idea of wavelet analysis. Wavelet analysis is a tool for studying the long-term changes in hydrometeorological elements, and it can not only clearly reveal the various time-varying periodic oscillations hidden in the time series but also reflect their changing trends (Wang et al., 2014; Feng et al., 2017). The wavelet transform of time series f(t) is defined as follows (Feng et al., 2017):

where W_f(a, b) is the transformation coefficient; a is a scale factor representing the cycle length; b is a time factor representing the time location; t is the time variable; ψ(t) is the mother wavelet; and $ \overline{\mathrm{\Psi }} $ is the conjugate function of ψ(t).

Since the observed time series are mostly discrete, the discrete wavelet transform can be defined as follows:

where f (kΔt) is the time series; k = 1, 2, …, N; and Δt is the time interval.

Wavelet variance is defined as follows:

The wavelet transform reflects the characteristics of the time domain and frequency domain of f(t) simultaneously (Cazelles et al., 2008). If a is decreased, then the frequency resolution declines, but the time resolution improves; if a is increased, then the frequency resolution becomes better and shifts toward low frequencies, but the time resolution becomes poorer.

Some spatial interpolation methods, such as spline, ordinary kriging (OK) and inverse distance weighting (IDW), have been used to interpolate hydro-meteorological variables, including evapotranspiration (Ashraf et al., 1997; Dalezios et al., 2002; Zhang et al., 2010). The spline method is characterized by fitting a smooth and continuous surface with the observation points and does not need a preliminary estimate for the structure of temporal variance and statistical hypothesis (Zhu et al., 2012). The interpolated surface passes exactly through the data points and has a minimum curvature (Zuo et al., 2012). Previous investigators have demonstrated that spline is a better interpolation method in climatic and meteorological research (Hutchinson and Gessler, 1994; Price et al., 2000; Zhu et al., 2012; Shan et al., 2015). Hutchinson and Gessler (1994) found that the spline method was approximately as accurate for interpolation as kriging, but avoided initial estimation of the covariance structure. There were no singularities at the data points. Price et al. (2000) compared the inverse distance weighting and spline method by calculating the relative values of the difference and showed that the latter produced better results. Zhang et al. (2010) used the spline method based on the lower square root of cross-validation errors it produced compared with some other methods and settings based on previous work. Shan et al. (2015) concluded that the results from spline were more favorable around topographic features than other interpolation methods by summarizing previous studies and therefore chose this method. Li et al. (2016) firstly interpolated the values to map the spatial distribution of the evapotranspiration based on the three methods; comparisons were made between the accuracy of the results calculated using three spatial interpolation techniques using cross-validation. The results indicated that the accuracy of the different methods was satisfactory on the Loess Plateau. Besides, the spline method fits a smooth and continuous surface and has a minimum curvature. Therefore, the spline method was applied in this study.

To evaluate the relative contribution of meteorological variables to the ET₀ trends, a stepwise linear regression analysis and a generalized linear model were utilized in this study, with ET₀ as the dependent variable, and sunshine duration (n), mean air temperature (T), wind speed (U₂), mean relative humidity (RH), precipitation (Pr), maximum air temperature (T_max), minimum air temperature (T_min), vapor pressure deficit (VPD), solar radiation (R_s), and net radiation (R_n) as the independent variables. Stepwise linear regression is a multiple linear regression analysis method. It was applied to remove the variables that cause multicollinearity and select effective predictors (Cristea et al., 2013; Shan et al., 2015), and the selected predictors were considered the main meteorological variables that affected the changes and variations in ET₀. The generalized linear model, which is a direct extension of the general linear model, is an effective statistical model for processing discrete observational data as well as an effective method of separating the relative contributions of individual independent variables from that of the dependent variable (Wang et al., 2019). Prior to applying these methods, the original data and the meteorological variables were normalized to eliminate the effects of the dimension and magnitude of the data.

To quantify the contributions of changes in meteorological variables to ET₀, previous studies (Liu et al., 2010; Qin et al., 2017) provided an effective method of detrending based on the following three steps: 1) removing the change trends in the meteorological variables to make them stationary data series; 2) recalculating ET₀ using one detrended data series of meteorological variables and the other original meteorological variables; and 3) comparing the recalculated ET₀ with the original ET_0, and the difference is considered the influence on the trend by that meteorological variable, which could be quantified by an evaluating indicator R:

where n is the length of the data series and ET₀^O and ET₀^R are the original and recalculated ET_0, respectively. When R > 0, the change in that meteorological variable has a positive effect on the ET₀ changes; when R < 0, the change in that meteorological variable has a negative effect on the ET₀ changes. The larger the value of |R| is, the greater the effect on ET₀ changes will be.

To remove the change trend in the meteorological variables, we adopted a detrend method based on the following steps: 1) regarding the first year 1960 as the baseline; 2) calculating the linear yearly trend in the meteorological data series; and 3) subtracting an arithmetic progression with the linear slope as the tolerance from the original data sequence to obtain the stationary time series.

The specific process of the present research and the data used are shown in the technology roadmap (Fig. 2).

Figure 2.  Technology roadmap. T is average air temperature; T_max is maximum air temperatures; T_min is minimum air temperatures; n is actual sunshine duration; U₂ is wind speed at 2 m height; RH is mean relative humidity; VPD is vapor pressure deficit; R_s is global solar radiation; and R_n is the net radiation

The temporal changes in annual ET₀ showed a fluctuating downtrend on the LP from 1960–2017 (Fig. 3a). The annual ET₀ was 984 mm on average and ranged from 898–1065 mm (Table 2). The dominant monotonic trend in the annual ET₀ showed that the ET₀ declined on the LP and in the semihumid region with tendency rates of –1.80 mm/10 yr and –15.68 mm/10 yr, respectively, while it increased in the semiarid region and the arid region with tendency rates of 1.79 mm/10 yr and 3.57 mm/10 yr, respectively (Fig. 3). The results show that the change trends of the annual ET₀ were inconsistent in the three subregions from 1960–2017 and that the annual ET₀ in the semihumid region exhibited a significant downward trend, which dominated the ET₀ changes on the LP.

Figure 3.  Change points and segmented trends of the annual ET₀ from 1960–2017 in (a) the Loess Plateau, (b) the semi-humid region, (c) the semi-arid region, and (d) the arid region. The dots denote the observed values; the red line denotes the trend by the segmented regression model; the blue dashed line denotes the trend by the simple linear regression model; the gray line denotes the five-year moving average; and the vertical lines indicate the change points

To thoroughly analyze the change patterns of ET₀, a segmented regression method was performed to accurately detect change points for the annual ET₀ series from 1960–2017, and a linear regression was performed to identify the change trends in different periods. Finally, four change points were found in 1972, 1990, 1999, and 2010, and an increasing-decreasing-increasing-decreasing-increasing change pattern of the annual ET₀ was found in the three subregions and on the LP (Table 2 and Fig. 3). It indicated that the annual ET₀ showed different change trends in each subperiod on the LP, instead of a monotonic trend over a long period. The results agreed with the previous research conclusions demonstrating that ET₀ or pan evaporation had abrupt changes in similar years of 1972, 1990, and 1999 in China (Liu et al., 2011; Li et al., 2016; Xing et al., 2016; She et al., 2017a). A new change point was detected in 2010 because the research time series was longer than that in previous studies. In the five subperiods, the highest value of ET₀ was found in the arid region, while the lowest was found in the semiarid region (Table 2). The average annual ET₀ from 1972 to 1990 was very close to the average value from 1990 to 1999, and the average annual ET₀ between 1999 and 2010 was very close to the average value between 2010 and 2017 (Table 2).

Fig. 4 shows the spatial distribution of the annual ET₀ trends from 1960–2017. A negative trend was found in the southeast, west and north of the LP. Only 24 stations (approximately 34.8% of the total number of stations) had significant downward trends at the 95% confidence level, and they were mainly distributed in the north and west of the semihumid and semiarid regions. Baotou in Inner Mongolia showed the largest decreasing trend (–38.2 mm/10 yr, P < 0.01). In contrast, the western part of the arid region and most of the semiarid region showed an increasing trend. The stations showing positive trends were located in semiarid and arid regions. However, only 17 stations (approximately 24.6% of the total number of stations) achieved a significant level (P < 0.05), and the Wutaishan station in Shanxi had the largest increasing trend (43.0 mm/10 yr, P < 0.01).

Figure 4.  Spatial distribution of the ET₀ trends on the Loess Plateau from 1960–2017. The bar chart denotes different significance station numbers

The annual ET₀ trends on the LP indicated that there was a certain periodicity of interannual variability in ET₀. This situation seemed to be due to changes in relevant meteorological factors. From 1960–2017, T, T_max, T_min, and VPD showed a significantly increasing trend, while n, RH, U₂, Pr, R_s, and R_n of most stations exhibited a significantly decreasing trend (Table S1). In general, higher air temperature and vapor pressure deficit and lower relative humidity augment ET₀, while lower wind speed and solar radiation lessen ET₀ (Allen et al., 1998; Li et al., 2012). Finally, the counteracting effects of meteorological variables resulted in a zigzag change pattern of ET₀ on the LP.

The four seasons were distinct on the LP, and the ET₀ change trends varied considerably in different seasons (Fig. 5). The spring and winter ET₀ had increasing trends at rates of 2.34 mm/10 yr and 0.85 mm/10 yr, respectively, but the summer and autumn ET₀ showed decreasing trends at rates of –4.18 mm/10 yr and –0.54 mm/10 yr, respectively (Fig. 5). The summer ET₀ had the largest values (405 mm on average) among the four seasons, and it varied from 371 to 445 mm and had the largest decreasing rate, while the downward trend of ET₀ was gentler in autumn (Table 3 and Fig. 5b). Table 3 shows that the change trends of ET₀ in autumn and winter were diverse in different climate regions. The winter ET₀ in the semihumid region showed an insignificant downward trend, and the autumn ET₀ in the semiarid and arid regions increased insignificantly, which is different from the seasonal variation trends on the entire LP.

Figure 5.  Time evolution of ET₀ on the Loess Plateau from 1960–2017 in (a) spring, (b) summer, (c) autumn, and (d) winter. The dots denote the observed values; the blue dashed line denotes the trend by the simple linear regression model; and the gray line denotes the five-year moving average

To explore the change cycle of the annual ET₀ on the LP, a Morlet wavelet analysis was performed to analyze the scale and oscillation characteristics of the annual ET₀ series from 1960–2017. The annual ET₀ had periods of 6 yr and 11–13 yr on the LP (Fig. 6). The largest peak corresponded to the time scale of 6 yr, indicating that the periodic oscillation at approximately 6 yr was the strongest; therefore, 6 yr was the first major cycle of the annual ET₀. The second peak corresponded to the time scale of 11–13 yr, which was the second major cycle. The positive and negative phases of the wavelet coefficients alternately changed at time scales of 6 yr and 11–13 yr, respectively, and the catastrophic points of the wavelet coefficients occurred approximately every 2–3 yr and 5–6 yr, which means that the change trend of the annual ET₀ changed once. The Morlet wavelet analysis shows that the annual ET₀ exhibited a short cycle of 6 yr and a medium cycle of 11–13 yr, but no long cycle was discovered on the LP.

Figure 6.  Morlet wavelet results of the annual ET₀ from 1960–2017 on the Loess Plateau: (a) variance of the wavelet coefficient, (b) wavelet coefficients, and (c) real part of the wavelet coefficients. The solid lines in Fig. 6c are values of wavelet coefficients and represents flow periods

The spatial distribution of the long-term average annual ET₀ on the LP from 1960–2012 is shown in Fig. 7. The average annual total ET₀ varied between 693 mm and 1185 mm. The spatial distribution of the annual ET₀, which was similar to that of VPD (Fig. S1h), had strong gradients with higher values in the southeast and northwest and lower values corresponding to the middle, southwest and northeast of the LP (Fig. 7a).

Figure 7.  Spatial variability of the (a) annual, (b) spring, (c) summer, (d) autumn, and (e) winter ET₀ from 1960–2017 on the Loess Plateau. Ⅰ, Ⅱ and Ⅲ are the same as Table 2

Figure S1.  Spatial distribution of (a) sunshine duration, (b) temperature, (c) wind speed at 2m height, (d) relative humidity, (e) precipitation, (f) maximum air temperature, (g) minimum air temperature, (h) vapor pressure deficit, (i) solar radiation, and (j) net radiation. Ⅰ is the semihumid region; Ⅱ is the semiarid region; Ⅲ is the arid region

From spring to autumn, the spatial distribution of ET₀ was close to that of annual ET₀ and initially decreased and then increased from northwest to southeast (Fig. 7b, c, and d). In spring, the maximum ET₀ (372.7 mm) appeared in the arid region, and the low value shifted to the semiarid region (Fig. 7b). The spatial gradient of ET₀ was strongest in summer. Almost all ET₀ values in the arid region were above 450 mm, while the lowest value (only 203.7 mm) occurred in the semiarid region (Fig. 7c). The highest autumn ET₀ appeared in the semihumid region (Fig. 7d). Compared with the spatial distribution of other seasons, the ET₀ value of the whole area in winter was lower (Fig. 7e). The mean seasonal ET₀ values from high to low were 41.2%, 31.9%, 117.9%, and 9.0% of the annual ET₀ in summer, spring, autumn, and winter, respectively.

ET₀ is mainly affected by changes in meteorological variables. Table 4 illustrates the regression coefficients of the meteorological variables on the LP and in the three subregions. The mean R² value of the regression equations was 0.946, which indicated a good fit for the model. Among the ten meteorological variables, the regression coefficients of T, U₂, VPD, R_s, and R_n were greater than 0, which indicated that the changes in ET₀ were positively correlated with these variables but negatively correlated with n; while RH, Pr, T_max, and T_min presented variations in different climatic regions. The VPD in all regions was selected in the regression model, indicating that VPD was a common factor affecting the ET₀ changes on the LP. However, for the model, only n and T_min were selected in the semihumid region; Pr was selected in the semiarid region; T, RH, and R_n were selected in the arid region; U₂ was selected in the semihumid and arid regions; and R_s was selected in the semihumid and semiarid regions, demonstrating that the above meteorological variables were regional impact factors.

To evaluate the relative impacts of dominant factors on ET₀ trends in different regions on the LP, a generalized linear model was used to calculate the relative contributions of the selected meteorological variables to the ET₀ changes (Table 5). The results revealed that the relative contributions of T, U₂, VPD, and R_s to the variability in ET₀ on the LP were 12.08%, 15.82%, 35.25%, and 36.85%, respectively. VPD and R_s may be considered the dominant meteorological factors mainly driving the changes in ET₀ on the LP, followed by U₂ and T. The increases in T and VPD led to an increase in ET₀, and the decreases in R_s and U₂ resulted in a reduction in ET₀. Finally, the counteracting effect of these major meteorological factors led to a decrease in ET₀ at a rate of –1.80 mm/10 yr from 1960 to 2017. Therefore, these changes explain why the ET₀ declined slightly from 1960–2017 on the LP. R_s and U₂ showed a negative effect, whereas VPD and T exerted positive influences. The combined effects of four meteorological variables resulted in negative effects that dominated the variability in ET₀. Moreover, among the four major meteorological factors on the LP, R_s, VPD and U₂ presented the highest relative contributions in the semihumid, semiarid, and arid regions, respectively.

Among the five major driving factors in the semihumid region, the decreased U₂ and R_s resulted in a negative trend of ET₀, while the reduction of n and the increases in T_min and VPD caused increases in ET₀. Finally, the changes in the five driving factors decreased ET₀ at a rate of 15.68 mm/10 yr. The results suggested that the decreased R_s was the main driving factor for the decreasing ET₀ in the semi-humid region from 1960–2017. This effect was intensified by the decreased U₂, although it was counteracted by the decreased n and increased T_min and VPD to certain degrees. In the semiarid region, the combined effect of VPD, R_s and Pr resulted in a downward trend of ET₀, with an increase in VPD representing a major driver with a contribution rate of 52.03%. Moreover, this effect was aggravated by the increased Pr and offset by the decreased R_s. In the arid region, ET₀ declined with decreasing U₂ and R_n. However, this negative effect was counterbalanced by the positive effects of the decreased RH and the increased T and VPD. Finally, the joint effect of the five meteorological variables augmented ET₀ at a rate of 3.57 mm/10 yr.

The temporal evolution of the relative contribution rates of the main meteorological variables to the ET₀ trend was determined by applying the generalized linear model for moving windows of 10 yr width (i.e., 1960–1969, 1961–1970, …, 2008–2017). The selected meteorological factors on the LP and in the three subregions were used as independent variables, and ET₀ was used as a dependent variable. The time series of the relative contributions of the major meteorological variables to the ET₀ changes is shown in Fig. 8.

Figure 8.  Temporal evolution of the contribution rates of the meteorological factors to ET₀ in the (a) semihumid region; (b) semiarid region; (c) arid region, and (d) the Loess Plateau. n, T, U₂, RH, Pr, Tmax, Tmin, VPD, R_s and R_n are same as in Table 4

According to Figs. 8a–8c, the relative contribution of VPD in the three subregions presented a fluctuating downtrend from 1969–2017 and played a dominant role in the semiarid region. R_s complemented the fluctuation trend of VPD in the semihumid and semiarid regions and jointly affected the ET₀ changes. The contribution rate of U₂ had stable trends in the semiarid and arid regions and played a relatively subordinate role in the arid region. The relative contributions of T_min and Pr showed stable upward trends, RH and R_n had fluctuating uptrends, and n and T exhibited stable trends.

Fig. 8d shows the evolution of the relative contributions of the major meteorological factors to ET₀ on the LP. The relative contribution of VPD exhibited an overall downward trend with two rapid decline phases, reaching a value of approximately 80% before 1974, which was followed by the first sharp decline to approximately 50%; then, it fluctuated in a small range, which was followed by the second rapid decline since 2012, reaching the valley bottom (approximately 4%) in 2014. The relative contribution of R_s showed high instability and fluctuated between approximately 4% and more than 50%, and it reached minimum values in 1991 and 2014 and showed two processes from rising to falling from 1969–2017. A stable upward trend was presented in the relative contribution of T, which became more significant after 2008. By 2017, the relative contribution of T exceeded VPD and R_s and reached 56.8%. The U₂ relative contribution fluctuated in the range of 0–29% from 1969–2017.

The specific contributions of the meteorological variables to ET₀ were further analyzed with the detrending method and the evaluating indicator R in the five subperiods of the entire LP (Fig. 9). From 1960–1972, U₂ increased at a rate of 0.18 (m/s)/10 yr (P < 0.1) (Table S1) and was a dominant factor resulting in the uptrend of ET₀ with an R value of 0.247. The upward trend of R_n played a relatively subordinate role in the ET₀ changes with an R value of 0.232. Whereas T_max had the largest negative influence on ET₀, with an R value of –0.129. The combined effects of the meteorological variables rendered ET₀ increase at a rate of 18.32 mm/10 yr (Table 2). From 1973–1990, U₂ was still a dominant factor influencing the ET₀ changes, with an R value of –0.653, but the significantly declining trend of U₂ led to a decline in ET₀ (–43.39 mm/10 yr, P < 0.01) (Table 2). R_n, VPD, RH, and R_s exerted certain impacts on the decrease in ET₀ from high to low, with R values of –0.481, –0.337, –0.196, and –0.165, respectively. After 1990, VPD played a dominant role in the ET₀ changes. From 1991–1999, VPD increased at a rate of 0.09 kPa/10 yr (P < 0.1) (Table S1), mainly driving the increase in ET₀, with an R value of 0.376. This effect was intensified by the decrease in RH (R = 0.290) and weakened by the decline in R_n (R = –0.166). However, VPD was the factor that had the greatest impact on the ET₀ changes, but its insignificant upward trend did not completely dominate the changes in ET₀ from 2000–2010. The positive effects of VPD and RH on ET₀ were offset by the negative effects of factors such as U₂, T_max, and R_s. Finally, the joint effect of these variables resulted in the reduction in ET₀ at a rate of –30.89 mm/10 yr (Table 2). From 2011–2017, the increases in VPD and U₂ increased ET₀, while the reductions in R_n and RH lessened ET₀. Together, their counteracting effect increased ET₀ at a rate of 79.99 mm/10 yr (P < 0.1) (Table 2).

Figure 9.  Absolute contribution of the meteorological variables to ET₀ in different time periods on the Loess Plateau. T, U₂, RH, T_max, T_min, VPD, R_s and R_n are same as in Table 4. T, U₂, RH, T_max, T_min, VPD, R_s and R_n are same as in Table 4

In this study, ET₀ was calculated with the daily meteorological data from 69 meteorological stations on the LP for the period from 1960–2017 by the FAO56-PM model. Uncertainties in the FAO56-PM ET₀ estimation and its contributing meteorological factors primarily relate to the model itself and the accuracy of the input data. The distribution of the meteorological stations is sparse and uneven, and the measurement dates of all the meteorological stations are inconsistent. The meteorological stations with as complete data as possible from 1960 to 2017 were selected. For the missing date value, we selected the daily data of 15 d before and after the missing dates and used the cubic spline interpolation method for supplementation. Moreover, uncertainty may come from spatial interpolation. The spline method was applied to obtain the spatial distribution of ET₀. The spline method is considered a reliable and robust method to interpolate hydro-meteorological variables (Price et al., 2000). However, elevation was not considered in the interpolation, which may introduce uncertainty. The FAO 56-PM ET₀ models need additional calculating equations to estimate R_n, which is not typically measured from meteorological stations. Previous research showed that the calculated R_n with the FAO 56-PM models was accurate enough to replace the measured R_n to calculate ET₀. The elevation is an input of the calculation of R_n (Allen et al., 1998). However, because of the complexity of the underlying surface and the uneven distribution of meteorological stations, other influencing factors cannot be ignored even when considering elevation, such as the land use change and interpolation of meteorological variables (Dong et al., 2020; 2021).

To verify the representativeness of our research results, we randomly selected nine meteorological stations and removed data in missing and discontinuous years to compare and verify the annual ET₀ with pan evaporation. The results showed that ET₀ had a positive linear relationship and large correlation coefficients (r) with pan evaporation (Fig. 10). Correlation coefficients were all above 0.7, and R² values were mostly above 0.6, with only two sites between 0.5 and 0.6. The result proved that there was a considerable linear correlation between ET₀ and pan evaporation. In addition, to verify the representativeness of meteorological stations, we randomly selected 90% of the stations and calculated the absolute contribution of meteorological variables to ET₀ in five subperiods with these stations and repeated three sample surveys. The results are shown in Fig. 11. Although the index R value of each meteorological variable was slightly different each time, the results of the three samples were consistent with the original conclusion. Therefore, despite the above uncertainties, this study confirmed the spatiotemporal variations in ET₀ and indicated the dominant factors contributing to ET₀.

Figure 10.  Relationships between the annual ET₀ calculated by the FAO56-PM model and pan evaporation on the Loess Plateau

Figure 11.  Absolute contribution of the meteorological variables to ET₀ of the original and three random sampling stations on the Loess Plateau in the subperiods. T, U₂, RH, T_max, T_min, VPD, R_s and R_n are same as in Table 4

Under the background of global climate change, downward trends in ET₀ have been found worldwide (Roderick and Farquhar, 2004; Irmak et al., 2012; Hosseinzadeh Talaee et al., 2014; Zheng and Wang, 2015), and more complex trends in ET₀ have been observed in recent years (Papaioannou et al., 2011; Shadmani et al., 2011). In this study, the annual ET₀ showed a slightly downward trend with obvious differences in various periods on the LP. Besides, the annual ET₀ decreased in the semi-humid region at a rate of –15.38 mm/10 yr and increased in the semiarid and arid regions. The annual ET₀ exhibited significant negative trends (P < 0.01) at 24 stations (approximately 34.8% of the total number of stations), which were largely located in the north and west of the semihumid and semiarid regions of the LP. The stations with an upward trend in the annual ET₀ were mainly distributed in the semiarid and arid regions of the LP, and 17 stations (approximately 24.6% of the total number of stations) reached the significance level (P < 0.01). The results are consistent with the conclusions described by Wang et al. (2012). Wang et al. (2012) found that the significant decreases in ET₀ were mostly located in the southeast, north and midwest of the Yellow River Basin, and the significant increases in ET₀ mainly occurred in the middle and southwestern corners.

Previous studies showed a change point of ET₀ in 1990, with ET₀ declining prior to 1990 and increasing after 1990 on the LP (Li et al., 2016). However, we found that ET₀ decreased significantly in all regions from 1960–1990, but ET₀ increased significantly in the semiarid region and decreased in the semi-humid and arid regions from 1991–2017. These findings indicated that the change trend of ET₀ was mainly attributed to the ET₀ changes in the semiarid region rather than the entire LP from 1991–2017 (Table S2). Thus, the trend of first increasing and then decreasing may not fully summarize the general trend of ET₀ in the LP in the past 58 years. The present study revealed that ET₀ showed a zigzag fluctuating trend of ‘increasing-decreasing-increasing-decreasing-increasing’, with change points at the year of 1972, 1990, 1999, and 2010. The results are consistent with the conclusions reported in some studies in similar regions in China that ET₀ or pan evaporation had abrupt changes in similar years of 1972, 1990, and 1999 (Liu et al., 2011; Li et al., 2016; Xing et al., 2016; She et al., 2017a). A new change point was detected in 2010 because the research time series was longer in our study.

Changes in meteorological variables caused by climate change and human activities are the essential influencing factors of ET₀ changes (Yin et al., 2010; Zhang et al., 2011). The analysis of the relative contribution rate revealed that the ET₀ changes on the LP were mainly affected by R_s (36.85%) and VPD (35.25%) from 1960–2017, followed by U₂ and T. The increase in VPD and T caused the rising trend of ET₀, while this effect was offset by the decrease in R_s and U₂, which are consistent with previous research in other regions. Qin et al. (2017) found that the reduction in wind speed was the key factor resulting in the downward trend of annual ET₀ from 1961–1987, while the increase in VPD and the decrease in relative humidity were the primary reasons for the augmentation in annual ET₀ from 1988 to 2012 in the Qinhuai River Basin in the southern China. She et al. (2017a) investigated the characteristics and driving factors of ET₀ changes in the middle of the Yellow River in China and revealed that solar radiation and wind speed controlled the changes in ET₀ at most stations and in most subperiods. Maček et al. (2018) analyzed the impact of meteorological variables on the ET₀ trends in different climate types of weather stations in Slovenia, Europe. They denoted that the relative impact of wind speed on ET₀ is much smaller than the relative contribution of solar radiation, and the latter generally controls the change in reference evapotranspiration, followed by temperature and VPD.

The major driving factors of ET₀ are diverse in different subregions of the LP. Solar radiation, VPD, and wind speed are the key factors affecting the ET₀ changes in the semihumid, semiarid, and arid regions, respectively. The results of other reports provided support for our conclusions. Some studies indicated that the reduction in surface solar radiation became apparent in many observational records up to 1990, a phenomenon called global dimming, which may substantially affect surface temperature, the hydrological cycle, and ecosystems and cause a decline in the essential energy of evaporation (Roderick and Farquhar, 2002; Wild et al., 2005). However, some indications show that dimming did not persist into the 1990s and the amount of sunlight at the surface has increased, which may profoundly influence the surface climate (Wild et al., 2005). Such a change trend in solar radiation has also been confirmed in China (Che, 2005). Some studies stated that the phenomenon of ‘from dimming to brightening’ could be partially caused by the increased aerosol loading from pollutant emissions and the changes in the mix of fuels and consumption technologies since the 1980s in China (Qian et al., 2007b). Moreover, the absorption of R_s could be affected by water vapor in the atmosphere, and R_s reaching the land surface may be potentially influenced by variations in RH via the hygroscopic effect of aerosols and absorption at solar wavelengths (Qian et al., 2007b). This finding may partly explain why solar radiation has become the dominant factor affecting ET₀ in the semi-humid region of the LP. McVicar et al. (2012) evaluated the trend of global wind speed and its effect on the changes in ET₀ and revealed a worldwide reduction in the observed wind speed at the land surface (called ‘global stilling’). They also found that ET₀ was extremely sensitive to wind speed and ‘global stilling’ was considered the primary factor in the decline of ET₀. Studies have shown that the increasing land surface roughness caused by vegetation growth and rapid urbanization may explain this phenomenon to a certain extent (She et al., 2017a). The wind speed has dropped at a rate of –0.17 (m/s)/10 yr in recent decades in China (Chen et al., 2013), although the influence of wind speed on ET₀ was gradually decreasing, and the effects of other meteorological factors were becoming more important (Liu et al., 2011).

Wind speed was the major driving factor for the changes in ET₀ on the LP prior to 1990, which is supported by the conclusions of some studies on the same region. Li et al. (2016) believed that the change in wind speed was the main factor causing the downward trend of ET₀ and solar radiation played a relatively subordinate role between 1960 and 1990. Our study found that the changing trend of wind speed from 1960–1972 and 1973–1990 changed, which had the opposite effects on the changing trend of ET₀. The upward trend of ET₀ was primarily attributed to the increase in wind speed on the LP from 1960–1972, and ET₀ showed a downward trend because of the significant reduction in wind speed from 1973–1990 (Fig. 8). The upward and downward trends of net radiation from 1960–1972 and 1973–1990 became the second leading factors for the increase/decrease in ET₀, while the index R of solar radiation in the two subperiods was only 0.046 and –0.165 and thus did not play a dominant role (Fig. 8). After 1990, VPD was the dominant factor for ET₀ changes on the LP, which is inconsistent with previous studies in this region. Li et al. (2016) suggested that the rapidly rising temperature drove the changes in ET₀ after 1990 because previous research mainly focused on the influence of temperature and generally ignored the influence of VPD changes. VPD is the difference between the saturation vapor pressure and actual vapor pressure (Allen et al., 1998). Saturation vapor pressure is a function of the maximum temperature and minimum temperature, and actual vapor pressure is a function of the relative humidity. Therefore, VPD is a comprehensive function of temperature and relative humidity. Generally, higher air temperature and lower relative humidity cause higher VPD. It is difficult to completely distinguish the effects of relative humidity and temperature on ET₀ change trends on the LP because the changes in relative humidity are closely related to the warming of the study area. The increase in temperature, the intensification of drought and the decrease in soil moisture have led to a decline in atmospheric water content as well as a significant reduction in relative humidity. Relative humidity controls the loss of air moisture from plant leaves by increasing VPD (Liu and Feng, 2012), and the air in the study area becomes drier, which results in higher VPD and ET₀. Therefore, our research shows that VPD provides a more powerful explanation for the changes in ET₀ compared with the individual influences of temperature and relative humidity on ET₀. This finding is consistent with the research conclusion of Qin et al. (2017). The upward trend of ET₀ from 1988–2012 was mainly affected by decreased relative humidity and increased VPD; that is, the gradually dry environment of the Qinhuai River Basin was the primary reason for the increase in ET₀. Overall, VPD increased significantly throughout the study period (Table S1), which is consistent with the results of a global study (Bais et al., 2011).

The contribution of the meteorological variables to the ET₀ change trends has been found to have strong temporal and spatial variability on the LP, which may be caused by the strong gradients of the meteorological variables and the variation range of the variables in the spatial pattern. These variables have strong spatial changes due to climatic factors, topographic features and land cover/use on the LP. In the past few decades, the LP has undergone engineering measures, such as the Grain for Green Project, grassland reconstruction, terraces and dams. Although vegetation coverage has increased, large-scale vegetation restoration has also caused excessive consumption of soil water and affected the evolution of local meteorological variables, which may further explain the complex temporal and spatial patterns observed in the ET₀ changes.

In this study, daily data from 69 meteorological stations on the Loess Plateau (LP) from 1960 to 2017 were applied to calculate the reference evapotranspiration (ET₀) with the FAO-56 Penman-Monteith method. The spatiotemporal pattern evolution of ET₀ and the contribution of related meteorological variables to the ET₀ changes were analyzed using the Mann-Kendall test, segmented regression model, simple linear regression, wavelet analysis, and generalized linear model. Compared with previous studies on the ET₀ change trends and its attribution on the LP, this study emphasized the importance of the change points in the long-term ET₀ and quantitatively calculated the contributions to the changes by the detrending method.

(1) In the past 58 yr, the annual ET₀ on the LP showed strong gradients with an initial decrease and then increase from northwest to southeast and had a short cycle of 6 yr and a medium cycle of 11–13 yr but did not show a long-term cycle.

(2) From 1960–2017, R_s and VPD were the primary factors driving the annual ET₀ changes on the LP. With complex interactions of the decreases in R_s and U₂ and the increases in T and VPD, ET₀ declined slightly at a rate of –1.80 mm/10 yr. R_s, VPD, and U₂ were the dominant factors contributing to the ET₀ changes and resulted in downward, upward, and upward trends in annual ET₀ in the semihumid, semiarid, and arid regions, respectively.

(3) Four change points at the year of 1972, 1990, 1999, and 2010 were detected by the segmented regression model, which can divide the study period into five subperiods and lead to a zigzag change trend of ‘increasing-decreasing-increasing-decreasing-increasing’ in ET₀ in the three subregions. Furthermore, the dominant meteorological factors affecting ET₀ underwent significant changes. The ET₀ changes were mainly attributed to U₂ before 1990 and VPD after 1990. In general, the variations in ET₀ were influenced by the combined effects of the meteorological driving factors, and U₂ and VPD were the two major contributing factors that control the ET₀ variations from 1960–1990 and from 1991–2017, respectively.