Variations in Reference Evapotranspiration and Associated Driving Forces in the Pearl River Delta of China during 1960 – 2016

Recent climate warming and rapid urban development in the Pearl River Delta (PRD) of China exerted great impacts on the reference evapotranspiration (RET), which in turn affects the management of water resources and the quality of the urban environment. The objectives of this study are to examine (i) the temporal variability of the RET in the PRD and (ii) the underlying causes responsible for the temporal variation in the RET across space inside the PRD. The results indicate the following. (1) The RET in the PRD had an overall increasing trend caused by the increase of construction land during 1960 – 2016. (2) The increase of surface albedo caused by land cover conversion from woodland to grassland played an important role in the noticeable decline of the RET in Guangzhou and Zengcheng. (3) The dominant factors triggering RET variation varied across space in the PRD. In detail, the decline of sunshine duration (SD) decreasing Rn and the decline of wind speed (WS) weakening energy exchange were the dominant factors in decreasing RET in Guangzhou and Zengcheng. In contrast, the daily maximum temperature, daily minimum temperature, and relative humidity (RH), which were the factors causing the increase of vapor pressure deficit (VPD), were responsible for the RET increase in Taishan, Zhongshan, and Shenzhen. Overall, our results indicated that the RET in PRD exhibited a strong spatial heterogeneity due to differences in land use change and climatic conditions. Therefore, the improvement of water resources management and urban environment in the PRD should consider the spatial variation and underlying forces of RET changes.


Introduction
The Pearl River Delta (PRD) is one of three major regions in China that gather a large population size and has a strong innovation ability and rapid economic growth. It plays a very important role in the social and economic development of China and even the world. Economy and urbanization in the PRD have been growing rapidly since China's implementation of reform and opening up in the 1980s. The consequent effects of urbanization on the climate and hydrological processes have become more and more remarkable . For example, evapotranspiration (ET) can change because of the alteration in land surface condition under rapid urbanization. This will exacerbate the pressure on the available water resources because of the augmentation of water loss resulting from higher ET, eventually impeding urban development in the PRD where industry, agriculture, and population are highly concentrated (Thomas 2015;Zhang et al. 2013). Therefore, it is necessary to perform a thorough analysis on ET changes under the background of rapid urbanization and economic growth. Such study can provide scientific bases for water resources management in the PRD.
In recent decades, anthropogenic climate warming has become an indisputable fact, whose effects on ET differ across regions and even within the same region. As a result of climate warming, the average ET in China experienced a downward trend over the past few decades (Gao et al. 2006;Li et al. 2017;She et al. 2017;Thomas 2015;Wang et al. 2017), especially in southwestern, northwestern, southeastern, northern, south-central, and eastern China (Feng et al. 2012Huo et al. 2013;Li, X. et al. 2014). For example, a list of studies indicated that ET in major basins of China, including the Laohahe River basin, the Haihe River basin, the Jinghe River basin, the Poyang Lake basin, the Pearl River basin, the Liaohe River basin, the Yangtze River basin, and the Huaihe River basin, trended downward obviously Liu et al. 2012;Liu et al. 2015;Xu et al. 2015;Ye et al. 2015), whereas it trended upward in the Upper Heihe River basin, Songhua River basin, and Yellow River basin (Liu et al. 2012;Luo et al. 2017). To a large extent, it is widely recognized that the "evaporation paradox" was prevalent in most parts of China under climate warming. However, existing studies on ET variations were mostly concentrated in northern China and a few were carried out in southern China. This inadequacy might prevent us from under standing the spatial heterogeneity of ET variations because human activities in southern China, such as in the PRD, are more intense and thus more profoundly impact the natural environment, which in turn affects ET in these regions.
Although the ET variation in most regions of China showed high consistency, its sensitivity to changes in natural environmental factors varies across regions. For example, ET is very sensitive to wind speed (WS), relative humidity (RH), and temperature in northwestern, southwestern, and northern China (Feng et al. 2012Huo et al. 2013;Liu et al. 2012). The spatial heterogeneity of ET's sensitivities to dominant environmental factors across regions implies that ET may still vary inside a region because of the spatial variation in environmental conditions. However, because the dominant factor in triggering the ET variation in a small area may not change significantly or differ only slightly across space, it is likely that shifts in such factor may not lead ET to change significantly. By contrast, the environmental factors that experienced significant changes are likely to be responsible for triggering ET variations. For example, WS was not only a sensitive factor in northwestern China, but also the dominant factor triggering the decrease of ET (Huo et al. 2013;Liu and Zhang 2013;Thomas 2015;Zheng and Wang 2015). ET is very sensitive to RH in southwestern China, whereas the decrease in WS is the dominant factor actual causing the ET to decrease in this region (Li, Z. et al. 2014). Similarly, ET is sensitive to temperature variations in northern China, while the decrease in RH is responsible for the decrease of ET in this region (Ye et al. 2015). Differences in sensitive and actual controlling factors in inducing ET changes suggest that they may still be applicable in the same region owing to the complexity of interaction among environmental factors.
The rapid urbanization in China has been found to affect the actual ET at both local and regional scales. For example, the overall downward trend of ET in China in recent decades was partially a result of the Citation Liu, Y., G. Tang, L. Wu, Y. Wu, and M. Yang, 2019: Variations in reference evapotranspiration and associated driving forces in the Pearl River Delta of China during 1960-2016. J. Meteor. Soc. Japan, 97, 467-479, doi:10.2151/jmsj.2019 continuous expansion of urbanization (Liu et al. 2008). Although land cover change might have little impacts on the variation of regional ET in Guangdong, China, where precipitation is greater than potential ET, it can still greatly influence ET when the precipitation is less than ET (Chen et al. 2015). It has been indicated in an earlier study ) that the reduction of forests can bring about greater effects on ET than changes in other land cover types. Besides, the topography, saturated vapor pressure, rainfall days, soil water content, and WS can also affect ET (Li, X. et al. 2014;Li, Z. et al. 2014;Shang et al. 2017;Tong et al. 2017;Wu et al. 2017). Although it has been indicated in existing studies that either climatic or human factors had the most pronounced impact on ET, the effects of their interactions on ET change are often overlooked because the two factors were treated independently in most of the earlier studies.
The objectives of this study are to (i) analyze the variations of reference evapotranspiration (RET) in the PRD of China and (ii) examine the dominant factors in controlling RET variations in the PRD. To achieve these goals, the FAO Penman-Monteith equation was used to estimate the time series of daily RET during 1960 -2016. Based on derived RET data, statistical methods such as linear regression, the Mann-Kendall trend test, and spatial statistics were used to analyze the spatiotemporal variations of RET in the PRD. The Canonical Correlation Analysis (CCA) was used to determine the dominant factor responsible for the RET change. The results of this study can provide us with a better understanding of the historical RET variation and underlying mechanisms in the PRD and can eventually provide a scientific basis for better managing water resources and urban development in the PRD.

Data source
The 500 m and eight-day composite MOD09A1 and MOD15A2 products for the period of 2000 -2016 were used in this study. These data were obtained from the Land Processes Distributed Active Archive Center (LP DAAC) (https://lpdaac.usgs.gov/). In specific, the MOD09A1 products were used to calculate the albedo by referring to Liang (2000), and the MOD15A2 products were used to calculate the leaf area index (LAI).
Climate data for the period of 1960 -2016 were obtained from the Chinese Daily Terrestrial Climate Dataset 3.0, provided by the China Meteorological Data Sharing Service Network (http://data.cma.cn/). The climatic variables used in this study include the time series of daily total precipitation, mean tempera-ture (T mean ), maximum temperature (T max ), minimum temperature (T min ), sunshine duration (SD), RH, and WS. We deleted the climate stations which data missing of any variable are more than 10 days in a row. As a result, climatic variables from five national benchmark climate stations ( Fig. 1) were selected and used in this study. For the five stations, we fixed missing values using the Inverse Distance Weighted approach. The resultant climate data were used to estimate the RET by the FAO Penman-Monteith equation and analyze the effects of variation in climatic factors on the RET in PRD. Although there are only five national benchmark climate stations, their effective coverage includes most of the PRD because each benchmark station can represent the climatic condition in an area of circle with a radius of 45 km according to the current Chinese regulations (China Meteorological Administration 2007) (Fig. 1). Thus, these five stations can to a large extent represent the climatic condition of the whole PRD.
The land use/cover data for the period 1980 -2010 were provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn). These data are at a 30 m spatial resolution and annual scale. The types of land use/cover in a year were mainly based on the satellite data of that year but were adjusted by referring to the satellite data before and after that year (e.g., the types of land use/cover in 2010 were obtained from 673 graphs of Landsat TM for the period 2009 -2011) (Xu et al. 2014). We classified the land use/cover in the PRD into six types: cropland, woodland, grassland, water, construction land, and unused land. The bound- ary of the PRD (Fig. 1) was provided by the National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China (http://lake.geodata.cn/).

Methods a. Estimate of RET
The RET was estimated by the FAO Penman-Monteith equation. The so-called reference crop ET was based on a list of assumptions, including the reference surface being a grass reference crop of 0.12 m height, a fixed surface resistance of 70 m s −1 , and an albedo of 0.23. In reality, the reference surface resembles an extensive surface of green, well-watered, and actively growing grass of uniform height. In addition, the uniform grass was assumed to completely shade the ground (Allen et al. 1998). The FAO Penman-Monteith equation is expressed as follows: The meaning of each parameter is listed in Table 1. R a depends on the latitude and the number of the day. R SO is relatively stable according to Eq. (6). Among them, T mean for the 24-hour period is defined as the arithmetic mean of T max and T min . The detailed calcu-lation of other parameters is introduced by Allen et al. (1998).

b. CCA
The CCA and canonical variables were put forward originally by Hotelling (1936). This approach can be used to analyze the correlations between two sets of variables. In this study, the CCA was used to explore the relationship between the RET and a climatic factor. The data used in the CCA were conducted on an annual basis. According to Liu et al. (2016), the specific processes for conducting the CCA are as follows: (1) First, the RET is assigned as a dependent variable Y, and the climatic factors (SD, T mean , T max , T min , WS, and RH) are treated as independent variables X.
(2) Second, calculate the canonical variables and conduct significance tests and redundancy analyses.  The combination of the sign of C L and the slope of a climatic factor help determine whether a climatic factor acts to enlarge or reduce the trend of RET. For a given climatic factor, when C L is positive (or negative), it is positively (or negatively) correlated with RET. Correspondingly, any increase/decrease in the climatic factor causes the RET to increase/decrease (or decrease/increase).
In addition, the absolute values of C L and the significance level of the slope can suggest the strength of the effects on RET. C L of a climatic factor can be regarded as its influencing weight on RET and the coefficient of multivariate regression model. The larger the value of C L is, the greater the influence of an independent variable on the fitting value of a multivariate regression model is. Further, because C L is dimensionless, it can be used not only for withingroup comparisons, but also for between-group comparisons.

Evaluation of estimates of RET
We compared RET estimates on the basis of the Penman-Monteith equation with those based on either small evaporators (ET L ) or large evaporators (ET B ) in the PRD. This is because the ET L data spanned from 1960 to 2001, after which the ET L approach was no longer used and shifted to the ET B approach. However, because data missing of ET L and ET B were severe for the Zhongshan station, this station was excluded from the RET estimates. As a result, field measurements of ET L and ET B data from four weather stations (i.e., Guangzhou, Zengcheng, Taishan, and Shenzhen) were used in evaluating the RET estimates. Figure 2 shows the relationship between RET and ET L for the period 1960 -2001 and that between RET and ET B for the period 2002 -2012 in Guangzhou, Zengcheng, Taishan, and Shenzhen. According to the scatterplots and calculated R-squared values (Fig. 2), although RET is more strongly positively correlated with ET L (0.76 ≤ R 2 ≤ 0.79) than with ET B in all four  Figure 3 shows the trends and variations of the estimated RET in each of the five stations and that averaged for the whole PRD. According to Fig. 3, the RET in Guangzhou and Zengcheng decreased by a rate of 0.85 mm yr −1 , whereas it increased by a rate of 1.23 mm yr −1 in Taishan, 2.47 mm yr −1 in Zhongshan, and 1.15 mm yr −1 in Shenzhen. As the increasing rates of RET in Taishan, Zhongshan, and Shenzhen were greater than the decreasing rates of RET in Guangzhou and Zengcheng, the average RET for the whole PRD showed an overall upward trend by a rate of 0.63 mm yr −1 during 1960 -2016 (Fig. 3). Table 2 and Fig. 4 show the variations of land use/ cover change (LUCC) in the PRD during 1980 -2010.

Analyzing the driving forces for RET a. Land use and land cover change
According to the temporal variation of LUCC, the PRD was dominated by woodland and cropland during 1980 -2010, accounting for 43.0 % and 33.7 % of the total area, respectively. However, the area of cropland had been greatly reduced by the end of 2010 and was 10.6 % lower than that in 1980. Also, the area of woodland decreased by 2.3 % during 1980 -2010. The area of grassland and unused land was very small, accounting for 1.5 % and 0.1 % of the PRD, respectively. In contrast, the coverage of waters accounted for 10.3 % of the PRD and increased only slightly. Construction land accounted for 11.4 % of the PRD but increased consistently during 1980 -2010. For example, the construction land increased rapidly up to 19.7 % of the PRD in 2010, which was almost triple that (6.9 %) in 1980.
There existed three major conversion types of LUCC in the PRD: conversion of cropland and woodland to construction land, conversion of cropland to waters, and conversion of waters to cropland ( Table 2). The increase in the construction land was mainly from the cropland, followed by woodland, suggesting that urbanization in the PRD has been increasing rapidly so that a large number of croplands and woodlands were used for urban construction, especially in Guangzhou and Shenzhen. For example, a large number of croplands were distributed around Guangzhou in 1980 (Fig. 4a), but most of them were converted into construction lands in 2010 (Fig. 4d). Similarly, a large number of woodlands were distributed around  Shenzhen in 1980 (Fig. 4a), but most of them were converted into construction lands in 2010 (Fig. 4d).
The conversion of a large number of croplands to construction lands caused RET to increase in the PRD during 1980 -2010. The more the construction land increased, the more the RET increased (Fig. 5). For example, the construction land increased slowly in the 1990s, during which the average RET was only 25 mm higher than that in the 1980s. In the 2000s, when the construction land rapidly increased, the RET was 44 mm higher than that in the 1990s. This suggests that RET in the PRD was obviously affected by land use change that affects land surface albedo (see Section 4.1). Figure 6 shows trends and variations in T mean , T max , and T min in each station and the whole PRD. The calculated statistics suggest that these climatic indicators had significant upward trends ( p < 0.05) during 1960 -2016 in Zengcheng, Taishan, Zhongshan, and Shenzhen. As a result, when averaged for the PRD, T mean also increased rapidly in the PRD. However, the increasing rate of T min was greater than that of T max in both Taishan and Shenzhen, indicating that the diurnal temperature difference was reduced.

c. Effects of climatic factors on RET
The C C value calculated for each station between the RET and each of the climatic factors was approximately 1.0 and was significant at the 0.01 significance level (Table 3), indicating that the RET was significantly correlated with the variation of climatic factors, although the strength of correlation differs among climatic variables.
(1) Guangzhou and Zengcheng The SD was the dominant factor in triggering the RET variation. The C L values of T max and T min were smaller than those of other climatic factors, indicating that the temperature had a little influence on the RET. In addition, the RET was positively correlated with SD, T max , and WS, whereas the correlation between T min and RH was negative. The decreasing RH played the main role in promoting RET with a combination of effects from other climatic factors. In contrast, the decreasing SD and WS played the main role in decreasing the RET. The forces that accelerate RET to decrease were obviously stronger than that causing the RET to increase. As a whole, the RET showed an overall decreasing trend.
(2) Taishan and Zhongshan T max was the dominant factor affecting the RET variation, whereas WS had the least influence on the RET. The RET was positively correlated with each of SD, T max , and T min . In contrast, RH was negatively correlated with the RET. The increase in T max and T min and the decrease in RH jointly played a dominant role in enlarging the RET. Although WS had an obviously downward trend, the influencing weight of WS was so small that it had the least impact on the RET variation. In addition to WS, all other climatic factors experienced obvious changes, which to some degree are contributive to the RET increase. Thus, the magnitude of RET increase in Taishan and Zhongshan was greater than in the other stations.
(3) Shenzhen T max was the dominant factor controlling the RET variation, and SD had the least influence on the RET, although the RET was positively correlated with SD, T max , T min , and WS. In contrast, it was negatively cor- related with RH. The increase in T max and T min and the decrease in RH were the main factors in speeding up the RET increase, whereas the decrease in WS was the main factor in decreasing the RET. The influencing weights of T max , T min , and RH that played the major roles in speeding up the RET decreases were high, whereas the influencing weights of WS that played the major inhibiting roles were low. Therefore, the role of climatic factors in promoting RET was obviously remarkable.

Effects of LUCC on RET
It was suggested in Fig. 5 that the increase in the RET was affected by the increase in the construction land. By calculating the average albedo of each land use type in the PRD, we found that the surface albedo (in an descending order) was about 27 % for unused land, 21 % for grassland, 19 % for cropland, 17 % for construction land, 14 % for waters, and 13 % for woodland. A comparison of land use types between 2000 and 2010 showed that construction land in areas surrounding Taishan, Zhongshan, and Shenzhen stations remained unchanged, whereas woodlands in areas surrounding Guangzhou and Zengcheng stations were converted to grasslands (Fig. 4). Because the albedo for woodland was comparatively the lowest among the five land use types, the conversion of woodland to grassland increased the surface albedo in Guangzhou and Zengcheng. In addition, the conversion of woodland to grassland in Guangzhou and Zengcheng decreased LAI significantly ( p < 0.01) by a rate of 0.03 and 0.02 units per year, respectively (Fig. 7a). In fact, our results suggest that LAI was negatively correlated with albedo in most areas of the PRD. The smaller the LAI was, the greater the albedo was (Fig. 7b). This explains that the decrease in LAI in Guangzhou and Zengcheng due to the conversion of woodland to grassland increased surface albedo by a rate of 0.34 ( p < 0.01) in Guangzhou and 0.27 ( p < 0.01) in Zengcheng. Because such land use conversion was not significant in Taishan, Zhongshan, and Shenzhen, the resultant albedo change was not significant in these areas. Overall, the decline of the RET in Fig. 7. Variation of (a) LAI and albedo and (b) the correlation coefficient between LAI and albedo during 2000 -2016. *, **, and *** are represented as p < 0.1, p < 0.05, and p < 0.01, respectively. "0" means the trend of LAI is zero. The slope units of LAI and albedo are yr −1 and % yr −1 , respectively.
Guangzhou and Zengcheng was closely related to the albedo increase resulting from the LAI decrease under LUCC.

Mechanisms for climate change to affect RET a. Guangzhou and Zengcheng
The SVP increased significantly by a rate of 0.010 kPa yr −1 ( p < 0.05) in Guangzhou and 0.008 kPa yr −1 ( p < 0.1) in Zengcheng owing to the obvious increase of T max and T min (Eq. 8). However, the increase in RH was not so large that actual vapor pressure (AVP) decreased by a rate of 0.004 kPa yr −1 ( p < 0.05) in Guangzhou and 0.003 kPa yr −1 ( p < 0.1) in Zengcheng according to Eq.(7).
According to Eq.(4), the increase of T max and T min and the decrease of the AVP will cause R nl to increase. However, because R nl decreased by a rate of 0.007 MJ m −2 day −1 ( p < 0.05) in Guangzhou and 0.007 MJ m −2 day −1 ( p < 0.01) in Zengcheng, this indicates that the decrease of R s was the only reason for the decrease of R nl . In fact, R s decreased significantly by a rate of 0.033 MJ m −2 day −1 ( p < 0.01) in Guangzhou and 0.027 MJ m −2 day −1 ( p < 0.05) in Zengcheng because of the obvious decrease in SD (Table 3). Likewise, according to Eq.(3), the increase in albedo (Fig.  7) and the decrease in R s can cause R ns to decrease significantly by a rate of 0.025 MJ m −2 day −1 ( p < 0.05) in Guangzhou and 0.021 MJ m −2 day −1 ( p < 0.1) in Zengcheng.
This suggested that both incoming and outgoing radiation were decreasing. However, according to radiation balance (Eq. 2), because the incoming radiation decreased more than that of the outgoing radiation, the net radiation was thus diminishing by a rate of 0.018 MJ m −2 day −1 ( p < 0.05) in Guangzhou and 0.014 MJ m −2 day −1 ( p < 0.05) in Zengcheng. In addition, the stabilization of air mobility and the weakening of energy exchange due to the significant reduction in WS reduced the RET in Guangzhou and Zengcheng, respectively.

b. Taishan and Zhongshan
On the one hand, the mechanisms triggering R n changes in Taishan and Zhongshan were similar to that in Guangzhou and Zengcheng: the reduction of R s due to the obvious decrease in SD resulted in the reduction of R n . However, because the increase in SD was much smaller than that in Guangzhou and Zengcheng, the reduction of R n was relatively small (slope: −0.003 MJ m −2 day −1 [ p < 0.1] in Taishan and −0.005 MJ m −2 day −1 [ p < 0.1] in Zhongshan), and the consequent effects on the RET were not as remarkable as in Guangzhou and Zengcheng. On the other hand, the high air humidity reduces the ET demand in humid tropical regions (Allen et al. 1998). Because the RH decreased significantly (Table 3) and the air became drier, namely, the vapor pressure deficit (VPD) became larger (slope: 0.019 kPa yr −1 [ p < 0.05] in Taishan and 0.026 kPa yr −1 [ p < 0.05] in Zhongshan), more water can be evaporated into the atmosphere (Allen et al. 1998). This played an important role in promoting the RET.

c. Shenzhen
The mechanisms for triggering RET changes in Shenzhen was characterized by the integration of the above two mechanisms. On the one hand, like Guangzhou and Zengcheng, the sharp decrease in the SD decreased R n by a rate of 0.019 MJ m −2 day −1 ( p < 0.05) in Shenzhen. On the other hand, similar to Taishan and Zhongshan, the large increase in T min and high C L of T max and T min led to the decrease in RH and hence in the VPD by a rate of 0.025 kPa yr −1 ( p < 0.05) in Shenzhen. Although the absolute value of C L of WS was the 2nd smallest among the five climatic factors in Shenzhen, however, it was larger than those in the other four stations, and the absolute value of C L for WS at Shenzhen was ranked 14th among all 25 values listed in Table 3. Therefore, the decrease in WS caused the RET to decrease significantly in Shenzhen. In addition, the decrease in R n and WS and the increase in VPD counteracted each other to affect the RET in Shenzhen. The absolute value of C L for VPD was larger than WS and R n (Table 4), indicating that the increase of VPD was the main factor responsible for the RET increase in Shenzhen. As a whole, the increased RET was induced by the decreased RH, and hence the increased VPD resulted from the increase in temperature.

Conclusions
In this study, we quantified the variations of daily RET in the PRD during 1960 -2016 on the basis of the FAO Penman-Monteith equation and examined the dominant factors triggering RET variations. Our results indicated the following. (1) The RET in the PRD had an overall increasing trend during 1960 -2016 along with the increase of construction land. The more the construction land increased, the greater the RET increased.
(2) The noticeable decline of RET in Guangzhou and Zengcheng was closely related to surface albedo increase due to the conversion of woodland to grassland.
(3) The dominant factors responsible for RET variations varied across space inside the PRD. SD and WS were the dominant factors in decreasing the RET in Guangzhou and Zengcheng through decreasing R n and energy exchange. In contrast, T max , T min , and RH were the three main factors increasing the RET in Taishan, Zhongshan, and Shenzhen owing to their effects on the VPD. Based on the mechanism analysis for land use change and climate change to affect the RET, we can see that it will be more conducive to reducing RET by increasing the area of grassland and decreasing the area of construction land. The in-depth mechanisms for the changes of the natural environment caused by land use change need to be analyzed in future work based on this study.