Attribution analysis of the spatial variations in potential evapotranspiration on the Loess Plateau of China by a total differential equation

As one of the major components of the hydrologic cycle, evapotranspiration (ET) drives energy and water exchanges among the hydrosphere, atmosphere, and biosphere (Wang et al. 2007). An accurate estimate of ET is important for improving the efficient use of water resources and for planning future water resource management (Gavilan & Castillo-Llanque 2009). However, direct measurement of ET is time-consuming and expensive. Thus, potential evapotranspiration (ET₀) has been widely used in water resource studies and management practices (Li et al. 2014) because it can be directly estimated by meteorological data (Allen et al. 1998).

In recent decades, global air temperature has become increasingly significant. However, ET₀ and pan evaporation have decreased worldwide over the past decades, known as the ‘evaporation paradox’ phenomenon (Brutsaert & Parlange 1998). Subsequently, attribution analysis has been widely conducted to understand this phenomenon. Three main approaches including the linear stepwise multivariate regression, detrended method, and sensitivity analysis method are commonly used. The linear stepwise multivariate regression method generally considers the major climatic factors as independent variables; a higher correction associated with ET₀ for a certain variable suggests that the variable has a larger impact on the ET₀ variations (Liu et al. 2004, 2015; Shan et al. 2015). The steps of the detrended method for quantifying the contributions of a certain climatic variable to the trend of ET₀ include the following: (1) removing the trend from a time-series climatic variable to make it stationary; (2) recalculating ET₀ by using this detrended variable and other climatic variables; (3) the differences between the original ET₀ values and the recalculated ones indicate the influence of this variable imposed upon ET₀ (Xu et al. 2006; Liu & Yang 2010). Compared with the former two methods, the total differential method has solid theories, because it calculates the contributions of major variables based on the total differential equation. Specifically, the ET₀ change induced by a certain variable is evaluated by multiplying the long-term trend of the variable with its partial derivative. Further, many studies have demonstrated that the attribution errors induced by this method are relatively small and the sum of the contributions of each variable fit well with the actual trend of ET₀ (Roderick et al. 2007; Zheng et al. 2009; Liu et al. 2011; Ning et al. 2016).

However, the above three approaches have rarely been applied to conduct quantitative attribution analysis of the spatial ET₀ variations. Current studies are more likely to focus on qualitative analysis to solve the most related climatic variables that contribute to the spatial distribution of ET₀; for example, only by qualitatively analyzing the spatial patterns of ET₀ and related climatic variables, did Li et al. (2012) detect the main cause of the lowest ET₀ values in the southwest region and the highest ET₀ values in the northwest region. Furthermore, methods such as clustering or zonality analysis have also been used to assist the attribution study over the spatial variations of ET₀ and other climatic variables. Based on clustering analysis, Guangdong Province in China was divided into four parts by He et al. (2015) to calculate the correlation between pan evaporation and certain climatic variables in order to detect the underlying driving factors that cause pan evaporation changes in subregions. The same method was also used to study the spatiotemporal variations of pan evaporation over all of China by Zhang et al. (2015). Even though clustering analysis can be regarded as a quantitative method operated on the spatial variations of ET₀, it is still a temporal analysis that essentially explores the attributes of ET₀ changes in subregions. The zonality analysis method, building the multiple regression equations for climatic variables with longitude, latitude, and altitude information, can indicate how the climatic variables vary spatially with these geographical factors (Li et al. 2016a; Sun & Zhang 2016). Nevertheless, the problem with this method is that it ignores the interaction among the climatic variables.

Therefore, in this study, we extended the total differential method based on the FAO 56 Penman–Monteith formula to conduct an attribution analysis of the spatial ET₀ variations over the Loess Plateau, China. In addition, the attribution results were compared with those of the multiple regression analysis. Three main directions with distinct ET₀ gradients were selected to operate the attribution analysis.

The Loess Plateau is located in the upper and middle reaches of the Yellow River in North China (Figure 1), covering a total area of 6.4 × 10⁵ km². Most areas are dominated by the semi-arid and subhumid climate, with decreased precipitation along the southeast–northwest direction, ranging from 200 to 750 mm. Characterized by highly erodible wind-deposited loess soil, sparse vegetation, but heavy rainstorms in the summer, the plateau faces the severe problems of erosion and sedimentation. Several soil and water conservation practices have been implemented since the 1950s. Other human activities, such as coal mining and urbanization, have also seriously altered the local landscapes and hydrological cycle.

Daily meteorological data composed of daily maximum and minimum temperature at an elevation of 2 m, atmospheric pressure, wind speed at an elevation of 10 m, mean relative humidity, and sunshine duration were used to calculate the daily ET₀. They were obtained from the China Meteorological Administration (http://data.cma.cn/); the records were collected from 96 stations during the period of 1960–2013. There are a few missing data points that were replenished based on the regression relationship with their neighboring station records. DEM (digital elevation model) data were provided by Geospatial Data Cloud of China.

The root mean square error (RMSE) and the relative error (RE) were employed to assess the performance of the cokriging interpolation method (Table 1). The results showed that the cokriging method that considers the elevation information can improve the accuracy of the interpolation for all variables.

Considering the spatial distribution of ET₀, three transects along different directions were chosen to conduct the attribution analysis of the spatial variation of ET₀. The ET₀ samples were taken at an interval distance of 10 km along each transect. Then, the mean ET₀ at each point with a 5-km buffering area was extracted to represent the spatial ET₀ variations along the transects. Similar operations were performed on the four variables of U₂, R_s, T_mean, and e_a. Subsequently, the spatial variation sequences for ET₀ and those for the four climatic variables were obtained along the three transects.

To prove the attribution results of the total differential method, a multiple stepwise regression method was also applied to identify the primary and the leading climatic variables. Using SPSS, a multiple regression analysis was carried out by considering the spatial series of ET₀ as a dependent variable and the four climatic variables (U₂, R_s, T_mean, and e_a) as independent variables. Specifically, the standard regression coefficient served as the basis to determine the contribution of each climatic variable to the ET₀ variation. The larger the standard regression coefficient is, the larger the contribution.

Figure 2 shows the spatial distribution of mean annual ET₀ over the Loess Plateau for the period 1960–2013. The mean annual ET₀ was 987.3 mm, ranging from 700 to 1,226 mm. ET₀ distinctly decreased in three directions, i.e., from the northwest to east, south, and southwest of the Loess Plateau. To perform the quantitative attribution analysis for the spatial variations in ET₀ over the Loess Plateau, three transects were set along these three directions, namely, northwest–east, northwest–south, and northwest–southwest. High ET₀ values (more than 1,050 mm/a) were found in the Ordos Plateau and the Ningxia Plain to the northwest and in Sanmenxia to the southeast. The lowest ET₀ values (less than 900 mm/a) were discovered in the west of Liupan Mountain, the southwest and the north of Wutai Mountain, and the northeast of Lvliang Mountain.

The spatial distributions of U₂, R_s, T_mean, and e_a over the Loess Plateau are shown in Figure 3. The wind speed was not substantially varied in space. Except for the regions to the west and to the east, the U₂ values in most of the regions are in the range of 1.5–2.5 m/s. Solar radiation decreased from northwest to southeast, with the highest R_s values (more than 16.5 MJ /m²) found in the northwest and the lowest values (less than 15.0 MJ /m²) found in the southeast. The spatial distribution of the mean temperature was generally in line with the changes in the actual vapor pressure, decreasing from southeast to northwest.

The mean annual sensitivity coefficients for wind speed (S_U₂), solar radiation (S_R_s), mean temperature (S_T_mean), and actual vapor pressure (S_e_a) were calculated for each station, and then they were interpolated to show the spatial pattern (Figure 4). S_U₂ increased from southeast to northwest in the range of 0.13–0.28, which means a 10% increase in U₂ would result in a larger change in ET₀ in the northwest than in the southeast. Conversely, S_R_s decreased from southeast to northwest in the range of 0.25–0.50. The highest values of S_T_mean were found in the southeast (more than 0.6), and then it decreased from southeast to southwest, northwest, and north. The absolute value of S_e_a decreased from south to north ranging from −0.75 to −0.34.

The mean annual sensitivity coefficients for the whole Loess Plateau were the largest for e_a (−0.51), intermediate for R_s and T_mean (0.38 and 0.38), and the smallest for U₂ (0.20). These results indicated that a 10% decrease in R_s, T_mean, and U₂ would result in a 3.8%, 3.8%, and 2.0% decrease in ET₀, respectively, while a 10% decrease in e_a would result in a 5.1% increase in ET₀.

Figure 2 shows that there was a substantial decrease in ET₀ from the Ertuoke Banner in the northwest of the Loess Plateau to Fanzhi County in the east. Figure 5(a) presents how ET₀ varied with different sample points along the northwest–east transect. Except for the relatively high values captured around Shuozhou County in Shanxi Province, ET₀ exhibited a clear longitudinal zonality with the change rate of −4.9 mm/10 km from northwest to east. In the same direction, U₂ and e_a increased and R_s and T_mean decreased. Except for the U_2, trends in the other three climatic variables all were significant (p < 0.01) (Table 2).

Equations (5)–(9) were used to quantify the contribution associated with the four individual factors to the spatial variation in ET₀. The results implied that the downward trend in the mean annual ET₀ along the direction of 39°N was potentially contributed by a decrease in R_s and an increase in e_a, with the contribution values identified as −2.65 and −1.32 mm/10 km, equal to 52.6% and 26.3%, respectively. Lower impacts were calculated for U₂ (5.3%) and T_mean (15.8%) on the decreased ET₀ along the northwest–east transect. The attribution results of the multiple regression analysis showed that the contribution of R_s to the ET variation along this transect was largest with a stand regression coefficient of 0.55 (p < 0.01), followed by e_a, T_mean, and U₂, which agrees with the results of the total differential method (Table 3). However, the multiple regression analysis could not calculate the specific contribution rate of each variable.

The ET₀ variation from northwest to south across the Loess Plateau is shown in Figure 5(b). This transect is at the longitude of 108°E. ET₀ decreased overall from northwest to south at the rate of −4.5 mm/10 km. Along this transect, U₂ and R_s decreased significantly (p < 0.01), and their contributions to the decreased ET₀ were calculated as −1.50 and 0.13 mm/10 km, equal to the relative contributions of 34.4% and −3.04%, respectively. In contrast, T_mean and e_a showed an increasing trend (significance test: p < 0.01); their contributions to the decreased ET₀ were identified as 7.97 and −10.95 mm/10 km, which are equivalent to a relative contribution of −183.2% and 251.9%, respectively. In summary, even though an increasing T_mean could result in an increase in ET₀, such a trend might be offset by the function of e_a and U₂. The attribution results of the multiple regression also proved that e_a was the dominant factor that influenced the ET₀ variation along the northwest–south transect.

It is worth noting that ET₀ strongly decreased from Dingbian County to Wuqi County in the Shaanxi Province at the rate of 15.8 mm/10 km along the northwest–south transect. The decreasing U₂ made the maximum contribution to the decreased ET₀ (58.8%), followed by increasing T_mean and e_a (23.2% and 15.7%), and the contribution of decreasing R_s was minimal.

ET₀ also decreased significantly from northwest to southwest at the rate of −4.0 mm/10 km (Figure 5(c)). Except for e_a, the other three variables all exhibited a decreasing trend, while the four climatic variables all decreased ET₀ in this direction. The maximum relative contribution to the decreased ET₀ along this transect was U₂, with the values computed as 46.6%, followed by R_s (37.7%) and e_a (15.4%), compared with the minimum percentage of T_mean (Table 4). The order of importance of the climatic variables to the ET₀ variation along this transect detected by the sensitivity method agrees with that by the multiple regression method.

Even though ET₀ showed an overall decreasing trend along this direction, it displayed the following regional differences: (1) from Ertuoke Banner to Lingwu County in Ningxia Province, ET₀ decreased only at the rate of −2.9 mm/10 km; (2) from Lingwu to Zhongning County in Ningxia Province, ET₀ increased up to 1,170 mm/a by a rate of 10.6 mm/10 km; (3) from Zhongning to Dongxiang County in Ningxia Province, ET₀ decreased at a rate of −9.4 mm/10 km to only 870 mm per year. U₂ seems to be the dominant factor for the decreasing trend in ET₀ in the first two regions, with 194.6% and 84.9% in the relative contribution; however, U₂ decreased in the first region but increased in the second region. In the third region, the impacts of U₂ on ET₀ were degraded to only 24.5%; ET₀ seems to be strongly associated with R_s (64.6%).

Due to little precipitation and cloud coverage in the northwest of the Loess Plateau, where the region is dominated by a semi-arid climate (Qian 1991), the actual vapor pressure tends to be low but high in solar radiation (see Figure 3). This collectively contributes to the high ET₀ in this region. However, the southeastern part of the Loess Plateau is affected by a semi-humid continental monsoon climate, which is jointly controlled by summer monsoon and westerlies (Liang et al. 2015; Yan et al. 2016); thus, the local precipitation is richer than other regions of the Loess Plateau. As a consequence, the high actual vapor pressure and low solar radiation in the southeast result in relatively low values of ET₀. Combining the relatively low elevation in the southeast region with a valley terrain (Wu et al. 1982), the temperature in this area is high, and it makes a large positive contribution to ET₀. Furthermore, there are several abrupt points of ET₀ on the three transects in our study (see Figure 4); they are highly related to the variation in terrain. For example, from the Ordos Plateau in the northwest to the loess hilly-gully region in the south (Figure 1), the terrain became more complex and the land surface elevation roughly decreased, which may explain the decrease in wind speed in the loess hilly-gully region (i.e., from Dingbian County to Wuqi County in Figure 5(b)), and finally resulted in a rapid decrease in ET₀. To further quantify the impacts of topography on the spatial change in ET₀, the correlation between the mean annual ET₀ and elevation for 96 stations was calculated. The calculation results showed that the determining coefficient was 0.25 (p < 0.01), indicating that regional variation in topography would significantly impact the spatial pattern of ET₀ (Figure 6(a)).

Moreover, the change in vegetation cover is another potential impact on ET₀. It was reported that the vegetation cover has been effectively improved in the loess hilly-gully region due to the implementation of the ‘Grain to Green Project’ since 1999 (Zhang et al. 2013; Sun et al. 2015), which may increase the surface roughness, decrease the wind speed (Vautard et al. 2010; Cowie et al. 2013) and, finally, impact ET₀. With the Global Inventory Modeling Studies 3 g database (GIMMS3 g), the mean annual NDVI was calculated from 1982 to 2012, which can reflect the vegetation condition of the whole Loess Plateau. Then, the NDVI value for each station was extracted and its relationship with ET₀ was evaluated. The results showed that ET₀ was significantly related to NDVI (p < 0.05), suggesting that the improvement of the vegetation condition would result in the ET₀ change (Figure 6(b)). However, the revegetation project was not carried out for the whole Loess Plateau, and it mainly focused on the hilly-gully region (Xin et al. 2008). Thus, the heterogeneity of human interference would influence the spatial pattern of ET₀.

There is no denying that the spatial interpolation accuracy of ET₀ and the four related climatic variables could impact the final imposed uncertainty of our attribution analysis. Here, we will try to discuss this uncertainty. The spatial distribution of these variables generally follows patterns similar to those of previous studies (Li et al. 2012; Li et al. 2016b). Thus, it can be determined that the interpolation error in this study is acceptable and it will not influence the global spatial distribution characteristics of climatic variables. Following the same steps of interpolating mean annual values of ET₀ and the four climatic variables, the annual values of these five variables from 1961 to 2012 were first interpolated. Then, their values at each corresponding station were extracted. Finally, the attribution analysis of the extracted ET₀ change was conducted and the results were compared with that of the actual ET₀ for each station (Figure 7). The results showed that the contributions of the interpolated U₂, R_s, T_mean, and e_a variables to ET₀ were similar to those of the four actual variables to ET₀, with a high determination coefficient R² (0.85, 0.88, 0.91, and 0.80). Thus, it can be concluded that the influence of the spatial interpolation error on the attribution result was acceptable.

Based on the FAO 56 Penman–Monteith formula, this study used the total differential equation method to quantify the contributions of each impact factor to the spatial ET₀ variation across the Loess Plateau. Although this method is effective in presenting the contribution of a certain climatic variable to the spatial variation in ET₀, errors exist between the calculated and observed total ET₀ changes. The errors exist in all three transects (see Ɛ values in Tables 1–3). The maximum Ɛ appeared in the direction of 39°N, accounting for 2.4% of the observed ET₀ change. Similar phenomena have been detected in the attribution analysis for the temporal variation of ET₀ (Zheng et al. 2009; Liu et al. 2011, 2013; Feng et al. 2014).

Except for the impact of the interpolation accuracy, several other aspects can also potentially contribute to the errors. First, only a limited number of climatic variables (four in this study) are considered in such kinds of studies; therefore, the ET₀ changes calculated from the limited climatic variables cannot represent those from the climate (Zheng et al. 2009). In addition, the errors could come from the assumption of the total differential equation that the climatic variables are independent from each other, which is not true in reality; for example, the increasing temperature would result in a decrease in actual vapor pressure (Liu et al. 2011; Liu & Zhang 2013). Furthermore, the method used for contribution analysis is actually a Taylor expansion. Considering only the first-order approximation to evaluate their contributions would underestimate the calculated total changes in ET₀.

In this study, the total differential equation was applied to analyze the contribution of four climatic factors to spatial variation in ET₀ over the Loess Plateau. First, the annual ET₀ records collected from 96 stations in the study area during 1960–2013 were calculated using the FAO 56 Penman–Monteith equation. Second, the annual sensitivity coefficients of ET₀ to U₂, R_s, T_mean, and e_a were calculated for each station. Then, the mean annual ET₀, four climatic variables and sensitivity coefficients were mapped using the cokriging interpolation method to show their spatial patterns. Three transects with distinguishable ET₀ gradients were selected and sampled with an interval of 10 km for maps of ET₀ and the four climatic variables. Finally, using distance as an independent variable, the contributions of each climatic variable to the spatial ET₀ variations along the three transects were assessed by the total differential equation method. In addition, the attribution results were compared with those of the multiple regression analysis. The analysis results showed the following: the mean annual ET₀ roughly decreased from northwest to east, south, and southwest over the Loess Plateau. The mean annual sensitivity coefficients for the whole Loess Plateau were largest for e_a, intermediate for R_s and T_mean, and smallest for U₂. The dominant factors that caused ET₀ to decrease from northwest to east and to southwest were the decreasing solar radiation and increasing actual vapor pressure, respectively. As another key contributor, the decreasing wind speed induced a decreasing trend in ET₀ from northwest to southwest. Except for the three directions given in our study, the total differential equation method can also be used for any direction. Compared with the result of the multiple regression analysis, the total differential method can not only obtain the dominant factor but also accurately evaluate the contribution rate of each climatic factor.

Current studies are more likely to focus on qualitative analysis to determine the most related climatic variables that contribute to the spatial distribution of ET₀. This study provides a quantitative method to conduct the attribution analysis of the spatial ET₀ variations. In addition, this method can also conduct an attribution analysis of the spatial variation in other climatic variables, which should have a clear functional relation with its impact factors.

This study was supported by the National Key Research and Development Program of China (No. 2016YFC0501602), the Opening Fund of State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau (A314021402-1804), the National Natural Science Foundation of China (No. 41571036), and the Public Welfare Industry (Meteorological) Research Project of China (No. GYHY201506001).