Spatiotemporal Variability of Actual Evapotranspiration and the Dominant Climatic Factors in the Pearl River Basin, China

Abstract

Evapotranspiration is a vital component of the land surface process, thus, a more accurate estimate of evapotranspiration is of great significance to agricultural production, research on climate change, and other activities. In order to explore the spatiotemporal variation of evapotranspiration under global climate change in the Pearl River Basin (PRB), in China, this study conducted a simulation of actual evapotranspiration (ETa) during 1960–2014 based on the variable infiltration capacity (VIC) model with a high spatial resolution of 0.05°. The nonparametric Mann–Kendall (M–K) test and partial correlation analysis were used to examine the trends of ETa. The dominant climatic factors impacting on ETa were also examined. The results reveal that the annual ETa across the whole basin exhibited a slight but not significant increasing trend during the 1960–2014 period, whereas a significant decreasing trend was found during the 1960–1992 period. At the seasonal scale, the ETa showed a significant upward trend in summer and a significant downward trend in autumn. At the spatial scale, the ETa generally showed a decreasing, but not significant, trend in the middle and upper stream of the PRB, while in the downstream areas, especially in the Pearl River Delta and Dongjiang River Basin, it exhibited a significant increasing trend. The variation of the ETa was mainly associated with sunshine hours and average air pressure. The negative trend of the ETa in the PRB before 1992 may be due to the significant decrease in sunshine hours, while the increasing trend of the ETa after 1992 may be due to the recovery of sunshine hours and the significant decrease of air pressure. Additionally, we found that the “paradox” phenomenon detected by ETa mainly existed in the middle-upper area of the PRB during the period of 1960–1992.

1. Introduction

Evapotranspiration (ET) is a vital component of the hydrological cycle over land surface and energy balances [1,2,3]. Over 60% of precipitation on land is returned to the atmosphere via ET, thus, ET greatly influences the water availability at the land surface [4]. As a key energy flux, ET is also associated with more than half of the radiation energy absorbed by land surfaces [5], thereby helping to adjust the land surface temperature, with crucial implications for regional and global warming [6,7]. Accordingly, ET is not only an essential element of water balance, but also plays a key role in the energy budget of the earth-atmospheric system. The discussion on the variability of ET is, thus, a hot topic recently under the background of global climate change.

Accurate estimation of ET is essential in various fields, such as analysis of drought, water resources utilization, climate change, ecological protection, and agricultural production [8,9,10,11,12,13]. However, ET is a complex process influenced by land surface factors (e.g., soil properties and vegetation) and climatic parameters (e.g., temperature) [10,14]. Therefore, precise and quantitative estimation of ET remains a major problem [15] and the terrestrial ET has long been a hot topic in meteorology, geography, and hydrology worldwide. ET can be generally expressed by pan evaporation (PE), potential ET (PET), reference ET (ET₀), and actual ET (ETa). PE refers to estimates of ET observed from an evaporation pan. PET is the amount of water transpired in a given time by a short green crop, completely shading the ground, of uniform height and with adequate water status in the soil profile. ET₀ was proposed to make PET more controlled and standard and it can be determined by the Penman–Monteith method recommended by the Food and Agriculture organization of the United Nations (FAO), i.e., ET₀ can be regarded as a special case of PET. ETa is used to designate the ET under natural conditions with often inadequate water supply when it cannot be measured or computed directly [16]. To date, studies on PET have been relatively abundant and many formulas and methods for estimating PET have been widely used [16,17,18]. However, PET describes the circumstance that water is abundantly available at the underlying surface, usually significantly larger than the ETa [19,20,21]. In this sense, ETa would more realistically represent ET and provide a better understanding of terrestrial ET under natural conditions.

Relative to PET and ET₀, the research on ETa and its applications are relatively few and limited due to the complexity of ETa process [21,22]. Methods to determine ETa mainly depend on indirect approaches, including conversion from PET, instrument measurements, remote sensing observations, water balance systems, simple hydrological models, and distributed hydrological models [6,15]. Among these indirect approaches, state-of-the-art distributed hydrological models have successfully improved the completeness of the process representation and are regarded as a promising approach for estimating ETa [23]. For example, the variable infiltration capacity (VIC) model has been validated in a wide range of climatic regions and achieved good results on ETa estimation, including humid and arid areas [24,25,26,27,28,29,30]. As this model can comprehensively consider the complicated interactions and hydrological feedback among atmosphere, soil and vegetation, we are able to simulate the land surface ET process more realistically when the underlying surface and meteorological data are adequate [23].

In addition, there is an interesting phenomenon in the field of geography and hydrology related to ET change, referred to as “evaporation paradox”. As the average global temperature rises due to global warming, it is generally expected that evaporation from terrestrial water bodies will increase [31]. However, contrary to widely held expectations, the observed pan evaporation across the world shows a decreasing trend, and such unexpected phenomenon is called the “evaporation paradox” [7,32,33,34,35]. The evaporation paradox was also found in some parts of China and its cause and influential factors were also discussed [36,37,38,39,40,41]. However, such studies were mainly based on PE or ET₀ data, and the use of ETa data to examine the paradox is still scant. The PE, ET₀ and ETa can reflect the evaporation, but the spatiotemporal features of these indices may present some difference due to the diverse definitions, assumptions, and calculations. In this instance, whether ETa indeed increases or the occurrence times of the ETa “paradox” are consistent with that of PE or ET₀ are rarely reported and then is worth exploring [42,43].

Although some studies on ETa have been conducted in China, most of them were in the North China Plain, Yangtze River Basin, Haihe River Basin, etc., [22,44,45,46,47] and studies on ETa in the Pearl River Basin (PRB), the largest and most vital basin in Southern China, are still lacking. Additionally, studies on a national spatial scale have preliminarily found that the PE in the PRB shows a declining trend as the air temperature rises continuously [37,48,49]. However, as mentioned before, whether the paradox still exists in the PRB when using ETa estimation requires further verification. Moreover, what the spatiotemporal features and variability of ETa will present in PRB under global warming, and which factors dominate the trends of ETa and the “paradox” are still unknown.

Therefore, this study uses the VIC model for simulating the hydrological cycle in the PRB, and its objectives include: (1) to analyze the spatiotemporal variation of simulated ETa from 1960–2014 under global warming; (2) to explore the correlations between ETa and related influencing climatic factors; (3) to clarify whether a “paradox” phenomenon still exists when using ETa estimation. This study can serve as reference for further studies about ETa and “evaporation paradox” under global warming, and contribute to agricultural management and production, prevention, and reduction of droughts, and overall ecological restoration in the PRB.

2. Data and Methods

2.1. Study Area

The PRB (Figure 1) is mainly located in Southern China and covers parts of the Yunan, Guizhou, Guangxi, Guangdong, Hunan, and Jiangxi provinces, as well as Vietnam (102°14′–115°53′ E, 21.6–26.9° N). Its elevation ranges from −18 to 2795 m, and its drainage area is approximately 45.4 × 10⁴ km². The Pearl River is the second largest river (according to streamflow magnitude) in China with an annual streamflow of 34.9 × 10¹⁰ m³/a. The Xijiang River, Beijiang River and Dongjiang River are three major tributaries [50,51,52]. Located in the subtropical humid monsoon zone, the PRB has an annual mean air temperature varying between 14 and 22 °C and the annual mean precipitation ranges from 1200 to 2200 mm. The Pearl River Delta, including cities like Hong Kong, Guangzhou, and Shenzhen, is an important economic region in China, accounting for approximately 10% of the Chinese Gross Domestic Product (GDP) [48,53]. The upper region of the PRB mainly features karst landform and has suffered from frequent drought events [48,54,55,56].

2.2. Data

Basic input data for the VIC model include the meteorological forcing data and underlying surface data. The underlying surface data consist of soil data and vegetation coverage data. In this study, we collected daily meteorological data (including precipitation, daily maximum and minimum temperature, wind speed, sunshine hours, shortwave radiation, and relative humidity (for details please see

Table 1. Meteorological and hydrological data used in this study.

Data Category	Data Type	Source	Unit
Meteorological	Precipitation	The Chinese Meteorological Administration	mm
	Max daily air temperature		°C
	Min daily air temperature		°C
	Wind speed		m/s
	Sunshine hours		h
	Shortwave radiation		W/m2
	Relative humidity		%
Hydrological	Observed streamflow	Hydrological stations	m3/s

) spanning the period from 1959 to 2014. These data, from 63 meteorological stations in or around the PRB, were provided by the Chinese Meteorological Administration (CMA), an organization taking charge of observing, collecting, compiling, and releasing meteorological data of high quality in China. Monthly streamflow observations were collected from three hydrological stations (i.e., Gaoyao, Shijiao, and Boluo) in the PRB. The three hydrological stations correspond to the three main streams of the PRB respectively, i.e., the Xijiang River, Beijiang River, and Dongjiang River. Then, the meteorological forcing data are interpolated to produce the same spatial scale as the VIC model using the Kriging interpolation method. The soil data were obtained from the global soil classification and texture dataset based on the Harmonized World Soil Database (HWSD) published by the FAO [57]. The soil hydraulic characteristic parameters required by the VIC model were calculated through the soil conversion formula summarized by Saxton et al. [58]. The vegetation coverage data were extracted from the global land cover data with a resolution of 1 km × 1 km [59]. The digital elevation model (DEM) data at a resolution of 90 m were derived from the Shuttle Radar Topography Mission (SRTM). The above data were checked by quality control procedures to ensure their reliability [60,61] and then were converted into the input data files through ArcGIS 9.3 and the R programming language. Since the VIC model can only output runoff yield data of each grid, it is also necessary to generate the streamflow values at the outlets of the basins by coupling a routing model developed by Dag Lohmann et al. [62] for calibration.

At the spatial scale, the study area was divided into 13 sub-regions according to the three-level zoning of water resources and the administrative division of the region (see Figure 1). The grid data output by the VIC model are aggregated into each sub-region for analysis. At the temporal scale, the data were aggregated into monthly, seasonal and annual scales for analysis, in which the annual data were generated through the natural year method. Seasonal data were calculated according to the practice of dividing a year into spring (from March to May), summer (from June to August), autumn (from September to November), and winter (from December to February).

2.3. VIC Hydrological Model

The VIC model, a macro-scale distributed hydrological model based on grid-cell representations, was developed by Liang et al. [63] and maintained by three universities, including Washington University (Saint Louis, MO, USA). The VIC model can be operated in both water and energy balance modes in each grid cell. Additionally, various elements in a grid cell, including different vegetation types, air, soil, and terrain, are fully considered by the model. Thus, this model makes up for the deficiency of the traditional hydrological model in energy exchange and has great regional applicability. The model simplifies the soil structure in the grid and divides it vertically into several layers, and the VIC model with three layers (VIC-3L) is the most widely used mode. The research results of previous studies indicate that the VIC-3L model is applicable in southern China [25,26,27]. In this study, the VIC-3L model was established in the PRB with each grid at 0.05° × 0.05° spatial resolution. The ET simulation consisted of three parts: canopy evaporation, soil water evaporation, and vegetation transpiration. Integrating water and energy balance, the calculation of the ETa includes energy and water exchange processes of evaporation, infiltration, and runoff. Since, at present, the VIC model does not support sunshine hours as an input parameter and this study lacks thermodynamic data, such as short-wave radiation. Therefore, we used the Angstrom formula [64] to estimate the short-wave radiation data required by the VIC model based on sunshine hours, air temperature, and water vapor pressure; then we investigated the simulation of the ETa.

Three performance indices, the Nash–Sutcliffe coefficient of efficiency (NSCE), correlation coefficient (R) and relative error (RE), were used to assess the hydrological simulation performance of the model. The NSCE represents the agreement between the simulated and observed values for the time series of streamflow. The value of the NSCE ranges from −∞ to 1. A value of 1.0 represents a perfect match between simulated and observed streamflow:

$$NSCE=\frac{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-S_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(O_i-\overline{O}\right)}^2}$$

(1)

where $O_i$ and $S_i$ are the observed and simulated values of streamflow time series, and $\overline{O}$ is the mean value of $O_i$.

In this study, during the calibration period (1960–1980), the VIC model was calibrated at the three hydrological stations (i.e., Gaoyao, Shijiao, and Boluo) by means of optimizing the NSCE of the simulated streamflow. Then, the calibrated model was run in the validation period of 1981–2000 to validate the calibrated results for evaluating the applicability of the model.

2.4. Methodology

2.4.1. Mann–Kendall Trend Test (M–K)

To analyze the temporal variability of the ETa and related influential factors in the PRB, in this study we used the non-parametric Mann–Kendall (M–K) trend test [65,66]. The M–K test can effectively detect trends in a time series and does not require data to obey normal distribution. Moreover, this method is not easily affected by outliers. Thus, the M–K test is strongly recommended by the World Meteorological Organization (WMO) and has been widely used in trend analysis of meteorological and hydrological data [65].

The advantage of the M–K analysis is that the series does not require a certain sample distribution, thereby avoiding the potential interference of outliers. However, the use of M–K test requires that the sample data are serially independent. Climatic and hydrologic series typically display serial correlation, making them not well suited for the analysis [67]. To eliminate the effect of a serial correlation on the M–K test, Yue and Wang proposed a modified pre-whitening procedure called trend-free pre-whitening (TFPW) [68]. In this study, we used TFPW to eliminate the effect of serial correlation. Then, two crucial indicators of the trends of time series can be obtained through the M–K test. The size of the statistical variable M is used to assess the significance of the trends based on a certain significance level. Negative values of M represent downward trends, while positive values of M indicate increasing trends. The slope reflects the average rate of change of the time series. In addition, the M–K test can also be used for detecting the abrupt change points [62]. In this study, we used the forward sequence of M to establish a curve called UF. Then, the curve called UB was produced based on the reverse sequence. The intersections of the two curves were the points where the abrupt changes occurred [69].

2.4.2. Partial Correlation Analysis

In order to investigate the causes of ETa changes, the partial correlation method [70,71] was used in this study to evaluate the degree of correlation between each climatic factor and the ETa in the basin. Partial correlation analysis can exclude the influence of the other variables and describe the unique correlation between the two considered variables. In multivariate correlation analysis, if the interaction among variables is so complex that the simple correlation coefficient cannot truly reflect the correlation between variables, partial correlation analysis is a suitable method. To date, partial correlation coefficient has been widely used in the field of hydrology and water resources [44,72,73]. In this study, the partial correlation coefficient is calculated using the “corpcor” package (https://cran.r-project.org/web/packages/corpcor/) in the R platform [74].

Partial correlation coefficient can be iteratively calculated in terms of the following formula:

$$r_{XY|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_n}=\frac{r_{XY|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_{n-1}}-r_{X{\mathrm{Z}}_n|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_{n-1}}·r_{Y{\mathrm{Z}}_n|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_{n-1}}}{\sqrt{1-r_{X{\mathrm{Z}}_n|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_{n-1}}^2}\sqrt{1-r_{Y{\mathrm{Z}}_n|{\mathrm{Z}}_1\cdots {\mathrm{Z}}_{n-1}}^2}}$$

(2)

where $r_{XY|Z_1\cdots Z_n}$ is an N-th order partial correlation coefficient between X and Y that removes $Z_1\cdots Z_n$. The N-th order partial correlation coefficient is to remove N variables. The 0-th order partial correlation coefficient is the simple correlation coefficient and the calculation formula is as follows:

$${\mathrm{r}}_{\mathrm{XY}}=\frac{{{\displaystyle \sum }}_{i=1}^n\left(X-\overline{X}\right)\left(Y-\overline{Y}\right)}{{{\displaystyle \sum }}_{i=1}^n{\left(X-\overline{X}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(Y-\overline{Y}\right)}^2}$$

(3)

To perform the significance test of the partial correlation coefficient, the Student’s t test method was used. The t-statistics are calculated as follows:

$$\mathrm{t}=\frac{r\sqrt{n-q-2}}{\sqrt{1-r^2}}$$

(4)

where r is the partial correlation coefficient, n is the sample number of the factor sequence, q is the order of r, and the statistic t satisfies the t distribution with n-q-2 degree of freedom. By calculating the statistical t and taking the bilateral probability of its t distribution, we can evaluate its significance by comparing it with the given significance level.

3. Results

3.1. Spatiotemporal Characteristics of ETa

3.1.1. Assessment of the VIC Model Performance

As shown in Figure 2, the monthly NSCEs of the three hydrological stations (Boluo, Gaoyao, and Shijiao) are all higher than 0.85 in both calibration (1960–1980) and validation (1981–2000) periods; all the correlation coefficients are higher than 0.95 with the absolute values of relative errors less than 15%. The results indicate that the VIC model performs satisfactorily for the hydrological simulation and maintains a reasonable water balance [25], demonstrating that the calibrated model can be used for the estimation of ETa in the PRB.

3.1.2. Annual ETa in the PRB

The daily ETa data of all grids in the entire basin exported by the VIC model from 1960 to 2014 were statistically analyzed at the spatial and temporal scales. From the statistical analyses, we obtained the annual sequence of the ETa over 55 years in the basin (Figure 3) and the spatial distribution of the annual average ETa (Figure 4).

The data shown in Figure 3 reveal that the annual average ETa in the PRB was 626.1 mm/a and varied from 585.5 to 679.3 mm over the 55 years. In terms of temporal changes, the annual ETa first decreased and then increased during the 55 years. Before the 1990s, the annual ETa generally showed a downward trend and reached the minimum value of 585.5mm in 1992. Subsequently, the annual ETa increased to a maximum value of 679.3 mm in 1997.

The prominent spatial heterogeneity in the spatial distribution of the ETa is displayed in Figure 4. The annual average ETa in the PRB was within the range of 368.4–951.1 mm. The high ETa values were mainly distributed around the estuary of the Pearl River Delta, and downstream of the Zuojiang, Youjiang, and Hongshui rivers as well. The low values mainly distributed in western region, including the sub-region of the Nanpan River, Beipan River, and the upper reaches of the Youjiang as well as Hongshui rivers.

3.1.3. Spatial-Temporal Variation of the ETa

The change trends of the ETa were analyzed by using the M–K trend test method at the monthly, seasonal, and yearly scales (

Table 2. M–K values for the ETa in the Pearl River Basin and its sub-regions at the monthly scale.

Region	January	February	March	April	May	June	July	August	September	October	November	December
1	0.25	0.43	0.41	−0.18	0.57	0.82	0.01	−0.72	−1.6	−1.81△	−1.74△	−1.95△
2	−1.44	0.31	0.52	−0.83	−0.24	1.2	1.33	−0.73	−1.43	−1.47	−1.36	−2.45▲
3	−1.6	0.85	0.7	−1.87△	−0.38	1.75△	0.57	−1.69△	−0.91	−1.63	−1.71△	−1.34
4	−1.81△	0.94	0.72	−1.46	−0.18	3.13★	1.4	−1.1	−1.55	−1.28	−2.19▲	−1.28
5	−1.53	0.38	0.3	−0.83	−0.24	2.34▲	1.07	0.23	−1.91△	−2.03▲	−1.17	−1.44
6	−1.27	−0.14	0.17	−0.91	−1.1	2.98★	2.3▲	2.17▲	−0.86	−2.05▲	−1.66△	−0.56
7	−1.71△	0.59	−0.17	−0.2	−0.05	2.71★	1.62	−1.53	−1.66△	−2.19▲	−1.44	−0.21
8	−0.53	0.31	1.1	0.81	1.4	3.36★	3.33★	1.07	0.69	−1.11	−1.39	−0.17
9	−0.08	2.5▲	0.49	1.31	0.73	2.93★	1.89△	0.02	0.21	−2.33▲	−1.49	−0.08
10	−0.56	0.78	0.52	1.44	1.37	3.81★	2.01▲	1.04	1.5	−1.36	−0.46	0.41
11	0.2	1.04	1.21	1.33	1.23	2.46▲	3.3★	3.13★	2.72★	−1.39	−1.02	0.28
12	0.47	0.63	0.68	0.59	2.46▲	3.42★	4.28★	3★	2.81★	−0.53	−1.07	1.02
13	0.5	−0.05	0.86	1.14	3.48★	4.45★	4.07★	3.46★	3.68★	−0.97	−0.69	1.2
PRB	−0.86	0.66	0.88	−0.52	0.82	3.87★	2.3▲	0.34	−0.79	−1.89△	−1.62	−0.37

Notes: the names of the sub-regions refer to those in Figure 1; △, ▲, and ★ represent significance level of 0.1 (1.65), 0.05 (1.96), and 0.01 (2.58), respectively.

and

Table 3. M–K values for the ETa in the Pearl River Basin and its sub regions at the seasonal and annual scale.

Region	Spring	Summer	Autumn	Winter	Year
1	−0.05	0.38	−2.53▲	−0.16	−0.94
2	−0.24	0.24	−2.42▲	−0.68	−0.98
3	−0.24	0.81	−2.36▲	−0.96	−1.55
4	−0.85	1.52	−2.69★	−1.07	−1.47
5	−0.3	2.11▲	−2.53▲	−1.14	−0.54
6	−0.72	3.72★	−2.4▲	−1.32	−0.76
7	−0.04	1.04	−2.55▲	−0.9	−0.83
8	2.04▲	3.46★	−1.59	−0.43	1.52
9	1.62	2.11▲	−2.07▲	0.74	1.33
10	1.84△	3.71★	−1.11	0.1	2.56▲
11	1.49	3.7★	−0.54	0.16	2.42▲
12	1.65	4.78★	−0.15	0.72	3.03★
13	2.42▲	5.02★	−0.08	0.89	3.58★
PRB	0.43	3.16★	−2.21▲	−0.69	0.21

Notes: the name of the sub-regions refer to those in Figure 1; △, ▲, and ★ represent significance level of 0.1 (1.65), 0.05 (1.96), and 0.01 (2.58), respectively.

). The results of the M–K test show that the M–K value of the annual ETa in the PRB during the period of 1960–2014 is 0.21, which means that the ETa generally showed a slight, but not significant increasing, trend over the past 55 years. The average change rate was approximately 0.042 mm/a, that is, it increased by about 2.34 mm in the past 55 years. At the seasonal scale, the ETa in the PRB showed a significant upward trend in summer (p < 0.05) and a significant downward trend in autumn (p < 0.01). At the monthly scale, the ETa across the whole basin significantly decreased in October (p < 0.1) and significantly increased in June (p < 0.05) and July (p < 0.01). Additionally, at the sub-region scale, the ETa in the upstream region displayed a more obvious downward trend, while the ETa in the downstream region exhibited a more obvious upward trend.

The data presented in Figure 5 reveal that the UF curve only comes in contact with the UB curve in 1973 and there was no significant intersection, which means that there is no obvious abrupt change over the 55 years in the PRB. According to the shape of the UF curve, the ETa roughly showed an obvious downward trend in the period of 1960–1992, with a M–K value of about −1.89 in the year 1992, passing the significant test at the 0.1 significance level (−1.65) and close to 0.05 significance level (−1.96). This means that from 1960 to 1992 the ETa had a decreasing trend in the PRB although there is no obvious abrupt change in the 55-year series; subsequently, the UF curve increased, indicating that the ETa in the basin gradually turned to an upward trend in the 1990s.

At the spatial scale (Figure 6a), the ETa in the upper stream of the PRB generally showed no significant increasing or decreasing trend. Some regions, including the sub-region of the Liujiang River, Zuojiang River and the upstream of the Hongshui River, exhibited a significant decreasing trend at the 0.05 significance level. Additionally, in the downstream region of the PRB, especially the Pearl River Delta, the Dongjiang River basin and the downstream of the Beijiang River, the ETa displayed a significant increasing trend at the 0.01 significance level. The change rate of the ETa in the PRB is shown in Figure 6b. The largest reduction of the ETa occurred in the south part of the PRB and the area with the largest increase was mainly located in the Pearl River Delta.

3.2. Relationships between the ETa and Climatic Factors

3.2.1. Partial Correlation Coefficient

In order to investigate the causes of the spatial distribution and temporal variation trend of the ETa in the PRB, a partial correlation coefficient analysis was conducted to evaluate the influence of various climatic factors on the ETa. With the aim of removing the interference of other factors as much as possible, ten climatic factors used in the relevant studies were selected to calculate the 9th order partial correlation coefficient of the ETa, and then we conducted a significance test by the Student’s test. The ten factors included: water quantity factors, such as precipitation (PREC), relative humidity (RH), water vapor pressure (VP); thermal factors, such as average temperature (TEMP), daily minimum and maximum temperature (TMIN, TMAX), daily temperature range (TMRG), sunshine hours (SH); dynamic factors, such as average pressure (PRESS), and mean wind speed (WIND). All the ten climatic factors are statistically independent from each other.

Accordingly, from the foregoing, we can learn that ETa in the PRB exhibited two different trends, namely before and after 1992. After analyzing the time series of climatic factors, we found their trends also changed around the year 1992. Thus, the ETa and climatic factors in the basin were divided into two periods (i.e., 1960–1992 and 1992–2014) for the convenience of calculation. The partial correlation coefficients of the influence of climatic factors on the ETa at different time scales in the two time periods are shown in

Table 4. Partial correlation coefficients between climatic factors and the ETa during the two periods of time in the Pearl River Basin.

	1960~1992						1992~2014
	Monthly	Spring	Summer	Autumn	Winter	Yearly	Monthly	Spring	Summer	Autumn	Winter	Yearly
PREC	0.55★	0.66★	0.09	0.65★	0.39▲	0.44▲	0.48★	0.89★	0.74★	0.80★	0.70★	0.70★
RH	–0.06	0.11	0.24	0.09	0.27	–0.04	–0.11△	0.42	–0.36	0.05	0.23	0.24
VP	0.38★	0.05	–0.14	0.15	0.05	0.53★	0.35★	–0.54▲	0.34	0.02	–0.03	–0.25
TEMP	0.18★	–0.16	–0.02	–0.20	–0.03	–0.47▲	0.18★	–0.60▲	–0.78★	–0.15	0.08	–0.59▲
TMAX	0.05	0.34△	0.04	0.44▲	0.27	–0.06	0.06	0.00	0.76★	–0.04	0.32	–0.10
TMIN	–0.06	–0.34△	–0.04	–0.44▲	–0.27	0.06	–0.06	–0.00	–0.76★	0.04	–0.32	0.10
TMRG	–0.06	–0.34△	–0.04	–0.44▲	–0.27	0.06	–0.06	–0.00	–0.76★	0.04	–0.32	0.10
SH	0.24★	0.38△	0.14	0.49▲	0.02	–0.25	0.20★	0.83★	0.26	0.19	0.68★	0.66★
PRESS	–0.21★	0.46▲	–0.09	0.20	0.31	–0.06	–0.35★	0.82★	0.17	0.14	0.18	0.48△
WIND	0.06	–0.35△	0.47▲	0.51★	0.09	–0.13	–0.01	–0.69★	0.28	–0.11	0.21	0.05

Notes: △, ▲, and ★ represent significance level at 0.1, 0.05, and 0.01, respectively.

.

During the 1960–1992 period, except for PREC, the VP and TEMP showed a significant correlation with the ETa both at the yearly scale, and the strong correlation of VP with ETa passed the significance test at the 0.01 level. At the monthly scale, the ETa showed significant correlations with the SH and PRESS. In addition, thermodynamic factors, such as WIND, TMAX, TMIN, and TMRG showed a high partial correlation with the ETa at the monthly scale in spring and summer.

During the 1992–2014 period, the SH had the most notable impact on the ETa (except for PREC), exhibiting a very significant correlation with a significant level of 0.01 in spring, summer and at the yearly scale. TEMP also had a very significant correlation with the ETa in spring, summer, and at the monthly and yearly scales. According to the data in Table 4, the partial correlation coefficients of the ETa and different climatic factors showed significant differences in the different time scales, which is consistent with the findings reported by Li et al. [44].

3.2.2. Temporal Variation of the Ten Climatic Factors

The change trends of the ten climatic factors in the PRB during the period of 1960–2014 are shown in Figure 7. The regression analysis results for each of the factors in the two periods of 1960–1992 and 1992–2014, revealed that factors including RH, TMAX, SH, PRESS, and WIND showed different change trends before and after 1992.

From 1960 to 1992, PREC, TEMP and VP, significantly correlated with the ETa at the yearly scale, but none showed obvious declining or rising trend. SH, which exhibited a very significant correlation with the ETa at the monthly scale, showed a significant declining trend. Thus, it can be speculated that the negative trend of the ETa in the PRB before 1992 may be mainly due to the significant decrease of SH. From 1992 to 2014, at the monthly and yearly scales, PRESS, TEMP, and SH showed significant correlations with the ETa. Compared with the period of 1960–1992, PRESS turned from a no apparent change trend to a significant decreasing trend; TEMP rose more significantly; SH turned from a significant decline to no obvious change trend, which indicate that the upward trend of the ETa in the PRB after 1992 may be due to the significant decrease of PRESS and the less sharp decline of SH.

In order to further investigate the main climatic factors affecting the ETa in each sub-region of the PRB, the partial correlation coefficients between the ETa and different climatic factors in each sub-region were calculated and analyzed at the yearly scale. The results, which are shown in

Table 5. Partial correlation coefficients between climatic factors and the ETa at the sub-region scale in the Pearl River Basin.

	PREC	RH	VP	TEMP	TMAX	TMIN	TMRG	SH	PRESS	WIND
1	0.62★	0.14	−0.04	−0.02	0.02	−0.02	−0.02	0.06	−0.21	0.03
2	0.69★	0.00	0.07	−0.06	0.26△	−0.26△	−0.26△	0.33▲	−0.12	0.17
3	0.18	0.10	0.23	−0.22	0.09	−0.08	−0.09	0.07	0.02	0.14
4	0.07	−0.10	0.26△	−0.42★	0.14	−0.13	−0.13	−0.16	0.14	0.15
5	0.35▲	0.02	0.02	−0.16	0.08	−0.07	−0.08	0.22	−0.02	0.03
6	0.64★	−0.28△	0.29▲	−0.24△	−0.06	0.06	0.06	0.30▲	−0.13	0.12
7	0.45★	−0.34▲	0.38★	−0.60★	−0.15	0.15	0.15	−0.25△	0.08	0.62★
8	0.63★	−0.13	0.09	−0.38★	−0.04	0.04	0.04	−0.07	−0.02	0.40★
9	0.51★	−0.27△	0.27△	−0.47★	0.01	−0.01	−0.01	−0.10	−0.08	0.24
10	0.53★	−0.31▲	0.30▲	−0.41★	0.00	0.00	0.00	0.02	−0.07	0.39★
11	0.48★	0.03	0.05	−0.14	0.08	−0.08	−0.08	−0.09	−0.21	0.40★
12	0.49★	0.25△	−0.26△	−0.06	−0.29▲	0.29▲	0.29▲	0.00	−0.14	0.38★
13	0.69★	0.20	−0.21	−0.07	−0.33▲	0.33▲	0.33▲	−0.01	−0.22	0.38★

Notes: the name of sub-regions refer to those in Figure 1; △, ▲, and ★ represent significance level at 0.1, 0.05, and 0.01, respectively.

, revealed that PREC is an important factor affecting the ETa in most regions. In addition, TEMP and VP had great influence on the ETa in the middle and lower reaches of the basin (i.e., in the sub-region of the Liujiang River, Zuojiang River, Yujiang River, Hegui River, etc.). Additionally, WIND also had a significant influence on the ETa in downstream areas, such as the Dongjiang River basin and the Pearl River Delta.

4. Discussion

4.1. Comparison with Related Studies

In the PRB, we found that the ETa showed a significant declining trend during the period of 1960–1992 and then gradually turns into an increasing trend from 1992 onwards. Zhang et al. stated the ETa decreased gradually from east to the west and our study showed a consistency with theirs. They also found decreasing wind speed and decreasing sunshine hours greatly dominated the change of actual evapotranspiration [75], which verified that sunshine duration is a dominant climatic factor of ETa in the PRB. However, only four climatic factors (i.e., temperature, precipitation, wind speed, and sunshine hours) and a shorter time series (1960–2005) were used in their study, causing a few differences with our study. Zhang and Chen found that the ETa in PRB showed an increasing trend during 2000–2014 [76], which is similar with the trend after the year 1992 in our study. However, Li et al. used the advection-aridity (A-A) model to carry out research in the PRB and their results were completely different from those of our study [77,78]. According to Li et al., the ETa in the PRB exhibited a significant decreasing trend from 1960 to 2008, which may be explained by the fact that the A-A model only considers a relatively ideal condition. In fact, according to the Penman–Monteith and Priestley–Taylor methods, the A-A model only roughly considers the vegetation and ignores the underlying surface details, which contributes to its possible low performance [79]. Thus, this indicates that the A-A model may not adequately reflect the actual situation and may be unsuitable for an extremely wet environment [80,81].

In addition, we also compared our results with the ones based on PE, ET₀, and PET. Yao et al. found the annual PET in south China showed a decreasing trend then increased around 1995 [82]; Xing et al. also discovered the ET₀ had a similar turning point in 1995 in Southeast China [40]. This change trend was also verified in Northwest China revealing that PE exhibited an obvious decreasing trend until early 1990s and such a downward trend reversed in 1993 [39]. A similar upward trend of land ET at the global scale was also confirmed by Zhang et al. during the period of 1981–2012 [14]. Their studies revealed that ET had an upward trend at the global scale and are also consistent to a certain extent with our results in the PRB during 1992–2014. Additionally, our study also shared similar characteristics with some studies on a national or global scale [7,14,33,38,39,83,84,85]. Based on the above relevant studies and our study, we could also speculate that the ET in the PRB experienced a turning point in the 1990s, in which it initially decreased and later ETa increased.

The trends of the seasonal ETa in the PRB had obvious regularity that exhibited a significant upward trend in summer and a decreasing trend in autumn. Zhang and Chen reported that the ETa increased in both summer and autumn, which is different from our study [76]. This inconsistency may be explained by the difference in the length of time series between the two studies. On the spatial scale, low ETa was mainly found in the western region with high elevation and lower vegetation coverage, while the eastern and central areas were characterized by high ETa. The results of our study showed that the spatial distribution of low ETa in the PRB was associated with elevation and vegetation types and coverage [76]. Notably, the ETa showed a significant increasing trend in the eastern region, and a decreasing trend in other regions, including the sub-region in the Liujiang River, Zuojiang River and the upstream of the Hongshui River, which showed significant spatial difference due to different topographic and climatic features.

4.2. The Dominant Climatic Factor of ETa in the PRB

In this study, ten climatic factors were chosen to examine their influence on the variation of the ETa in the PRB. According to the aforementioned analysis of the correlations between climatic factors and the ETa, SH can be regarded as the main factor influencing the variation of the ETa in the PRB. Since the decreasing SH had a significant correlation with the decreasing ETa during the period of 1960–1992, we were able to infer that SH was the major factor promoting the decrease of the ETa. The results are consistent with many former related studies [41,45,86,87,88]. For example, one study found that SH was the most important variable influencing PET in Southeast China [86]. Through the analysis of the contribution of various factors to ET₀ variation, Wang et al. found that sunshine duration had the greatest impact on the variation of ET₀ in coastal regions with increasing level of urbanization and large population, including the Southern China region [41] where the downstream areas of the PRB are located. Study in Guangdong Province, which belongs to the PRB, also determined that decreasing sunshine duration was the crucial climatic factor leading to the decreasing PE [87]. To some extent, the reliability of the findings in this study was reinforced by the above studies. Most previous studies concluded that the increasing cloud coverage and aerosol concentration would be the main causes of the decreasing sunlight [32,86,89]. It is widely known that, in recent decades, large amounts of aerosol were released to the atmosphere due to human activities, resulting in decreasing SH. Additionally, high-rise buildings in cities hinder the diffusion of aerosols more or less, thereby intensifying the decrease of SH. Thus, these may be the two reasons that SH was the dominant climatic factor influencing the ETa in the PRB, especially in the Pearl River Delta.

In terms of the previous analysis, we found that the partial correlation coefficients between the ETa and climatic factors showed significant differences in different time scales, which is consistent with the study of Li et al. [44]. Compared with the period of 1960–1992, PRESS turned into a significant decreasing trend; TEMP rose more clearly; SH turned from having a significantly declining trend to no obvious changing trend. Since PRESS had a negative correlation with the ETa [44], we also inferred that the increasing trend of the ETa in the PRB after 1992 may be due to the significant decrease of PRESS and the less sharp decline of SH. In addition, at the sub-regional scale, it is somewhat interesting that WIND had a significant correlation with the ETa in the downstream of the PRB including the Dongjiang River Basin and the Pearl River Delta. Considering the variation trend of WIND, we surmised that to a certain degree WIND may contribute to the variation of the ETa. Meanwhile, TEMP and VP had a great influence on the ETa in the middle and lower reaches of the basin. To explain these results, we reasoned that each region had its unique climatic and topographic condition, which resulted in the various correlated factors of ETa in different areas in the same period of time [44,76]. However, in general, it was indisputable that the variation of the ETa in the PRB was dominated by SH.

In particular, the variation of the ETa in the PRB is a complex process, due to the synthetical effects of related climatic factors. During the period of 1960–1992, even when the temperature rose slightly, the markedly decreasing SH led to the total decrease of the ETa. After 1992, with the more markedly rising temperature, SH decreased less sharply when PRESS fell steeply, ultimately resulting in the ETa turning into an upward trend.

4.3. Does the “Paradox” Exist When Using ETa Estimation?

According to the above analysis, over the past 55 years, the annual ETa showed a significant downward trend during the period of 1960–1992. However, the average temperature in the PRB showed an increasing trend during 1960–1992 (Figure 7) and then we can speculate that the “paradox” phenomenon existed in this period, i.e., the “paradox” phenomenon in the PRB can also be detected by using ETa estimation. Subsequently, the ETa gradually turned to an upward trend later in the 1992–2014 period. At the spatial scale, the ETa in the middle and upper stream of the PRB generally exhibited a non-significant decreasing trend when some of its parts showed a relatively significant decreasing trend. In the downstream areas, the ETa showed as a significant increasing trend. According to the above-mentioned results, in general, the “paradox” existed in the variation of ETa during the period of 1960–1992 in the middle-upper area of the PRB and the decline in SH mainly controlled the decrease in the ETa. Furthermore, an obvious decline of SH before the early 1990s in our study was associated with the “global dimming” and the recovery of SH after the 1990s roughly coincided with the “global brightening” [90], consistent with some previous studies [39,40]. The phenomenon of “global dimming and brightening” described by Wild et al. suggested that the declining global solar radiation during the period of 1960–1990s did not continue into the following years [90]. Some researchers concluded that cloud coverage and aerosol, which were in a way influenced by human activities, were the most likely factors to explain the “dimming and brightening”, and the strong aerosol effects of rising air pollution on the decline of solar radiation in China have also been confirmed [40,91]. These ideas could to some extent explain the “paradox” detected by ETa in the PRB, suggesting that we could use ETa to characterize the “evaporation paradox”.

Undoubtedly, the variation of the ETa resulted from synthetical effects of different climatic factors. Although temperature is always considered as a key driver of the ETa changes, it is neither the only factor nor the most decisive factor influencing the variation of the ETa in the PRB. In reality, each region has its own decisive factors. The so called “evaporation paradox” simply associated the ETa with temperature and neglected the effects of other related climatic factors and human activities, which may lead to increasing or decreasing ET. Therefore, to gain a deeper understanding of the “evaporation paradox”, further studies should fully consider the spatiotemporal variation of different influential factors, and focus on how their changes result in the spatiotemporal variations of ET.

Although this study analyzed the spatiotemporal variability of ETa and identified the dominant climatic factors and finally achieved some results, there are inevitably many remaining deficiencies. For example, as the VIC model was used to estimate the ETa and assess the impact of different climatic factors on ETa, the results will be more reliable if a sensitivity and uncertainty analysis are carried out [92,93]. The driving mechanism of climatic factors on ETa in the PRB is still unclear and needs further discussion. Hopefully, these deficiencies could be further discussed in the future in order to provide more information for understanding the advantages of the VIC model as well as more details of ETa in the PRB.

5. Conclusions

This study established a VIC model with a high spatial resolution of 0.05° in the PRB, and conducted a simulation of the ETa during the 1960–2014 period. The spatiotemporal variation of the simulated ETa and its influential climatic factors in recent decades were analyzed. The following conclusions can be drawn from the main results of this study:

(1)

The analysis of the simulated ETa derived from the VIC model showed that the mean annual ETa of the PRB was 626.1 mm. Additionally, overall, the annual ETa showed a slight but not significant increasing trend in the PRB over the past 55 years, whereas it showed a negative trend during the 1960–1992 period. At the spatial scale, the ETa in the middle and upper stream of the PRB generally exhibited a non-significant decreasing trend, except for some parts which showed a relatively significant decreasing trend. In the downstream areas, especially in the Pearl River Delta and Dongjiang River Basin, the ETa exhibited a significant increasing trend.

(2)

According to the analysis of the partial correlation coefficients, the partial correlation coefficients between the ETa and climatic factors varied at different time scales. From 1960 to 1992, the ETa in the PRB had significant correlations with TEMP and VP at the yearly scale. Additionally, it also significantly correlated with SH and PRESS at the monthly scale. From 1992 to 2014, it had a significant correlation with TEMP, SH, and PRESS. In general, after considering the variations of these factors, the changes of the ETa across the PRB were mainly influenced by SH and PRESS.

(3)

Overall, mainly due to the decreasing SH, the “paradox” phenomenon existed in the middle-upper area of the PRB during the period of 1960–1992. The phenomenon could be explained by the “global dimming” event in China, resulting from the changes of cloud coverage and aerosol accumulation. The “paradox” phenomenon in the PRB can be detected by ETa estimation, which implies we could also use ETa to further verify the “evaporation paradox”.