Analysis of variation in reference evapotranspiration and its driving factors in mainland China from 1960 to 2016

Understanding the variation in reference evapotranspiration (ETo) is vital for hydrological cycles, drought monitoring, and water resource management. With 1507 meteorological stations and 130 radiation-measured stations, the annual and seasonal ETo were calculated at each site from 1960 to 2016 in mainland China. The phenomenon of coefficient ‘a’ being less than 0.25 and coefficient ‘b’ being greater than 0.50 in the Angstrom–Prescott model occurred in almost the whole country, except for a small area of western and northeastern China. Moreover, the Xiao’s method was more applicable to calculate the net longwave radiation (R nl) and then improve the estimation accuracy of ETo. The annual ETo varied from 538.8 to 1559.8 mm and had a high-value center located in the plateau and desert of northwestern China and a low-value center located in Northeast China and near the Sichuan Basin. The spatial distribution of seasonal ETo was roughly similar to that of annual ETo, except for that in winter when ETo was high in the south and low in the north. In mainland China, the annual ETo decreased by 21.2 mm decade−1 because of the declining sunshine duration before 1993 and increased by 21.1 mm decade−1 due to the decreased relative humidity (RH) after 1993. Generally, the abrupt change of ETo mainly occurred in the southern China rather than northern China (except for Qinghai Tibet Plateau). Basically, the dominant driving factors of annual and seasonal ETo were RH and/or T max after the abrupt change in most parts of China.


Introduction
Actual evapotranspiration (ET a ), as the nexus of the water and energy balances, is a combined process whereby water is lost from the soil surface by evaporation and from vegetation by transpiration, thus linking the soil, vegetation and atmosphere together in terrestrial ecosystems (Allen et al 1998). Moreover, ET a is the crucial component of the hydrological cycle and climate system and is broadly applied in regional water resource assessments, drought monitoring and early warning, and evaluations of the water use efficiency of ecosystems (Hanasaki et al 2008, Anderson et al 2011, Yang et al 2015. Therefore, a suite of measurements or estimation methods have been developed to quantify ET a , including weighing lysimeters, eddy covariance, scintillometers, crop coefficient methods, and conservation of mass or energy methods (Norman et al 1995, Allen et al 1998, Bastiaanssen et al 1998, Yee et al 2015, Soubie et al 2016. However, observations of ET a are not only costly but also have insufficient time series, and ET a is strongly affected by surface characteristics (e.g. land cover and land change), climatic factors and the amount of available water, which constrain its applications in the aspect of climate change (Arnell andLiu 2001, McVicar et al 2012).
In contrast, the atmospheric evaporative demand (AED) represents the maximum evapotranspiration from a surface with an unlimited water supply  under certain meteorological conditions, and it is mainly characterized by three indicators: potential evapotranspiration (ET p ), reference evapotranspiration (ET o ) and pan evaporation (ET pan ) Thomas 2018, Xiang et al 2020). Among these metrics, ET o and ET pan are extensively investigated in terms of temporal variability and attribution analysis because the former is independent of crop type, crop development and management practices and the latter provides a measurement of actual AED; both these variables are comparable in disparate climatic regions worldwide (Allen et al 1998, Wang et al 2017, Fan andThomas 2018). However, the consistency of the pan evaporation data is not high enough in China for the reason that the type of evaporation pan was changed from D20 pan to E601B around 2000 but the homogenization has not been applied to the data , National Standard of People's Republic of China, GB/T35230-2017 2017). Many studies published at the beginning of the 21st century verified that most regions or countries, such as Australia, New Zealand, Alberta, the United States, China, and India, have suffered significant decreases in annual ET pan or ET o in the past 50 years (Chattopadhyay and Hulme 1997, Thomas 2000, Roderick and Farquhar 2002, Hobbins and Ramírez 2004, Xu et al 2006, Rayner 2007. Obviously, the aforementioned phenomenon conflicts with the projected capacity of the atmosphere to hold more water vapor as the air temperature increased by 0.74 • C during the 20th century (Roderick andFarquhar 2002, McVicar et al 2012). Further investigation of synchronous climate information indicated that decreases in wind speed (termed 'stilling') and/or radiation (termed 'dimming') offset warming, inducing a decline in the AED (Roderick and Farquhar 2002, McVicar et al 2012, Jahani et al 2017. This is a convincing interpretation because the variability in AED is affected by the interaction of multiple meteorological elements, especially wind speed, atmospheric humidity, radiation and air temperature.
In recent years, as the length of the time series of meteorological observations has been expanding, abrupt change in annual ET o have been detected in some regions. In Greece, the abrupt change in annual ET o occurred in the early 1980s; before the 1980s, ET o decreased over time, while after the 1980s, it increased (Papaioannou et al 2011). In China, trend analyses have been implemented in different administrative, geographic and climatic regions, and monotonic and nonmonotonic changes have been reported in previous studies (Huo et al 2013, Lv et al 2016, Gao et al 2017, Wang et al 2019. For example, Wang et al (2016) depicted that the annual ET o in the three-river source region decreased significantly by −9.1 mm decade −1 from 1980 to 2012; Wang et al (2019) showed that the annual ET o first decreased and then increased in mainland China (MLC) based on the China Meteorological Forcing Dataset from 1979 to 2015. In addition, the timing of abrupt change of annual ET o generally occurred in the 1980s or 1990s in China because existing studies were performed at disparate spatial scales (table 1). Synchronously, the dominant driving factors of annual ET o have been extracted by various approaches, such as multiple linear regression, differentiation equations, and Gaussian geographic weighted regression models (Wang et al 2016, Gao et al 2017, Zhang et al 2019. Reviewing the published literature, several shortcomings regarding the ET o trend and its contributing factors must be addressed in MLC. First, strengthening the analysis of the seasonal ET o pattern is necessary in light of the variability in meteorological elements in different seasons. Second, it should be noted that markedly disparate timing of abrupt change of annual ET o was detected by several studies with the same spatial extent in MLC (table 1). Hence, it is crucial to adopt higher-quality data to further explore whether differences in data richness and data sources influence the assessment results. Finally, it is not clear whether the main driving factors of the ET o sequence change before and after the abrupt change.
Consequently, the objectives of this study were to (a) analyze the spatial distribution of annual and seasonal ET o in MLC from 1960 to 2016; (b) investigate the annual and seasonal ET o series and determine the corresponding abrupt change points based on observations from higher-density meteorological stations; and (c) implement attribution analysis to explore the dominant driving factors of ET o before and after the abrupt change. The results are expected to serve as a reference for the evaluation of crop water demand, drought risk assessment and water resource allocation in MLC.

Study area
China, the third largest country in the world, is characterized by various and complicated climates (e.g. temperate monsoon climate, temperate continental climate, subtropical monsoon climate, tropical monsoon climate and plateau alpine climate) and geographies (e.g. plain, plateau, basin, hill and mountain). Based on the specific combinations of climatic resources and terrains, MLC can be divided into nine agricultural regions, including the Northeast China Plain (I), northern arid and semiarid region (II), Huang-Huai-Hai Plain (III), Loess Plateau (IV), Qinghai Tibet Plateau (V), Sichuan Basin and surrounding regions (VI), middle-lower Yangtze Plain (VII), Yunnan-Guizhou Plateau (VIII), and southern China (IX) (figure 1(a), www.resdc.cn/ data.aspx?DATAID=275). The multiyear averaged conditions of climatic elements are summarized in table S1 (available online at stacks.iop.org/ERL/16/ 054016/mmedia) for each region.

Data sources and processing
The daily meteorological observation data of 1507 high-density synoptic sites (figure 1(b)) were collected by the National Meteorological Information Centre of the China Meteorological Administration (http://data.cma.cn) from 1960 to 2016, consisting of average temperature (T), maximum temperature (T max ), minimum temperature (T min ), wind speed at a height of 10 m (U 10 ), relative humidity (RH), sunshine duration (SD), and precipitation (Pre). Moreover, 130 radiation-measured sites recorded the surface net radiation (R n ) and shortwave radiation received (R s ) and reflected (R ref ) by the surface (only a few sites measured all elements) and were evenly distributed in MLC ( figure 1(b)). However, the available number of radiation-measured sites varied over time because some of the stations had missing observations in certain years ( figure S1).
An entire year is separated into four seasons: spring (March-May), summer (June-August), autumn (September-November) and winter (December-February). To estimate the relevant parameters during different seasons, screening criteria were adopted; for example, if a meteorological element is completely missing in a season, the relevant data for that season is no longer used. And the percentage of missing data against the entire sequence from 1960 to 2016 was controlled within 3.5% for each climatic factor. Subsequently, two gap-filling methods were implemented for the time series of climatic factors: (a) linear interpolation was used when the number of continuous missing data points was less than 5, and these missing values were substituted with the mean value of other years when the number of continuous missing data points was greater than 4 for all climatic factors except for precipitation; and (b) replacement with the value on the same day from the nearest meteorological stations was used for only the precipitation data (Wu et al 2019).

Estimation of ET o
In this study, the FAO (Food and Agriculture Organization of the United Nations) Penman-Monteith equation was employed to calculate the daily ET o with basic meteorological parameters, which were defined as a hypothetical grass surface with a crop height of 0.12 m, a fixed surface resistance of 70 s m −1 , an albedo of 0.23, and an adequate water supply (Allen et al 1998). The formula is expressed as follows: where ET o is the daily reference evapotranspiration (mm d −1 ); ∆ is the slope of the saturated vapor pressure curve (kPa • C −1 ); γ is the psychrometric constant (kPa • C −1 ); e s is the saturation vapor pressure (kPa); e a is the actual vapor pressure (kPa); T is the mean air temperature ( • C) at a 2 m height, which is computed as the mean of T max ( • C) and T min ( • C); G is the soil heat flux (MJ m −2 d −1 ) and is so small that it can be ignored for a one-day period; and U 2 is the wind speed at a height of 2 m (m s −1 ), and a logarithmic wind speed profile relationship is capable of converting U 10 to U 2 (Allen et al 1998): where U 10 is the measured wind speed at a height of 10 m, and z is the height of the measurement above the ground surface (i.e. 10 m). R n is the net radiation at the canopy surface (MJ m −2 d −1 ), which is the difference between net solar shortwave radiation (R ns , MJ m −2 d −1 ) and net longwave radiation (R nl , MJ m −2 d −1 ) (Allen et al 1998): where α is the albedo or canopy reflection coefficient, which is 0.23 for hypothetical reference grass; R s is the solar shortwave radiation (MJ m −2 d −1 ); n is the actual sunshine duration (i.e. SD, h); R a , R so and N are the extraterrestrial radiation (MJ m −2 d −1 ), clear-sky solar radiation (MJ m −2 d −1 ) and the maximum possible duration of sunshine (h), which are calculated by equations (21), (36) and (34) in FAO paper 56, respectively; σ is Stefan-Boltzmann constant (4.903 × 10 −9 MJ K −4 m −2 d −1 ); T max and T min are the maximum and minimum absolute air temperatures (K = • C + 273.16); and a, b, c, d and f are empirical regression coefficients (dimensionless), of which default values are 0.25, 0.50, 0.34, −0.14, and −0.35, respectively, recommended by FAO paper 56. The 'a' and 'a + b' variables represent the fraction of extraterrestrial radiation reaching the earth on overcast days and on clear days, respectively. Coefficients a and b change depending on atmospheric conditions and solar declination (Allen et al 1998). The vast extent, diverse climate and complex terrain in China may have a significant effect on the values of coefficients a and b. Consequently, the transformation of the Angstrom-Prescott model (i.e. equation (5)) was used to estimate 'a' and 'b' by the least square method: where R s and n are the observations from the radiation-measured stations; and R a is extraterrestrial radiation; and N is maximum sunshine hours. Here, (R s /R a ) and (n/N) are the dependent and independent variables, respectively. Then, the coefficients 'a' and 'b' were estimated using the daily radiation data for each season from each available radiationmeasured station from 1960 to 2016. Besides, different methods (i.e. different combinations of coefficients c, d and f ) to estimate R nl have great effects on ET o (Yin et al 2008, Xiao andKong 2021). In addition to equation (6), another basic form to simulate R nl is as follow (Penman 1948, Doorenbos and Pruitt 1977, Xiao and Kong 2021: for the empirical coefficients in equation (8), a series of reference values (table 2) have been provided by Penman (1948), Doorenbos and Pruitt (1977) and Xiao and Kong (2021) based on observations. And we tried to compare these methods with R nl observation in order to acquire the best R nl model. The actual R nl is calculated as follow (Xiao and Kong 2021): where R s , R n and R ref (MJ m −2 d −1 ) are observations from radiation-measured stations, of which R ref is the reflected solar radiation by the surface. Subsequently, daily ET o could be computed based on climatic factors, the optimal R nl model and the corresponding coefficients 'a' and 'b' derived from the nearest radiation-measured stations. And mean values of ET o from all meteorological stations within a certain region were used to represent the regional conditions.

Trend detection
The nonparametric Mann-Kendall test was adopted to quantify the trend of ET o at a prescribed significance level of α because its application is not limited by the distribution of sample data (Hirsch and Slack 1984, Dinpashoh et al 2011,  Jhajharia et al 2012). Note that the pre-whitening process for the ET o time series was implemented first to avoid the impact of serial correlation of the sample on the Mann-Kendall trend test (Yue and Wang 2002). Additionally, the nonparametric Sen's slope estimator, developed by Sen, was used to evaluate the slope of the trend of a time series (Sen 1968).
In the study, the significance level α of the Mann-Kendall trend test was set to 0.05 or 0.01. For the standard normal test statistic, Z, a significant ascending trend exists in the time series if Z is greater than 1.96 or 2.58, and a significant descending trend exists in the time series if Z is less than −1.96 or −2.58. Influenced by natural or anthropogenic factors, the variability of hydrological and meteorological time series are often periodic or oscillating but seldom monotonous (She et al 2017). The moment that the trend of the time series changes from increasing (decreasing) to decreasing (increasing) is termed abrupt change (Villarini et al 2012, Gu et al 2017. The Mann-Kendall abrupt test was widely applied to detect the abrupt change points of time series in the field of hydrology and meteorology (Yang and Tian 2009, Liang et al 2010, Zhao et al 2014. However, the Mann-Kendall abrupt test may induce the phenomenon of multiple abrupt change points (Fang et al 2016, Wu et al 2019. Therefore, the non-parametric Pettitt approach was employed to implement abrupt change analysis for the reason that it can identify one change point at a time (Pettitt 1979, Zuo et al 2012, Lv et al 2016, Wang et al 2016.

Multiple standardized stepwise regression
ET o , which consists of radiative and aerodynamic terms, is comprehensively affected by diverse climatic factors, such as T, U 2 , SD, and RH, among which complex direct or indirect relationships exist. Moreover, the contribution of a specific climatic factor to ET o depends on not only the sensitivity of ET o but also its degree of change (Zhao et al 2020). To extract the dominant driving factors of annual or seasonal ET o , a multiple standardized stepwise regression model was employed to establish the relation between dependent variables (i.e. ET o ) and independent variables (i.e. T max , T min , U 2 , RH, and SD) (Wang et al 2017). Specifically, the analysis process included: (a) normalizing the climatic factors and ET o (with the mean and variance equal to 0 and 1, respectively); (b) setting the maximum p value and minimum p value to 0.05 for an independent variable to be added to and removed from the linear model, respectively; (c) implementing the model fitting with MATLAB R2014a and R square representing the model performance; and (d) calculating the ratio of the regression coefficient of each independent variable against the sum of all the regression coefficients as the relative contribution of the climatic factor to ET o ; the larger the contribution rate is, the more significant the climatic variable to ET o is, and vice versa (Wang et al 2017).

Evaluation metrics
The coefficients of variation (CV) is used to depict the temporal variability of the relevant elements. The variability levels were limited by CV ⩽ 0.1, 0.1 < CV < 1.0 and CV ⩾ 1.0, as weak, moderate and strong, respectively (Ayantobo et al 2017).
The performance of the R nl models reported in table 2 were evaluated by the correlation coefficient (R), mean bias error (MBE) and root mean square error (RMSE) of the simulated values to measured values. The smaller the RMSE, the more accuracy of the R nl model. MBE represent the estimation error, positive and negative means higher and lower estimation, respectively (Itenfisu et al 2003, Yin et al 2008. The R, MBE and RMSE were given by: where S i and O i are simulated and observed values of the ith sample, respectively;S andŌ are the mean values of simulated and observed samples, respectively; and n is the sample size.

Spatiotemporal variability of coefficients 'a'
and 'b' in the Angstrom-Prescott model 3.1.1. Spatial distribution of coefficients 'a' and 'b' The multiyear average of coefficients 'a' and 'b' in the Angstrom-Prescott model, based on the seasonal fitted values, was calculated for 87 radiation-measured stations from 2000 to 2015 (relatively more stations were available within the period). Then, the ordinary kriging method was adopted to interpolate the 'a' and 'b' coefficients and analyze the spatial heterogeneity in MLC (figure 2). In most parts of MLC, coefficient 'a' was lower than 0.25, except in the northern zone I, northeastern zone II and western zone V ( figure 2(a)). The stripes extending from Northeast China to Western China and extending from Southeast China to Northwest China were roughly viewed as the high-value and lowvalue regions of coefficient 'a' , respectively. Furthermore, the whole middle-lower Yangtze Plain presented almost the lowest value of coefficient 'a' , with a value less than 0.19 ( figure 2(a)).
In most parts of MLC, coefficient 'b' was greater than 0.50, except in northeastern China and western zone II ( figure 2(b)). The subregion, constituted by zones I, II, III, and IV, had a coefficient 'b' below 0.54. However, coefficient 'b' in the rest of MLC was greater than 0.54. In addition, several high-value centers appeared in the northwestern zone II, southwestern zone V, western zone VI, and southeastern zone VII, with coefficient 'b' values greater than 0.56 ( figure 2(b)).

Temporal variation in coefficients 'a' and 'b'
Limited by the continuity of radiative data records, the temporal variation in coefficients 'a' and 'b' was analyzed from only 18 radiation-measured stations evenly distributed throughout MLC (figure S2). The annual and seasonal patterns of coefficient 'a' showed an approximately slight downward trend from 1960 to 1990, followed by an upward trend from 1990 to 2016 (figure 3). The multiyear average of seasonal coefficient 'a' gradually increased from spring (and summer) to winter (0.16, 0.16, 0.18, and 0.20 in sequence), and low variability with a coefficient of coefficient (CV) less than 0.10 was present in each season (table 3). The annual and seasonal coefficient 'b' fluctuated and decreased from 1960 to approximately 2005 but increased slightly in recent years (figure 3). Coefficient 'b' was largest in spring, with a multiyear average of 0.56, followed by that in autumn, summer and winter. A smaller variability (CV below 0.07) existed in the four seasons compared to that of coefficient 'a' . Moreover, a significantly moderate negative correlation was detected between coefficients 'a' and 'b' at different time scales (table 3). One of the main reasons for this result is that the sum of coefficients 'a' and 'b' is constrained to less than 1.

Comparison of estimation methods of R nl
Seventeen radiation-measured stations in MLC were selected to analyse the applicability of the four estimation methods of R nl because of R s , R n and R ref were observed in these sites. Table 4 depicted the site locations and sample size of the observations. The actual R nl could be obtained by equation (9) and the most suitable R nl model could be extracted by means of analysis of relevant statistical indicators. The correlations between four R nl methods and actual observations were basically similar, of which Xiao's method had a slightly stronger relationship with observations (i.e. mean value of R was 0.73 within 17 sites) ( figure 4(a)). Moreover, the accuracy of R nl simulation by FAO56 and FAO24 was extremely low and the RMSE were 3.77 and 3.49 MJ m −2 d −1 , respectively. However, the estimation methods of R nl developed by Penman and Xiao could calculate the net longwave radiation well and Xiao's method had the highest accuracy with RMSE equal to 2.14 MJ m −2 d −1 (figure 4(b)). In addition, estimation methods of R nl by FAO56 and FAO24 completely underestimated the net longwave radiation on the 17 stations and the mean values of MBE were −3.12 and −2.77 MJ m −2 d −1 , respectively. Nevertheless, the mean value of MEB between simulated R nl by Xiao's method and observations was approximately close to 0 (MBE was equal to −0.22 MJ m −2 d −1 ) (figure 4(c)). In conclusion, Xiao's method had a high accuracy to estimate R nl and was employed to compute ET o as a crucial sub procedure in the study.

Spatiotemporal variability of ET o 3.3.1. Annual and seasonal patterns of the spatial distribution of ET o
The multiyear average annual ET o ranged from 538.8 to 1559.8 mm in MLC. Generally, the inland ET o in the northwest was higher than that in the central and eastern parts of China. The high-value region with an ET o above 1100 mm mainly presented in the central and western parts of zone II and the western part of zone V, where the plateau and Gobi Desert are located. The low-value regions with an ET o below 700 mm were located in the northeastern part of the MLC and adjacent areas among zone VI, zone VII and zone VIII (figure 5).  The spatial distribution of ET o in spring, summer and autumn was basically similar to its annual pattern, except for the magnitude of ET o , the slightly shifted location and the diverse extent of high-value and low-value regions. For example, the locations of the relatively low-ET o region varied from South China (surrounded by a 200 mm contour) to Southwest China (surrounded by a 300 mm contour) to Central China (surrounded by a 150 mm contour). However, an obvious graded distribution of ET o was present in winter, gradually increasing from Northwest and Northeast China to Southwest China. Overall, the ET o was largest in summer, followed by that in spring, autumn, and winter (figure 5).

Temporal variation in ET o and its abrupt change
An abrupt change point that occurred in 1993 was detected in the time sequence of the annual ET o from 1960 to 2016 in MLC; before this point, the annual ET o decreased by 21.2 mm decade −1 , while after this point, the values displayed an increasing trend of 21.1 mm decade −1 (figures 6 and 8). Moreover, the ET o sequences during spring and summer were also separated into an increasing and a decreasing segments due to the abrupt change, and the timing of these abrupt change points were 2004 and 1996, respectively (figures 6 and 7).Compared with the ET o time series in summer, the ET o sequence in spring decreased more slowly before abrupt change and increased more faster after abrupt change (figure 8). However, there were no abrupt change points for the ET o sequences during autumn and winter and the insignificant increased trends were detected with 0.6 and 1.0 mm decade −1 , respectively (figure 8).
Generally, the abrupt change had not been detected for most annual and seasonal ET o time series in northern China, such as zone I, zone II, zone III and zone IV (except for zone V) (figure 7). However, a majority of the annual and seasonal ET o time series in the Qinghai Tibet Plateau (i.e. zone 1 9 9 3 1 9 9 2 2 0 0 4 1 9 7 4 1 9 9 7 1 9 9 2 2 0 0 3 1 9 9 3 1 9 9 3 1 9 9 6 1 9 9 0 1 9 7 5 1 9 9 2 1 9 7 4 1 9 9 8 1 9 9 5 1 9 9 2 1 9 9 6 1 9 9 5 1 9 9 3 1 9 9 7 1 9 9 7 1 9 9 7 1 9 9 4 MLC I II III IV V VI VII VIII IX Annual Spring Summer Autumn V) and the southern China were disturbed by the abrupt change, which caused a descending sequence and a subsequent ascending sequence (figures 7 and 8). The timing of the abrupt change of the annual ET o sequence in subregions was particularly early or late, such as that in zone VII, which had an abrupt change point appearing in 1974, and that in zone VI, which had an abrupt change point appearing in 2004. However, the abrupt change of the seasonal ET o time series generally occurred in approximately 1995, especially in autumn and winter ( figure 7).
The annual ET o in zone I and II and zone III and IV continuously but insignificantly increased or declined. Moreover, the annual ET o in zone V and IX decreased by more than 30 mm decade −1 before abrupt change, which was higher than that in MLC. And the annual ET o in zone VI and VIII increased by more than 40 mm decade −1 after abrupt change, which was far larger than in MLC or other subregions. Basically, the decrease in ET o in summer and the increase in ET o in summer and autumn played important roles in the corresponding trend of the annual ET o before and after the abrupt change in most subregions (figure 8).

Driving factors of annual and seasonal ET o
The annual ET o in MLC was regulated by RH, SD and T max , which contributed 48.5%, 27.1%, and 24.4%, respectively, to the variability in ET o from 1960 to 2016. However, the driving factors of annual ET o were completely disparate before and after the abrupt change, which were SD and RH respectively. For the whole seasonal ET o , the RH was the dominant driving factor in spring, summer, autumn, and winter, and the contributions were 39.7%, 100%, 34.8%, and 58.7%, respectively. Furthermore, the driving factors of seasonal ET o before and after abrupt change during spring and summer were completely disparate. Before abrupt change, the seasonal ET o were mainly affected by SD and RH in spring and summer, respectively. However, the dominant driving factors were U 2 and T max after abrupt change in spring and summer, respectively (figures 9 and 10). Affected by diverse climates and complex terrains, the drivers of annual and seasonal ET o were significantly distinct in the different subregions. For the annual ET o in each subregion, the RH in southern China (zone VI, VII, VIII, and IX) and zone II and IV, the T max in zone I and V, and the U 2 in zone III were the corresponding dominant driving factors of ET o . In spring, the seasonal ET o in central and eastern China (i.e. zone III, IV, VI, and VII) was mainly affected by RH. In summer, RH had an important D Wu et al 1960D Wu et al 1970D Wu et al 1980D Wu et al 1990D Wu et al 2000  impact on the variability of seasonal ET o in northern China (i.e. zone I, II, III, IV, and V). In autumn, the dominant driving factor was RH in northern China, while the T max was the main driver in southern China. However, the main driving factors were approximately opposite in winter compared to that in autumn (figure 11). Given that the abrupt change mainly occurred in southern China, the analysis of trend change was implemented in zone VI, VII, VIII, and IX. Both before and after abrupt change, the RH and SD were the dominant driving factors of annual ET o in zone VII and VIII and zone VII and IX, respectively. After abrupt change, the seasonal ET o was mainly influenced by SD in spring and winter in southern China. However, the T max was the dominant driving factor of seasonal ET o in summer and autumn after abrupt change. It should be noted that the SD was always the dominant driver of annual and seasonal ET o in zone IX after abrupt change (figure 11).

Discussion
The ET o is determined by two crucial components: radiative and aerodynamic terms that are comprehensively influenced by climatic factors, such as T, SD, U 2 , and RH (Allen et al 1998, McVicar et al 2007, Dinpashoh et al 2011, Jhajharia et al 2012. In MLC, the high-value regions of annual ET o were concentrated in northwestern and southwestern, where the annual ET o exceeded 1100 mm. However, northeastern China, the Sichuan Basin, and the middle and lower reaches of the Yangtze River were low-value regions of annual ET o , and some of these areas had an annual ET o less than 700 mm. Referring to table S1, zones II and IX and zones I and VI showed relatively  high and low annual ET o values, respectively. Basically, the district characterized by high air temperature, low atmospheric humidity, strong wind, and long sunshine duration (enhancing the holding capacity of water vapor in the atmosphere, improving the ability of water vapor transport, and meeting the energy requirements of the evaporation surface) had a higher ET o , and vice versa. Trend analysis of the annual ET o calculated with high-density meteorological stations in MLC again verified that the abrupt change in the ET o sequence induced a downward trend and a subsequent upward trend (Fan et  in Southwest China, which was basically consistent with the result (1997) of abrupt change for zone VIII in the study. In addition, there had few abrupt change of the ET o sequence by the Pettitt test in northern China. However, the abrupt change of ET o has been detected in some basins or provinces in northern China (Zuo et al 2012, Lv et al 2016. The difference of abrupt change point among the relevant studies is principally induced by four aspects that are the calculation procedure of ET o , the study area, data source, and the detection method of abrupt change. Compared to the original FAO 56 Penman-Monteith equation, using the corrected R ns and R nl procedure can obtain more accurate ET o (Yin et al 2008, Xiao andKong 2021). In addition, the difference in the aspects of the study area and data source (grids or sites, the quantities of meteorological stations) can also cause the disparity of ET o . Furthermore, a series of methods can be adopted to detect the abrupt change of ET o sequence, such as Mann-Kendall test ( The driving factors of ET o are determined by two aspects: the sensitivity of ET o to parameters and the degree of variability of these parameters (Zhao et al 2020). Local sensitivity analysis methods, such as partial derivatives, have been widely employed to detect the sensitivity of ET o to several climatic elements because of the simplicity of the algorithm in China (Huo et al 2013, Fan et al 2016, She et al 2017, Jiang et al 2019, Wang et al 2019, Zhao et al 2020. Nevertheless, a shortcoming of the local approach is that the effect of the variation in a factor on the dependent variable does not include its interactions with other factors. Therefore, a global method named the extended Fourier amplitude sensitivity test (eFAST) (Saltelli et al 1999(Saltelli et al , 2005) was adopted to analyze the first-order sensitivity (single effect) and total sensitivity (mutual effect) of 10 parameters to ET o . The most sensitive factors were the T max and RH, followed by the U 2 and T min . In contrast, the sensitivities of SD and other parameters were quite small (figure S3). Fan et al (2016), Wang et al (2019) and Zhao et al (2020) also indicated that the RH and T max presented similar high sensitivities.
Although the sensitivities of SD to ET o were relatively low, the very significant downward trend (−68.63 h decade −1 ) of SD indicated that the annual ET o decreased by 21.2 mm decade −1 from 1960 to 1993 ( figure S4). This so-called phenomenon of 'dimming' was consistent with the analyses by Roderick and Farquhar (2002) and Jahani et al (2017). After 1993, however, the decreasing trends of SD were relatively weak, and simultaneously, the annual ET o increased by 21.1 mm decade −1 as the decreasing RH which induced an increase in the vapor pressure deficit of the atmosphere ( figure S4). Moreover, the RH and/or T max , as dominant factors, always had a crucial effect on the annual and seasonal pattern of ET o in most parts of China. The increase in ET o may intensify the hydrological cycle in humid southern China and exacerbate the drought risk and pressure on water resources in the arid and semiarid regions of northern China. The result about trend change and driving factors of ET o are expected to serve as a reference for the evaluation of crop water demand, drought risk assessment and water resource allocation in MLC.

Conclusions
Coefficients 'a' and 'b' in the Angstrom-Prescott model were lower and higher than the reference values (i.e. 0.25 and 0.50) of FAO paper 56, respectively, in most parts of China, except for a small area of western and northeastern China. The temporal variations in the two parameters were all characterized by annual and seasonal patterns that first decreased and then increased from 1960 to 2016, but they had different change points. In addition, the Xiao's method to estimate R nl can greatly improve the accuracy of ET o .
In MLC, the multiyear average annual ET o was from 538.8 to 1559.8 mm, with a high-value region exceeding 1100 mm located in the plateau and desert areas of western and northwestern China and a lowvalue region below 700 mm located in northeastern China and near the Sichuan Basin. Moreover, the seasonal ET o was highest in summer, followed by that in spring, autumn, and winter. The spatial distribution of seasonal ET o was roughly similar to that of annual ET o , except for that in winter when ET o was low in northern China and high in southern China.
An abrupt change point that appeared in 1993 separated the annual ET o sequence into two segments: a decreasing interval with a slope of 21.2 mm decade −1 and an increasing interval with a slope of 21.1 mm decade −1 in MLC. Except for Qinghai Tibet Plateau, few abrupt change points of the ET o sequences have been detected in northern China. Although the abrupt change of annual ET o sequences occurred in most of the southern China, the timing of the abrupt change was not consistent. However, the abrupt change of the seasonal ET o in southern China generally occurred in approximately 1995, especially in autumn and winter. In MLC, the driving factors of annual ET o was SD before 1993 and RH after 1993. Basically, the variation in the annual and seasonal ET o was mostly controlled by RH and/or T max in most parts of China.

Data availability statement
The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.