Study on the Projection of Potential Evapotranspiration Trend in Yangtze River Delta and the Uncertainty

1. Introduction

Evapotranspiration is a key physical process related to the water and energy balance in climate and hydrological systems, which plays a vital role in studying climate change, water resources development and utilization, crop water demand management, drought warning and monitoring, etc. The actual evapotranspiration is often invisible and difficult to monitor [1]. Therefore, it is usually estimated from the potential evapotranspiration (Etp) through methods based on hydrology, micrometeorology, and et al. Etp is the evaporation from an extended surface of the short green crop which fully shades the ground, exerts little or negligible resistance to the flow of water, and is always well supplied with water [2]. Etp is not only an important measure of atmospheric evapotranspiration capacity, but also an important research content of land surface water and energy balance. It is of great significance for regional farmland irrigation, drought evaluation, and watershed runoff simulation and estimation. Previous studies show that Etp is closely related to air temperature, solar radiation, and wind speed [3,4,5,6,7]. However, the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) points out that human activities caused global warming [8]. In comparison with the period 1850~1900, the global mean surface temperature increased by 0.99 (0.84~1.10) °C during the period 2001~2020 and by 1.09 (0.95~1.20) °C during the period 2011~2020 [8]. Global warming may increase the evapotranspiration capacity of the atmosphere, thereby affecting the availability of clean water resources in a region. Therefore, projecting the Etp trend under future climate change scenarios is beneficial for adaptation to the adverse effects of climate change.

Predicting the response of Etp to climate change usually takes a ‘top-down’ approach which combines the outputs of global climate models (GCMs) with Etp calculation formulas. To describe the physical processes in the ocean-atmosphere-land-biosphere-cryosphere system, GCMs are constructed according to the Navi-Stoker equation of the rotating sphere [8,9]. It simulates the global climate system according to the mass continuity equation, energy conservation equation, state equation, and other static approximate equations. Then the global climate change in the next hundred years can be projected through GCMs according to the future greenhouse gas and aerosol emission scenarios, which is constructed according to the future economic development, population growth, science and technology development, and et al. The calculation formulas of Etp usually take air temperature, solar radiation, and wind speed etc. as inputs, and calculate the corresponding Etp via mathematical models which abstract and generalize the processes in evapotranspiration, making it possible to project the Etp response to climate change. However, the projected results of different GCMs and Etp calculation formulas usually differ from one another. This means that the Etp projections provided by this framework include multi-source uncertainties, which will reduce the reliability and robustness of the Etp projections to a certain extent [10]. The estimation and decomposition of uncertainty and a further analysis on its main sources can point out the direction for reducing the uncertainty of response projection and improving the robustness of this projection. At the same time, the adverse influence of uncertainty on the response projection can be evaluated, and a reference on the scientific formulation of adaptation policy can be provided.

Yangzi River Delta (114°53’E, 27°03’N ~ 122°57’E, 35°08’N) is located in the east coast of China, covering a land area of 359 thousand km². There are three provinces (Anhui, Jiangsu, Zhejiang) and a municipality (Shanghai) directly under the central government in this delta, formulating one of the ten largest city clusters in the word. It is one of the two economic deltas with the most active economic development, the highest degree of openness, and the strongest innovation capacity in China. It plays a pivotal strategic position in the overall situation of national modernization and all round opening-up. The developed economic of Yangzi River Delta is closely related to its pleasant climate (with the average annual temperature being 16 °C and the average annual precipitation being 1200 mm) and abundant water resources (with the first and third largest rivers in China flowing through it). However, global climate change has exerted some certain influence on the water resources in Yangzi River Delta. For instance, although the precipitation has increased slightly in the past 50 years, the total water resources in Zhejiang Province are obviously lower than the annual average since the 21st century due to the increasing evapotranspiration. Therefore, this study intends to evaluate the Etp trend in Yangzi River Delta under the background of global climate change through an impact modeling chain from 3 emission scenarios to 24 GCMs and then 8 bias correction methods and finally 4 Etp calculation formulas. Then the uncertainty of the Etp projections is estimated and decomposed to investigate its main source by using multi-way analysis of variance (ANOVA). Finally, the robustness of the projections is evaluated by using the signal-to-noise ratio (SNR).

2. Materials and Methods

2.1. Outputs of GCMs

In this study, the simulation and projection data of daily maximum temperature (TX) and daily minimum temperature (TN) of 24 GCMs (Table 1) under three emission scenarios in the Coupled Model Intercomparison Project phase 6 (CMIP6) [11] archive are used. These GCMs are forced by the historical forcing during the period 1950~2014 and by three shared socioeconomic pathways (SSPs) - representative concentration pathways (RCPs) matrix emission scenarios [12] during the period 2015~2100. The three emission scenarios include a low concentration scenario (SSP1-2.6), a medium stabilization concentration scenario (SSP2-4.5), and a high concentration scenario (SSP5-8.5).

2.2. Observed Meteorological Datasets

In this study, a 0.5°×0.5° grid dataset of daily surface temperature in China (V2.0) is used as observed meteorological dataset. This dataset is provided by the China Meteorological Data Service Center & National Meteorological Information Center (https://data.cma.cn/data/cdcindex/cid/00f8a0e6c590ac15.html). It is based on the basic meteorological elements data of China's high-density stations (2472 national meteorological stations), and interpolated into a common 0.5°×0.5° grid over the Chinese mainland via a thin plate spline method. This dataset covers the period 1961~2019, and there are 159 grid points in or around Yangzi River Delta, as shown in Figure 1.

To evaluate the accuracy of empirical Etp calculation formulas used in this study, the daily Etp calculated via Penman-Monteith formula based on the daily China surface climate dataset (V3.0) is compared with that calculated via the empirical formulas. The daily China surface climate dataset (V3.0) is provided by the China Meteorological Data Service Center & National Meteorological Information Center (http://data.cma.cn/data/cdcdetail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0.html) as well. In this study, the records of temperature, atmosphere pressure, relative humidity, wind speed and sunshine duration in 1961~2010 of 71 stations (as shown in Figure 1) are used, since they are relatively complete.

2.3. Bias Correction

As the outputs of GCMs are usually contaminated by large systematic biases, bias correction methods are widely used to correct the original climate outputs of GCMs. To include the uncertainty caused by different classes of bias correction methods, four single-site bias correction methods (Daily Bias Correction (DBC) [13], Quantile Delta Mapping (QDM) [14], Multivariate Bias Correction n-dimensional probability density function transform (MBCn) [15], Two-Stage multivariate Quantile Mapping (TSQM) [16]) and four multi-site bias correction methods (DBC-EC, QDM-EC, MBCN-EC, TSQM-EC) are used in this study. In DBC-EC, QDM-EC, MBCn-EC, and TSQM-EC, a post-processing approach named by Empirical Copular function [17,18,19,20] is used in resorting the outputs of DBC, QDM, MBCn, and TSQM to restore the correlation among sites and variables.

Before correcting the biases of GCM outputs, the performance of the eight bias correction methods is evaluated. Specifically, the periods 1961~1980 and 1981~2000 are selected as the calibration and validation periods, respectively. The outputs of each GCM in the calibration and validation periods under the historical forcing are taken as the simulation and projection data, respectively. Each bias correction method is applied to correct the bias of GCM outputs. Then the performance of the bias correction methods in the calibration and validation periods can be evaluated.

2.4. Etp Calculation

As the temperature projections of GCMs are generally much more reliable [21,22,23] than those of other climate variables (such as wind speed and humidity), using formulas based on temperature to project the response of Etp to climate change is relatively reasonable [24]. Therefore, four empirical formulas based on temperature variables are used in this study to calculate the daily Etp, including the Blaney-Criddle formula [25] (Eq. (1)), Hargreaves-Samani formula [6] (Eq. (2)), Makkink formula [7] (Eq. (3)), and Priestley and Taylor formula [3] (Eq. (4)).

$$\left\{\begin{array}{l}Etp=\left\{\begin{array}{l}{K}_1\times p\times (0.46\times T+8.13)\mathit{growing}\mathit{season}\\ K_2\times p\times (0.46\times T+8.13)\mathit{non}-\mathit{growing}\mathit{season}\end{array}\right.\\ p=2\times {\omega }_s/3.14\\ {\omega }_s=\mathit{arccos}(-\mathit{tan}\phi \times \mathit{tan}\delta )\\ \delta =0.409\times \mathit{sin}(2\times \pi \times J/365-1.39)\\ T=(TX+TN)/2\end{array}\right.$$

(1)

where $Etp$ is the daily potential evapotranspiration, mm/d; $p$ is the percentage of total daytime hours for the used period, dimensionless; ${\omega }_s$ the angle of sunrise, rad; $\phi$ is the latitude, rad; $\delta$ is the magnetic declination of the sun, rad; $J$ is the day order, dimensionless; $T$ is the daily mean temperature, °C; $TX$ is the daily maximum temperature, °C; $TN$ is the daily minimum temperature, °C; $K_1$ and $K_2$ are undetermined coefficients [monthly consumptive use coefficient, depending on vegetation type, location and season and for the growing season, varying from 0.5 for orange tree to 1.2 for dense natural vegetation, respectively 0.85 and 0.45 for the growing season (April to September) and the non-growing season (October to March) by default].

$$\left\{\begin{array}{l}Etp=K_3\times R_a\times \sqrt{(TX-TN)}\times (T+K_4)/\lambda \\ R_a=24\times 60\times 0.0820\times d_r\times ({\omega }_s\times \mathit{sin}\phi \times \mathit{sin}\delta +\mathit{cos}\phi \times \mathit{cos}\delta \times \mathit{sin}{\omega }_s)/\pi \\ d_r=1+0.033\times \mathit{cos}(2\times \pi \times J/365)\\ \lambda =2.501-T\times 2.361/1000\end{array}\right.$$

(2)

where $R_a$ is the extraterrestrial radiation, MJ/(m²×d); $d_r$ is the relative distance between the earth and the sun, dimensionless; $\lambda$ is the latent heat of vaporization, MJ/kg; $K_3$ and $K_4$ are undetermined coefficients (respectively 0.0023 and 17.8 by default).

$$\left\{\begin{array}{l}Etp=K_5\times \Delta \times R_s/[(\Delta +\gamma )\times \lambda ]+K_6\\ \Delta =4098\times 0.6108\times \mathit{exp}(\frac{17.27\times T}{T+237.3})/{(T+237.3)}^2\\ R_s=KT\times R_a\times {(TX-TN)}^{0.5}\\ KT=0.00185\times {(TX-TN)}^2-0.0433\times (TX-TN)+0.4023\\ \gamma =0.00163\times 101.3\times {(\frac{293-0.0065\times z}{293})}^{5.26}/\lambda \end{array}\right.$$

(3)

where $\Delta$ is the slope of the relationship between saturation vapor pressure and temperature, kPa/°C; $R_s$ is the solar radiation, MJ/(m²×d); $KT$ is the empirical coefficient, dimensionless; $\gamma$ is the psychometric constant, kPa/°C; $z$ is the altitude, m; $K_5$ and $K_6$ are undetermined coefficients (respectively 0.61 and -0.12 by default).

$$\left\{\begin{array}{l}Etp=K_7\times \frac{\Delta }{\Delta +\gamma }\times \frac{R_n}{\lambda }{+\mathrm{K}}_8\\ R_n=7.14\times {10}^{-3}\times R_s+5.26\times {10}^{-6}\times R_s\times {(T+17.8)}^{1.87}\\ -3.94\times {10}^{-6}\times R_s^2-2.39\times {10}^{-9}\times R_s^2\times {(T-7.2)}^2-1.02\end{array}\right.$$

(4)

where $R_n$ is the net radiation, MJ/(m²×d); $K_7$ and $K_8$ are undetermined coefficients (respectively 30.1 and 0 by default).

Usually, empirical formulas using the default coefficients are employed in previous studies [10,26,27,28,29,30]. However, large errors resulted for most empirical formulas when using the original constant coefficients and local calibration of these coefficients (hereafter undetermined coefficients) could reduce these errors to a large extent [31,32,33]. Therefore, undetermined coefficients of the four empirical formulas are locally calibrated in this study. Specifically, the periods 1961~1990 and 1991~2010 are selected as the calibration and validation periods, respectively. The Etp in the calibration period calculated via the Penman-Monteith formula [4] (Eq. (5)) using the daily China surface climate dataset is selected as the actual Etp. The Penman-Monteith formula is based on energy balance and water vapor diffusion, and has been proved to have high accuracy [5]. In this study, a revised version of this formula recommended by the Food and Agriculture Organization of the United Nations in 1998 is used. The undetermined coefficients are adjusted station by station to make the Etp calculated via each empirical formula using the same dataset as close as possible to the actual Etp. Then the four empirical formulas (respectively using the default and calibrated coefficients) and the Penman-Monteith formula are applied to calculated the Etp during the calibration and validation periods by using the daily China surface climate dataset as well. The mean R² [Eq. (6)] and mean absolute value of relative error [ARE, Eq. (7)] during the calibration and validation periods across all meteorological stations are calculated for the four empirical formulas. The performance of the four empirical formulas before and after locally calibration is evaluated in terms of R², relative error, inter-annual trend, and intra-annual distribution.

$$Etp=\frac{0.408\times \Delta \times (R_n-G)+\gamma \times \frac{900\times u_2\times (e_s-e_a)}{T+273}}{\Delta +\gamma \times (1+0.34\times u_2)}$$

(5)

where $R_n$ is the net radiation, MJ/m²/d; $G$ is the soil heat flux, MJ/m²/d; $u_2$ is the wind speed at 2m height, m/s; $e_s$ is the saturation vapor pressure, kPa; $e_a$ is the actual vapor pressure, kPa.

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n{({\widehat{y}}_i-\overline{y})}^2}}{{\displaystyle \sum _{i=1}^n{(y_i-\overline{y})}^2}}$$

(6)

where $y_i$ is the ith observed value; ${\widehat{y}}_i$ is the ith simulated value; $\overline{y}$ is the mean value of all observed values; $R^2$ is the coefficient of determination.

$$\mathrm{A}RE=100\times \left|\frac{\overline{\widehat{y}}-\overline{y}}{\overline{y}}\right|\%$$

(7)

where $\overline{\widehat{y}}$ is the mean value of all simulated values; $ARE$ is the absolute value of relative error.

After the evaluation, to take full advantage of the daily China surface climate dataset, the undetermined coefficients of the four empirical formulas are locally calibrated again with the Penman-Monteith formula being the benchmark while using this dataset during the period 1961~2010, then interpolated into the meteorological grids via the thiessen polygon method, and applied when simulating the historical and projecting the future Etp.

2.5. Uncertainty Estimation and Decomposition

In this study, the multi-way ANOVA is used to quantitatively estimate and decompose the uncertainty of response projections, so as to make it possible to calculate the relative contribution of each uncertainty component to the total uncertainty [34,35,36,37]. The basic idea of ANOVA is to decompose variation. Specifically, the explanatory variables are first set to several levels. Then the orthogonal combination of these variables is carried out and the corresponding dependent variable numerical matrix can be obtained. Finally, the total variance of the dependent variable is decomposed into the variance resulting from each explanatory variable individually and the variance resulting from two or more explanatory variables interactively. Taking the uncertainty of the increase projection of TX or TN after bias correction for instance, the uncertainty estimating and decomposing processes in the three-way ANOVA are illuminated as follows (those in other multi-way ANOVAs are similar).

Firstly, the increase of TX or TN in the projecting period relative to the reference period is calculated as follows:

$${\Delta }_{s,g,b}=Inde{x}_{s,g,b}^{pro}-Inde{x}_{s,g,b}^{ref};s=1,\mathrm{...},s\text{'};g=1,\mathrm{...},g\text{'};b=1,\mathrm{...},b\text{'}$$

(8)

where the subscript $s$ is serial number of emission scenario; $s\text{'}$ is the amount of emission scenarios; $g$ is the serial number of GCM; $g\text{'}$ is the amount of GCMs; $b$ is the serial number of bias correction method; $b\text{'}$ is the amount of bias correction methods; $Inde{x}_{s,g,b}^{pro}$ and $Inde{x}_{s,g,b}^{ref}$ are respectively the annual mean TX or TN during the projection and reference periods for the sth emission scenario, gth GCM, and bth bias correction method; ${\Delta }_{s,g,b}$ is the increase of TX or TN in the projection period relative to reference period under the sth emission scenario, gth GCM, and bth bias correction method.

Then the means of ${\Delta }_{s,g,b}$ are calculated as:

$$\{\begin{array}{l}\overline{{\Delta }_{\ast ,g,b}}=\frac{1}{s^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{s,\ast ,b}}=\frac{1}{g^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{s,g,\ast }}=\frac{1}{b^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{\ast ,\ast ,b}}=\frac{1}{s^{\prime }\times g^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}}{\displaystyle \sum _{g=1}^{g^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{\ast ,g,\ast }}=\frac{1}{s^{\prime }\times b^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}}{\displaystyle \sum _{b=1}^{b^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{s,\ast ,\ast }}=\frac{1}{g^{\prime }\times b^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}}{\displaystyle \sum _{b=1}^{b^{\prime }}}{\Delta }_{s,g,b}\\ \overline{{\Delta }_{\ast ,\ast ,\ast }}=\frac{1}{s^{\prime }\times g^{\prime }\times b^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}}{\displaystyle \sum _{g=1}^{g^{\prime }}}{\displaystyle \sum _b^{b^{\prime }}}{\Delta }_{s,g,b}\end{array}$$

(9)

Finally, the uncertainty components and total uncertainty are calculated as:

where $S$ is the uncertainty resulting from the main effect of emission scenario; $G$ is the uncertainty resulting from the main effect of GCM; $B$ is the uncertainty resulting from the main effect of bias correction method; $SG$ is the uncertainty resulting from the interaction effect between emission scenario and GCM; $SB$ is the uncertainty resulting from the interaction effect between emission scenario and bias correction method; $GB$ is the uncertainty resulting from the interaction effect between GCM and bias correction method; $SGB$ is the uncertainty resulting from the interaction effect among emission scenario, GCM, and bias correction method; and $T$ is the total uncertainty.

2.6. SNR

The influence of uncertainty on the change projections of temperature or Etp is that it interferes with the projections greatly, which makes it difficult for decision-makers to eliminate the influence of change noise and further identify the change signal accurately. In this study, the SNR proposed by Hawkins and Sutton [38] is used to quantify the influence of uncertainty on change signal. In this index, the mean of all increase projections is defined as the change signal, and the 90% confidence interval of all increase projections (approximately equal to 1.65 times the square root of total uncertainty) is defined as the change noise. Then the SNR is calculated as:

$$SNR=\overline{\Delta }/(1.65\times \sqrt{T})$$

(11)

where $\overline{\Delta }$ is the mean of all increase projections；$T$ is the total uncertainty.

2.7. Studying Route

In this study:

(1) The period 1971~2000 is selected as the reference period, and the periods 2021~2050 and 2061~2090 are selected to represent the near and far projecting periods. With the observed TX and TN as reference, the eight bias correction methods are used to correct the TX and TN outputs of the 24 GCMs under the historical forcing and three emission scenarios. Then the temperature trend in Yangzi River Delta under the background of global climate change is evaluated by using the corrected TX and TN. The uncertainty of temperature projections is estimated and decomposed by using three-way ANOVA, and its influence on temperature change signal is quantified by using SNR.

(2) The daily Etp at each grid point is estimated via the four empirical Etp calculation formulas by using the corrected TX and TN as inputs. Then the Etp trend during the two projection periods in Yangzi River Delta is evaluated. The uncertainty of Etp predictions is estimated and decomposed via four-way ANOVA, and the relative contribution of each uncertainty component to the total uncertainty is quantified. The influence of uncertainty on Etp increase signal is quantified by SNR.

3. Results

3.1. Performance of GCMs

Good simulation ability is a prerequisite for model projection. Therefore, the performance of GCMs is evaluated by comparing the simulated temperature to that in the 0.5°×0.5° grid dataset of daily surface temperature in China (V2.0) during the period 1961~2014. Since the systematic error in the outputs of GCMs is non-negligible and can be corrected by bias correction methods [8,13,14,15,16,17,18,19,20,21], the anomaly with the mean value in 1971~2000 as the reference is calculated for each meteorological grid and each climate variable. Then the across-grid mean inter-annual trend and intra-annual distribution is calculated for each GCM and the observation, and shown in Figure 2. The results reveal that the simulated inter-decadal trend and intra-annual distribution are generally consistent with those in the observation for both TX and TN. Specifically, for TX and TN, both the simulated and observed inter-decadal trends display a decline tendency in 1960s, a general stable tendency from the early 1970s to middle 1980s, a fast increasing tendency from the middle 1980s to 1990s, but a slowdown increasing tendency in the early 21st century. In addition, the simulated ensemble mean intra-annual distribution of the 24 GCMs is very close to the observed intra-annual distribution for both TX and TN. Therefore, the ensemble of GCMs used in this study is capable to reproduce the inter-decadal trend and intra-annual distribution, while it is still unable to reproduce the inter-annual trend of the observed TX and TN. Therefore, the inter-decadal rather than the inter-annual changes of temperature and Etp is analyzed in this study.

3.2. Performance of Bias Correction Methods

Prior to projecting TX, TN, and Etp, the performance of bias correction methods in the calibration and validation periods is evaluated, as shown in Figure 3, using six evaluation statistics from three perspectives: minimum value (5% quantile), mean value, and maximum value (95% quantile). The six statistics are: the errors in the minimum value of TX, the errors in the mean value of TX, the errors in the maximum value of TX, the errors in the minimum value of TN, the errors in the mean value of TN, and the errors in the maximum value of TN. For convenience of analysis, the mean absolute values of the six statistics across all meteorological grids are calculated. The probability distribution of the mean absolute values of the 24 GCMs is presented as a boxplot for each metric. The results show that the raw outputs of GCMs are considerably biased, with the mean absolute values ranging between 1.34°C and 2.52°C. However, the systematic bias of climate model simulations is largely reduced by the eight bias correction methods. Specifically, after bias correction, the mean absolute values generally range from 0.00°C to 0.01°C in the calibration period and from 0.22°C to 0.91°C in the validation period. For all bias correction methods, the absolute errors in the validation period are obviously larger than those in the calibration period, which is due to the inconsistency in the bias of GCM outputs [39,40]. In addition, for the two temperature variables, the performance of bias correction methods are generally comparable. Overall, the eight bias correction methods used in this study can effectively reduce the systematic bias of GCMs outputs.

3.3. Temperature Projections

Firstly, the mean increase of TX or TN across GCMs and bias correction methods in each projection period under each emission scenario is calculated for each grid, then the mean increase in each city is calculated via the thiessen polygon method and displayed in Figure 4. The results showed that, in relative to the reference period, both TX and TN obviously increase in the two projection periods, and the increase during the period 2061~2090 is visibly greater than that during the period 2021~2050. Both TX and TN exhibit the maximum increases in the SSP5-8.5 scenario during the two projection periods, and exhibit the minimum increases in the SSP2-4.5 scenario during the period 2021~2050 (but this increase is very close to that in the SSP1-2.6 scenario) and SSP1-2.6 scenario during the period 2061~2090. In addition, the maximum difference among temperature increases during the period 2061~2090 (2.11°C) is obviously larger than that during the period 2021~2050 (0.27°C), which means the uncertainty of temperature projections may increase with time. Moreover, the increase of TN is generally lower than that of TX, except for the SSP2-4.5 scenario during the 2021~2050 period. In addition, it seems almost impossible (extremely difficult) to limit the temperature rise over the Yangzi River Delta to 1.5°C (2.0°C) at the end of the 21st century, even in the low emission scenario (SSP1-2.6). Specifically, the increases of TX/TN under the emission scenarios of SSP1-2.6, SSP2-4.5, and SSP5-8.5 in Yangzi River Delta are respectively 1.74/1.67°C, 1.63/1.66°C, and 1.92/1.91°C during the period 2021~2050, and they are respectively 2.34/2.19°C, 2.98/2.88°C, and 4.41/4.32°C during the period 2061~2090. In terms of the spatial distribution, the mean increases of the two temperature variables are usually lower in southeastern coastal regions (e.g., Zhoushan City, Wenzhou City, and Taizhou City in Zhejiang Province) and higher in northwestern inland regions (e.g., Bozhou City, Fuyang City, Huaibei, and Suzhou in Anhui Province, and Xuzhou City in Jiangsu Province). This is expected, as the temperature increasing rate in oceans is much lower than that in lands [41].

3.4. Uncertainty of Temperature Projections

In this study, the total uncertainty of the increase projections in TX and TN is estimated via the three-way ANOVA for each grid point and each projection period and displayed in Figure 5. The results show that the total uncertainty observably increases with time for the two temperature variables. Specifically, the mean total uncertainties of TX and TN during the period 2021~2050 in Zhejiang Province are respectively 0.28°C² and 0.24°C², and they rise to 1.49°C² and 1.44°C² during the period 2061~2090. However, the rapid increase of total uncertainty over time does not alter its general spatial distribution patterns. Specifically, for the two temperature variables, the total uncertainties are lower in the southeastern coastal regions (e.g., Zhoushan City, Wenzhou City, and Taizhou City in Zhejiang Province) and higher in the northwestern inland regions (e.g., Fuyang City, Luan City, and Bozhou City in Anhui Province). In addition, the spatial distribution pattern of the total uncertainty is similar to that of the temperature increase. This is mainly because a large increase usually shows a large difference among increase projections, which means a large uncertainty.

Further, the total uncertainty is decomposed via the three-way ANOVA. Then the magnitudes and relative contributions of all uncertainty components are calculated, and their mean values in Yangzi River Delta are calculated via the thiessen polygon method and shown in Table 2. The results show that, for TX and TN, the uncertainty components resulting from the main effect of emission scenario, the main effect of GCM, and the interaction effect between emission scenario and GCM all increase visibly along with the projection lead times, and the increase of the uncertainty component resulting from the main effect of emission scenario is obviously larger than that resulting from the main effect of GCM. The increase of uncertainty component resulting from the interaction effect between emission scenario and GCM model is very small. Other uncertainty components are all close to 0 during the two projection periods. This indicates that the increase of total uncertainty is mainly due to the substantial increase of the uncertainty component resulting from the main effect of emission scenarios, followed by that resulting from the main effect of GCM, and lastly by the slight increase of uncertainty component resulting from the interaction effect between emission scenarios and GCM. The total uncertainty is dominantly contributed by the main effect of GCM (about 90%~91%) during the period 2021~2050, but by the main effects of emission scenario (about 51%~55%) and GCM (about 41%~44%) during the period 2061~2090. Therefore, further developing GCMs and optimizing emission scenarios are reasonable ways to reduce the uncertainty of temperature projections.

In order to evaluate the influence of uncertainty on temperature increase projections, the SNR values of TX and TN are calculated for each grid point and each projection period, then its mean value in each city is calculated by the thiessen polygon method and displayed in Figure 6. The results show that the SNR values for the two temperature variables all decrease visibly over time, indicating that the robustness of temperature increase projections decreases along with the projection lead times. In addition, the SNR of TN is slightly higher than that of TX during the period 2021~2050, and very close to that of TX during the period 2061~2090. To be more specific, the mean SNR values of TX and TN in Yangzi River Delta are respectively 2.03 and 2.19 during the period 2021~2050, while they decrease to 1.61 and 1.58 during the period 2061~2090. Even so, the SNR values are obviously higher than 1.0 for the two temperature variables during the two projection periods, indicating that the temperature increase projections in Yangzi River Delta are robust and reliable. In terms of spatial distribution, the SNR values of the two temperature variables show similar spatial distribution patterns in the two projection periods. To be more specific, the SNRs are generally higher in the eastern coastal regions (e.g., Yancheng City, Lianyungang City, and Nantong City in Jiangsu Province, Taizhou City and Zhoushan City in Zhejiang Province) and lower in the western inland regions (Anqing City, Chizhou City, Huangshan City, Luan City, and Tongling City in Anhui Province). This spatial distribution pattern is a result from the comprehensive effects of the temperature increase and uncertainty.

3.5. Performance of the Empirical Etp Calculation Formulas

To evaluate the performance of the locally calibration of undetermined coefficients in the four empirical Etp calculation formulas, the across-station mean R² and ARE during the calibration and validation periods are calculated for the four formulas respectively using the default and calibrated coefficients and displayed in Figure 7. The results show that the locally calibration of undetermined coefficients can improve the performance of empirical Etp calculation formulas to a large extent. To be more specific, for the four empirical formulas during the calibration and validation periods, the mean R² is all lower than 0.68 and the mean ARE is all higher than 7.0% when using the default coefficients. When using the locally calibrated coefficients, however, the mean R² is all higher than 0.82 and the mean ARE is all lower than 5.0% during the two periods. In addition, the mean R² and ARE during the validation period are generally close to those during the calibration period, indicating that the performance of the four formulas is steady.

Further, to evaluate the performance of Etp calculation formulas, the across-station mean inter-annual trends and intra-annual distributions of daily Etp during the period 1961~2010 are calculated for the four empirical formulas using the locally calibrated coefficients and the Penman-Monteith formula and presented in Figure 8. The results show that the inter-annual trends and intra-annual distributions of daily Etp calculated by the four formulas are generally consistent with those by the Penman-Monteith formula. Specifically, the across-station and across-year mean Etp values calculated by the Blaney-Criddle formula, Hargreaves-Samani formula, Makkink formula, and Priestley-Taylor formula are respectively 2.59mm/d, 2.53mm/d, 2.49mm/d, and 2.56mm/d, which are very close to that by the Penman-Monteith formula (2.54mm/d). In addition, the general increasing trend after the 1990s obtained by the Penman-Monteith formula is well captured by the four empirical formulas. Moreover, all the four empirical formulas show a general higher in summer while lower in winter intra-annual distribution pattern, which is consistent with that of the Penman-Monteith formula. However, it is worth mentioning that there is an obvious decreasing trend before the 1980s obtained by the Penman-Monteith formula, which is not completely captured by the Blancy-Criddle formula.

3.6. Etp Projections

The daily Etp is calculated via the four empirical formulas according to the corrected TX and TN. Then its mean increase in each projection period relative to the reference period is calculated for each grid point, each GCM, each bias correction method, and each Etp calculation formula. Finally, the mean increase of Etp in Yangzi River Delta under each emission scenario and each projection period across the 24 GCMs, 8 bias correction methods, and 4 empirical Etp calculation formulas is calculated via the thiessen polygon method and displayed in Figure 9. The results show that, in relative to the reference period, the annual mean daily Etp shows an obvious increasing trend for each projection period and each city in Yangzi River Delta, and the increase during the period 2061~2090 is obviously higher than that during the period 2021~2050. Similar to the two temperature variables, the increase in Etp is highest in the SSP5-8.5 scenario during the two projection periods, and lowest in the SSP2-4.5 scenario (but this increase is very close to that in the SSP1-2.6 scenario) during the period 2021~2050 but SSP1-2.6 scenario during the period 2061~2090. To be more specific, the mean increase of annual mean daily Etp in Yangzi River Delta during the period 2021~2050 is respectively 0.15mm/d, 0.14mm/d, and 0.17mm/d under the SSP1-2.6, SSP2-4.5, and SSP5-8.5 scenarios. During the period 2061~2090, they rise to 0.21mm/d, 0.27mm/d, and 0.41mm/d, respectively. In addition, the increase of Etp shows similar spatial distribution patterns under the three emission scenarios during the two projection periods. Specifically, the increase of Etp is usually lower in southeastern coastal regions (Zhoushan City, Wenzhou City, and Taizhou City in Zhejiang Province) and higher in northwestern inland regions (Bozhou City, Fuyang City, and Huaibei City in Anhui Province). In addition, the spatial distribution pattern of daily Etp increase is consistent with that of the two temperature variables in general. This is mainly because Etp and temperature usually show obvious positive correlation.

3.7. Uncertainty of Etp Projections

The total uncertainty of the mean increase projections in daily Etp is calculated via the four-way ANOVA and displayed in Figure 10. The results show that the total uncertainty increases with time as well. To be more specific, the mean value of total uncertainty in Yangzi River Delta is 0.004(mm/d)² during the period 2021~2050, and rises to 0.022(mm/d)² during the period 2061~2090. The total uncertainty is generally lower in the southeastern coastal regions (e.g., Zhoushan City, Wenzhou City, and Taizhou City in Zhejiang Province) and higher in the northwestern inland regions (e.g., Bozhou City and Fuyang City in Anhui Province) in the two projection periods. Moreover, the spatial distribution pattern of the total uncertainty for Etp is consistent with that of the two temperature variables in general.

After that, the total uncertainty is decomposed via the four-way ANOVA, and the magnitude and relative contribution of each uncertainty component are calculated for each grid point and each period. Then their mean values in Yangzi River Delta are calculated via the thiessen polygon method and displayed in Table 3. The results show that the uncertainty components resulting from the main effect of emission scenario, the main effect of GCM, and the main effect of Etp calculation formula all show obvious increasing trends along with the projection lead times. The increase in uncertainty component resulting from the main effect of GCM is lower than that resulting from the main effect of emission scenario and higher than that resulting from the main effect of Etp calculation formula. This indicates that the increase in the total uncertainty is mainly due to the increase in the uncertainty component resulting from the main effect of emission scenario, and then by that resulting from the main effects of GCM and Etp calculation formula. In addition, the total uncertainty is dominantly contributed by the uncertainty component resulting from the main effect of GCM (63%), followed by that resulting from the main effect of Etp calculation formula (24%) during the period 2021~2050. During the period 2061~2090, it is mainly contributed by the uncertainty component resulting from the main effect of GCM (36%) as well, followed by that resulting from the main effects of emission scenario (34%) and Etp calculation formula (18%). The magnitude and relative contribution of uncertainty component resulting from the main effect of bias correction method are almost negligible. Among the uncertainty components resulting from interaction effects, those resulting from the interaction effect between emission scenario and GCM and the interaction effect between GCM and Etp calculation formula are considerable, while those resulting from other interaction effects are quite small. Therefore, further developing GCMs and Etp calculation formulas and optimizing emission scenarios are reasonable ways to reduce the uncertainty of Etp projections.

In order to evaluate the influence of uncertainty on the Etp increase projections, the SNR is calculated for each grid point and each projection period, and then its mean value in each city is calculated via the thiessen polygon method and displayed in Figure 11. The results show that SNR during the period 2021~2050 is obviously higher than that during the period 2061~2090, indicating that the robustness of Etp projections decreases along with the projection lead times. To be more specific, the mean SNR in Yangzi River Delta is 1.48 during the period 2021~2050, and decreases to 1.22 during the period 2061~2090. However, the SNR values in the two projection periods are all higher than 1.0, indicating that the projections of Etp increase in Yangzi River Delta are relatively robust and reliable. In terms of spatial distribution, SNR is higher in the southeastern coastal regions (e.g., Taizhou City, Wenzhou City, Lishui City, and Ningbo City in Zhejiang Province) and lower in the northwestern inland regions (e.g., Fuyang City, Bozhou City, Huangshan City, and Chizhou City in Anhui Province) in the two projection periods.

3.8. Discussion

In this study, the Etp trend in Yangzi River Delta is projected by using an impact modeling chain from 3 emission scenarios to 24 GCMs and then 8 bias correction methods and finally 4 Etp calculation formulas. The results reveal that, in comparison with the reference period (1971~2000), the mean Etp will increase by 0.14~0.17mm/d during the period 2021~2050 and by 0.21~0.41mm/d during the period 2061~2090, respectively. The overall framework and conclusion are similar to those in most previous studies [10,24,42]. For instance, in the study by Shi, et al. [10], the Etp trend of eight meteorological stations in southeastern Australia is projected by using 2 emission scenarios (RCP4.5 and RCP8.5, which are respectively close to the SSP2-4.5 and SSP5-8.5 scenarios), 34 GCMs, and 7 Etp calculation models (four empirical formulas and three models based on random forest). And the results reveal that, the across-station mean Etp increases range from 0.09mm/d (2040s) to 0.35mm/d (2090s) under the RCP4.5 scenario and from 0.12mm/d (2040s) to 0.68mm/y (2090s) under the RCP8.5 scenarios, respectively. Therefore, the main conclusion that Etp will increase in the 21st century can be termed robust.

In addition, the uncertainty of Etp projections is estimated and decomposed via a framework based on ANOVA, which are widely used in previous studies [34,35,36,43,44,45], as well. However, many studies exhibit inconsistent conclusions on the main sources of uncertainty. For instance, this study reveals that the total uncertainty is mainly contributed by the main effects of GCM (63%) and Etp calculation formula (24%) during the period 2021~2050, and by the main effects of GCM (36%), emission scenario (34%), and Etp calculation formulas (18%) during the period 2061~2090. However, the study by Shi, et al. [10] shows that RCP-related uncertainty contributes the most to projected Etp uncertainty (around 40% for most stations) while GCM-related and Etp model-related uncertainties account for roughly equal amounts of projected Etp uncertainty (10%~30%). For the projection of future global Etp with various Etp models and GCMs, Kingston, et al. [46] claims that Etp model-related uncertainty is equal to or, in some circumstance, greater than GCM-related uncertainty. Therefore, much more efforts are still required to confirm the main sources of uncertainty in Etp projections.

4. Conclusion

In this study, the systematic biases in temperature outputs of 24 GCMs in CMIP6 under three emission scenarios are corrected via 8 bias correction methods, and then the Etp is projected via 4 empirical Etp calculation formulas by using the corrected temperature projections. After that, the temperature and Etp trends in Yangzi River Delta in the 21st century are illuminated. Subsequently, the uncertainty in the temperature and Etp increase projections is estimated and decomposed via the multi-way ANOVA frameworks, and then its main sources are explained. Finally, the influence of uncertainty on the change signal is quantified by SNR. The main conclusions are as follows:

(1) In comparison with the period 1971~2000, the increases in the annual mean TX in Yangzi River Delta are respectively 1.63~1.92°C and 2.34~4.41°C during the periods 2021~2050 and 2061~2090, and that for TN are respectively 1.66~1.91°C and 2.19~4.32°C. The total uncertainty of TX and TN increase projections is dominantly contributed by the main effect of GCM (about 90%~91%) during the period 2021~2050, and mainly by the main effect of emission scenario (about 51%~55%) and followed by the main effect of GCM (about 41%~44%) during the period 2061~2090. The increase projections of the two temperature variables in Yangzi River Delta is robust and reliable for the two projection periods.

(2) In comparison with the period 1971~2000, the annual mean daily Etp of Yangzi River Delta will increase by 0.14~0.17mm/d during the period 2021~2050 and by 0.21~0.41mm/d during the period 2061~2090, respectively. The total uncertainty of Etp increase projections is dominantly contributed by the main effect of GCM (63%) and followed by the main effect of Etp calculation formula (24%) during the period 2021~2050, and mainly contributed by the main effect of GCM (36%) and then by the main effects of emission scenario (34%) and Etp calculation formula (18%) during the period 2061~2090. The Etp increase projections are relatively robust and reliable in the two projection periods.

(3) The robustness of response projections usually decreases with the extension of the impact modeling chain (the SNR of Etp projections is obviously lower than those of TX and TN projections). Therefore, appropriate attention may be paid to the length of the impact modeling chain when making similar response projections to climate change.

Overall, when formulating the climate change adaptation policy in Yangzi River Delta, the influence of the rising Etp should be taken into consideration. In addition, the scientific community may further develop the GCMs and Etp calculation formulas and optimize emission scenarios to reduce the uncertainty of Etp projections.