## Abstract

A total differential equation was proposed to assess the driving factors for the spatial variations in potential evapotranspiration (*ET*_{0}). Using China's Loess Plateau as an example study area, three transects with distinct *ET*_{0} gradients in space, i.e., northwest–east, northwest–south and northwest–southwest, were chosen to sample spatially varied *ET*_{0} and four climatic variables (solar radiation, actual vapor pressure, wind speed, and mean temperature) at an interval of 10 km. Considered an independent variable, the distance was differentiated to quantify the contribution of each climatic variable to the spatial *ET*_{0} variations along each transect. A significant decrease in solar radiation and an increase in actual vapor pressure were identified as the dominant impact factors that led to a decreased *ET*_{0} in the northwest–east and northwest–south directions, respectively. As another key contributor, the decreasing wind speed induced a decreasing trend in *ET*_{0} from northwest to southwest. The above results implied that the dominant factor(s) for the spatial variations in *ET*_{0} differed among the regions. Therefore, the total differential equation is a powerful approach to determine the driving factors and to quantify their individual contribution to the spatial variations in *ET*_{0}.

## INTRODUCTION

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*_{0}) 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*_{0} 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*_{0} for a certain variable suggests that the variable has a larger impact on the *ET*_{0} 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*_{0} include the following: (1) removing the trend from a time-series climatic variable to make it stationary; (2) recalculating *ET*_{0} by using this detrended variable and other climatic variables; (3) the differences between the original *ET*_{0} values and the recalculated ones indicate the influence of this variable imposed upon *ET*_{0} (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*_{0} 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*_{0} (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*_{0} 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*_{0}; for example, only by qualitatively analyzing the spatial patterns of *ET*_{0} and related climatic variables, did Li *et al.* (2012) detect the main cause of the lowest *ET*_{0} values in the southwest region and the highest *ET*_{0} 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*_{0} 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*_{0}, it is still a temporal analysis that essentially explores the attributes of *ET*_{0} 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*_{0} 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*_{0} gradients were selected to operate the attribution analysis.

## DATA AND METHODS

### Study area

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^{5} km^{2}. 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.

### Data collection

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*_{0}. 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.

### Calculation of potential evapotranspiration

^{−1}). Reference evapotranspiration is a special case of potential evapotranspiration, because it also describes the maximum evaporative capacity from the surface, and even the surface is defined as a hypothetical grass reference crop with specific characteristics. Thus, reference evapotranspiration was used to represent the potential evapotranspiration in this study, and its specific calculation was as follows: where

*R*

_{n}is the net radiation at the crop surface (MJ m

^{−2}d

^{−1}); G is soil heat flux density (MJ m

^{−2}d

^{−1});

*T*

_{mean}is the mean daily air temperature measured at an elevation of 2 m (°C), which is defined as the mean of the daily maximum (

*T*

_{max}) and minimum temperatures (

*T*

_{min});

*U*

_{2}is the wind speed at an elevation of 2 m (m s

^{−1}); e

_{s}is the saturation vapor pressure (kPa);

*e*

_{a}is the actual vapor pressure (kPa); Δ is the slope of the vapor pressure curve (kPa °C

^{−1}); and is the psychrometric constant (kPa °C

^{−1}).

### Spatial interpolation

*ET*

_{0}of each station was calculated. Then, the cokriging method was chosen to spatially interpolate the mean annual

*ET*

_{0},

*U*

_{2},

*R*

_{s},

*T*

_{mean}, and

*e*

_{a}at a resolution of 1 km. Elevation information extracted from the DEM data was used as a correction factor. The correlation between climatic variables and elevation were described by the crossover-mutation function as: where

*z*(

*x*) and

*Z*are the climate data for the unknown point and station

_{ui}*i*, respectively;

*y*(

*x*) is the elevation data;

*n*is the number of stations;

*m*

_{y}and

*m*

_{z}are the mean elevation and climatic variable values, respectively; and and are the weights.

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.

ET_{0} (mm) | U_{2} (m/s) | R_{s} (MJ/m^{2}) | T_{mean} (°C) | e_{a} (kPa) | ||
---|---|---|---|---|---|---|

Not considering elevation | RMSE | 9.79 | 0.22 | 0.35 | 1.16 | 0.012 |

RE (%) | 10.3 | 5.5 | 1.7 | 12.1 | 4.2 | |

Considering elevation | RMSE | 5.82 | 0.18 | 0.28 | 0.33 | 0.007 |

RE (%) | 8.2 | 4.9 | 1.3 | 3.6 | 2.8 |

ET_{0} (mm) | U_{2} (m/s) | R_{s} (MJ/m^{2}) | T_{mean} (°C) | e_{a} (kPa) | ||
---|---|---|---|---|---|---|

Not considering elevation | RMSE | 9.79 | 0.22 | 0.35 | 1.16 | 0.012 |

RE (%) | 10.3 | 5.5 | 1.7 | 12.1 | 4.2 | |

Considering elevation | RMSE | 5.82 | 0.18 | 0.28 | 0.33 | 0.007 |

RE (%) | 8.2 | 4.9 | 1.3 | 3.6 | 2.8 |

### Transect sampling for spatial *ET*_{0} variations

Considering the spatial distribution of *ET*_{0}, three transects along different directions were chosen to conduct the attribution analysis of the spatial variation of *ET*_{0}. The *ET*_{0} samples were taken at an interval distance of 10 km along each transect. Then, the mean *ET*_{0} at each point with a 5-km buffering area was extracted to represent the spatial *ET*_{0} variations along the transects. Similar operations were performed on the four variables of *U*_{2}, *R*_{s}, *T*_{mean}, and *e*_{a}. Subsequently, the spatial variation sequences for *ET*_{0} and those for the four climatic variables were obtained along the three transects.

### Attribution analysis for the spatial *ET*_{0} variations

#### Total differential method

*ET*

_{0}for a certain transect with a unit of mm/10 km; C_(

*U*

_{2}), C_(

*R*

_{s}), C_(

*T*

_{mean}), and C_(

*e*

_{a}) are the contributions of

*U*

_{2},

*R*

_{s},

*T*

_{mean}, and

*e*

_{a}to the spatial

*ET*

_{0}changes, respectively. Furthermore, the relative contributions of each factor to the spatial change in

*ET*

_{0}for each transect can be calculated as:

### Multiple regression analysis

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*_{0} as a dependent variable and the four climatic variables (*U*_{2}, *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*_{0} variation. The larger the standard regression coefficient is, the larger the contribution.

## RESULTS

### Spatial variations in the mean annual *ET*_{0} and in the other four climatic variables

Figure 2 shows the spatial distribution of mean annual *ET*_{0} over the Loess Plateau for the period 1960–2013. The mean annual *ET*_{0} was 987.3 mm, ranging from 700 to 1,226 mm. *ET*_{0} 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*_{0} over the Loess Plateau, three transects were set along these three directions, namely, northwest–east, northwest–south, and northwest–southwest. High *ET*_{0} 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*_{0} 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*_{2}, *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*_{2} 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^{2}) found in the northwest and the lowest values (less than 15.0 MJ /m^{2}) 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.

### Spatial changes in the sensitivity of *ET*_{0} to the four climatic variables

The mean annual sensitivity coefficients for wind speed (*S*_*U*_{2}), 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*_{2} increased from southeast to northwest in the range of 0.13–0.28, which means a 10% increase in *U*_{2} would result in a larger change in *ET*_{0} 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*_{2} (0.20). These results indicated that a 10% decrease in *R*_{s}, *T*_{mean}, and *U*_{2} would result in a 3.8%, 3.8%, and 2.0% decrease in *ET*_{0}, respectively, while a 10% decrease in *e*_{a} would result in a 5.1% increase in *ET*_{0}.

*ET*_{0} attribution analysis along three transects

#### Along the northwest–east transect

Figure 2 shows that there was a substantial decrease in *ET*_{0} from the Ertuoke Banner in the northwest of the Loess Plateau to Fanzhi County in the east. Figure 5(a) presents how *ET*_{0} varied with different sample points along the northwest–east transect. Except for the relatively high values captured around Shuozhou County in Shanxi Province, *ET*_{0} exhibited a clear longitudinal zonality with the change rate of −4.9 mm/10 km from northwest to east. In the same direction, *U*_{2} 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).

ET_{0} | U_{2} | R_{s} | T_{mean} | e_{a} | Ɛ(%) | ||
---|---|---|---|---|---|---|---|

Northwest to east | Slope | −4.9** | 0.006 | −0.031** | −0.033* | 0.002** | |

Contributions (mm/10 km) | −0.27 | −2.65 | −0.80 | −1.32 | −0.12 | ||

Relative contributions (%) | 5.3 | 52.6 | 15.8 | 26.3 | −2.4 | ||

Northwest to south | Slope | −4.5** | −0.012** | −0.039** | 0.072** | 0.008** | |

Contributions (mm/10 km) | −1.50 | 0.13 | 7.97 | −10.95 | 0.10 | ||

Relative contributions (%) | 34.4 | −3.04 | −183.2 | 251.9 | 2.3 |

ET_{0} | U_{2} | R_{s} | T_{mean} | e_{a} | Ɛ(%) | ||
---|---|---|---|---|---|---|---|

Northwest to east | Slope | −4.9** | 0.006 | −0.031** | −0.033* | 0.002** | |

Contributions (mm/10 km) | −0.27 | −2.65 | −0.80 | −1.32 | −0.12 | ||

Relative contributions (%) | 5.3 | 52.6 | 15.8 | 26.3 | −2.4 | ||

Northwest to south | Slope | −4.5** | −0.012** | −0.039** | 0.072** | 0.008** | |

Contributions (mm/10 km) | −1.50 | 0.13 | 7.97 | −10.95 | 0.10 | ||

Relative contributions (%) | 34.4 | −3.04 | −183.2 | 251.9 | 2.3 |

*Note:* * and **suggest that the trend of a given variable is significant at the level of *p**=* 0.05 and *p**=* 0.01 by the F test, respectively; Ɛ is the error between the slope of *ET*_{0} and the sum of contributions of each climatic factor to the change in *ET*_{0}. (The same meanings apply to Tables 3 and 4.)

Equations (5)–(9) were used to quantify the contribution associated with the four individual factors to the spatial variation in *ET*_{0}. The results implied that the downward trend in the mean annual *ET*_{0} 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*_{2} (5.3%) and *T*_{mean} (15.8%) on the decreased *ET*_{0} 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*_{2}, 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.

Transects | U_{2} | R_{s} | T_{mean} | e_{a} |
---|---|---|---|---|

Northwest to east | −0.35** | 0.55** | 0.45** | −0.52** |

Northwest to south | 0.37** | −0.03 | 1.70** | −2.27** |

Northwest to southwest | 0.47** | 0.39** | 0.28** | 0.22 |

Transects | U_{2} | R_{s} | T_{mean} | e_{a} |
---|---|---|---|---|

Northwest to east | −0.35** | 0.55** | 0.45** | −0.52** |

Northwest to south | 0.37** | −0.03 | 1.70** | −2.27** |

Northwest to southwest | 0.47** | 0.39** | 0.28** | 0.22 |

#### Along the northwest–south transect

The *ET*_{0} 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*_{0} decreased overall from northwest to south at the rate of −4.5 mm/10 km. Along this transect, *U*_{2} and *R*_{s} decreased significantly (*p* < 0.01), and their contributions to the decreased *ET*_{0} 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*_{0} 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*_{0}, such a trend might be offset by the function of *e*_{a} and *U*_{2}. The attribution results of the multiple regression also proved that *e*_{a} was the dominant factor that influenced the *ET*_{0} variation along the northwest–south transect.

It is worth noting that *ET*_{0} 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*_{2} made the maximum contribution to the decreased *ET*_{0} (58.8%), followed by increasing *T*_{mean} and *e*_{a} (23.2% and 15.7%), and the contribution of decreasing *R*_{s} was minimal.

#### Along the northwest–southwest transect

*ET*_{0} 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*_{0} in this direction. The maximum relative contribution to the decreased *ET*_{0} along this transect was *U*_{2}, 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*_{0} variation along this transect detected by the sensitivity method agrees with that by the multiple regression method.

ET_{0} | U_{2} | R_{s} | T_{mean} | e_{a} | Ɛ (%) | ||
---|---|---|---|---|---|---|---|

Ertuoke-Lingwu | Slope | −2.9** | −0.021** | 0.004* | 0.116** | 0.011** | |

Contributions (mm/10 km) | −5.71 | 0.65 | 4.85 | −2.72 | −0.01 | ||

Relative contributions (%) | 194.6 | −22.2 | −165.3 | 92.9 | −0.4 | ||

Lingwu-Zhongning | Slope | 10.6** | 0.043** | 0.020** | 0.049** | −0.001 | |

Contributions (mm/10 km) | 8.99 | 1.56 | −0.14 | 0.18 | −0.003 | ||

Relative contributions (%) | 84.9 | 14.8 | −1.3 | 1.7 | −0.03 | ||

Zhongning-Dongxiang | Slope | −9.4** | −0.032** | −0.043** | −0.106** | −0.001* | |

Contributions (mm/10 km) | −2.30 | −6.06 | −0.89 | −0.14 | −0.01 | ||

Relative contributions (%) | 24.5 | 64.6 | 9.5 | 1.5 | −0.05 | ||

Slope | −4.0** | −0.017** | −0.018** | −0.001 | 0.002** | ||

Ertuoke-Dongxiang | Contributions (mm/10 km) | −1.86 | −1.51 | −0.01 | −0.62 | 0.02 | |

Relative contributions (%) | 46.6 | 37.7 | 0.31 | 15.4 | 0.5 |

ET_{0} | U_{2} | R_{s} | T_{mean} | e_{a} | Ɛ (%) | ||
---|---|---|---|---|---|---|---|

Ertuoke-Lingwu | Slope | −2.9** | −0.021** | 0.004* | 0.116** | 0.011** | |

Contributions (mm/10 km) | −5.71 | 0.65 | 4.85 | −2.72 | −0.01 | ||

Relative contributions (%) | 194.6 | −22.2 | −165.3 | 92.9 | −0.4 | ||

Lingwu-Zhongning | Slope | 10.6** | 0.043** | 0.020** | 0.049** | −0.001 | |

Contributions (mm/10 km) | 8.99 | 1.56 | −0.14 | 0.18 | −0.003 | ||

Relative contributions (%) | 84.9 | 14.8 | −1.3 | 1.7 | −0.03 | ||

Zhongning-Dongxiang | Slope | −9.4** | −0.032** | −0.043** | −0.106** | −0.001* | |

Contributions (mm/10 km) | −2.30 | −6.06 | −0.89 | −0.14 | −0.01 | ||

Relative contributions (%) | 24.5 | 64.6 | 9.5 | 1.5 | −0.05 | ||

Slope | −4.0** | −0.017** | −0.018** | −0.001 | 0.002** | ||

Ertuoke-Dongxiang | Contributions (mm/10 km) | −1.86 | −1.51 | −0.01 | −0.62 | 0.02 | |

Relative contributions (%) | 46.6 | 37.7 | 0.31 | 15.4 | 0.5 |

Even though *ET*_{0} 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*_{0} decreased only at the rate of −2.9 mm/10 km; (2) from Lingwu to Zhongning County in Ningxia Province, *ET*_{0} 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*_{0} decreased at a rate of −9.4 mm/10 km to only 870 mm per year. *U*_{2} seems to be the dominant factor for the decreasing trend in *ET*_{0} in the first two regions, with 194.6% and 84.9% in the relative contribution; however, *U*_{2} decreased in the first region but increased in the second region. In the third region, the impacts of *U*_{2} on *ET*_{0} were degraded to only 24.5%; *ET*_{0} seems to be strongly associated with *R*_{s} (64.6%).

## DISCUSSION

### The possible impacts of topography and vegetation change on the spatial pattern of *ET*_{0}

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*_{0} 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*_{0}. 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*_{0}. Furthermore, there are several abrupt points of *ET*_{0} 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*_{0}. To further quantify the impacts of topography on the spatial change in *ET*_{0}, the correlation between the mean annual *ET*_{0} 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*_{0} (Figure 6(a)).

Moreover, the change in vegetation cover is another potential impact on *ET*_{0}. 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*_{0}. 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*_{0} was evaluated. The results showed that *ET*_{0} was significantly related to *NDVI* (*p* < 0.05), suggesting that the improvement of the vegetation condition would result in the *ET*_{0} 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*_{0}.

### The possible impacts of the spatial interpolation accuracy on the attribution results

There is no denying that the spatial interpolation accuracy of *ET*_{0} 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*_{0} 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*_{0} change was conducted and the results were compared with that of the actual *ET*_{0} for each station (Figure 7). The results showed that the contributions of the interpolated *U*_{2}, *R*_{s}, *T*_{mean}, and *e*_{a} variables to *ET*_{0} were similar to those of the four actual variables to *ET*_{0}, with a high determination coefficient R^{2} (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.

### Other uncertainties

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*_{0} 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*_{0}, errors exist between the calculated and observed total *ET*_{0} 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*_{0} change. Similar phenomena have been detected in the attribution analysis for the temporal variation of *ET*_{0} (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*_{0} 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*_{0}.

## CONCLUSIONS

In this study, the total differential equation was applied to analyze the contribution of four climatic factors to spatial variation in *ET*_{0} over the Loess Plateau. First, the annual *ET*_{0} 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*_{0} to *U*_{2}, *R*_{s}, *T*_{mean}, and *e*_{a} were calculated for each station. Then, the mean annual *ET*_{0}, four climatic variables and sensitivity coefficients were mapped using the cokriging interpolation method to show their spatial patterns. Three transects with distinguishable *ET*_{0} gradients were selected and sampled with an interval of 10 km for maps of *ET*_{0} and the four climatic variables. Finally, using distance as an independent variable, the contributions of each climatic variable to the spatial *ET*_{0} 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*_{0} 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*_{2}. The dominant factors that caused *ET*_{0} 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*_{0} 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*_{0}. This study provides a quantitative method to conduct the attribution analysis of the spatial *ET*_{0} 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.

## ACKNOWLEDGEMENTS

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