Trend Projections of Potential Evapotranspiration in Yangtze River Delta and the Uncertainty

Abstract

Global warming may increase potential evapotranspiration (ETp), reducing the water resources in Yangzi River Delta. Therefore, it is important to investigate the trend of ETp there under the background of climate change. To this purpose, the systematic biases in temperature outputs of 24 global climate models (GCMs) under 3 shared socioeconomic pathways—representative concentration pathways (SSPs) emission scenarios (SSP1-2.6, SSP2-4.5, and SSP5-8.5)—are first corrected by using 8 bias correction methods. Then, the trend of ETp in Yangtze River Delta is projected by using 4 ETp calculation formulas (Blaney–Criddle, Hargreaves–Samani, Makkink, and Priestley–Taylor). The uncertainty of the projections is estimated and decomposed by using multi-way analysis of variance frameworks. The influence of uncertainty on the projected change signal is quantified by using the signal-to-noise ratio. The results show that all emission scenarios indicate robust increments of ETp. Specifically, relative to 1971~2000, ETp will increase by 0.14~0.17 mm d⁻¹ (5.7~6.8%) during 2021~2050 and by 0.21~0.41 mm d⁻¹ (8.5~16.7%) during 2061~2090, respectively. During 2021~2050, the uncertainty of ETp projections is dominantly contributed by the main effects of GCM (63%) and the ETp calculation formula (24%). During 2061~2090, it is mainly contributed by the main effect of GCM (36%), followed by the main effects of the emission scenario (34%) and the ETp calculation formula (18%). The ETp projections are generally reliable and robust during the two projection periods.

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 resource development and utilization, crop water demand management, drought warning and monitoring, etc. [1]. The actual evapotranspiration is often invisible and difficult to monitor [2]. Therefore, it is usually estimated from potential evapotranspiration (ETp) through methods based on hydrology, micrometeorology, etc. [1]. ETp is the quantity of water vapor which could be emitted by a surface of pure water, per unit surface area and unit time, under the existing atmospheric conditions [3,4]. 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 [5,6,7,8,9] show that ETp is closely related to air temperature, solar radiation, and wind speed. However, the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) points out that human activities caused global warming [10]. 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 [10]. 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 (referred to as the ‘impact modeling chain’ as well) 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 Navier–Stokes equation of the rotating sphere [10,11]. It simulates the global climate system according to the mass continuity equation, the energy conservation equation, the 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 future economic development, population growth, science and technology development, etc. The calculation formulas of ETp usually take air temperature, solar radiation, and wind speed 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. For instance, Song et al. [12] investigated the trend of ETp in South Korea by using 2 emission scenarios, 11 GCMs and 3 ETp calculation formulas, and found that the Penman–Monteith formula showed the greatest increase and the opposite was true for the Hargreaves–Samani formula. Kingston et al. [13] projected the ETp increase to a 2 °C rise in global mean temperature by using five GCMs and six ETp calculation formulas, and found that all formulas suggest an increase in ETp but the difference in the increase signal between formulas is over 100%. 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 [14]. 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. For instance, Shi et al. [14] projected the ETp tendency in southeastern Australia by using an impact modeling chain from two emission scenarios to 34 GCMs and finally 7 ETp calculation models, and found thatETp will increase by 2.1~2.8% during the 2040s and by 9.2~17.6% during the 2090s. In addition, the uncertainty of the ETp projections is estimated and decomposed via an analysis of variance (ANOVA) framework, and the result shows that the uncertainty is mainly contributed by the emission scenario. Lai et al. [15] projected the reference crop evapotranspiration tendency in the Chinese mainland by using an impact modeling chain of four emission scenarios, five GCMs, four statistical downscaling approaches and five variable selection schemes, and found its increment during 2041~2070 is approximately 0.25~0.27 mm d⁻¹ relative to 1976~2005. A four-way ANOVA framework is applied to investigate the uncertainty of the projections, and the result shows that GCM is the main uncertainty source.

Yangtze 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 world. 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 economics of Yangtze 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 Yangtze River Delta. For instance, although precipitation has increased slightly in the past 50 years, the total water resources in Zhejiang Province are obviously lower than the intra-annual average since the 21st century, which is mainly due to increasing evapotranspiration. Therefore, this study intends to evaluate the ETp trend in Yangtze 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) frameworks. Finally, the robustness of the projections is evaluated by using the signal-to-noise ratio (SNR).

This paper is organized as follows: Section 2 describes the materials and methods used in this study; Section 3 investigates the trend of ETp and further the robustness and uncertainty of ETp projections; Section 4 provides the overall conclusions.

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. Information of the 24 GCMs in the CMIP6 archive used in this study.

Number	Name	Horizontal Resolution	Organization/Country (Region)
		Longitude × Latitude
1	ACCESS-CM2	1.8750° × 1.25°	CSIRO-ARCCSS/Australia
2	ACCESS-ESM1-5	1.8750° × 1.25°	CSIRO/Australia
3	BCC-CSM2-MR	1.125° × 1.1213°	BCC/China
4	CanESM5	2.8125° × 2.7893°	CCCma/Canada
5	CMCC-CM2-SR5	1.25° × 0.9424°	CMCC/Italy
6	CNRM-CM6-1	1.4063° × 1.4004°	CNRM-CERFACS/France
7	CNRM-ESM2-1	1.4063° × 1.4004°	CNRM-CERFACS/France
8	EC-Earth3	0.7031° × 0.7017°	EC-Earth-Consortium/European Union
9	EC-Earth3-Veg	0.7031° × 0.7017°	EC-Earth-Consortium/European Union
10	FGOALS-g3	2° × 2.2785°	CAS/China
11	GFDL-ESM4	1.25° × 1°	NOAA-GFDL/America
12	HadGEM3-GC31-LL	1.8750° × 1.25°	MOHC/England
13	INM-CM4-8	2° × 1.5°	INM/Russia
14	INM-CM5-0	2° × 1.5°	INM/Russia
15	IPSL-CM6A-LR	2.5° × 1.2676°	IPSL/France
16	MIROC6	1.4063° × 1.4004°	MIROC/Japan
17	MIROC-ES2L	2.8125° × 2.7893°	MIROC/Japan
18	MPI-ESM1-2-HR	0.9375° × 0.9349°	MPI-M/Germany
19	MPI-ESM1-2-LR	1.875° × 1.8647°	MPI-M/Germany
20	MRI-ESM2-0	1.1250° × 1.1213°	MRI/Japan
21	NESM3	1.875° × 1.8647°	NUIST/China
22	NorESM2-LM	2.5° × 1.8947°	NCC/Norway
23	NorESM2-MM	1.25° × 0.9424°	NCC/Norway
24	UKESM1-0-LL	1.875° × 1.25°	MOHC/England

) under 3 emission scenarios in the Coupled Model Intercomparison Project phase 6 (CMIP6) [16] 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 [17] 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). The selection of the GCMs and emission scenarios is mainly due to their wide applications in previous studies [18,19,20,21]. It is worth noting that the selection of SSP5-8.5 is intended to elaborate the worst rather than the business as usual case [22]. However, this case is unlikely to occur [23], as the emission of fossil carbon dioxide is overestimated and increasingly implausible in this scenario [24,25].

2.2. Observed Meteorological Datasets

In this study, a grid (0.5° × 0.5°) dataset of daily surface temperature in China (V2.0) is used as the observed meteorological dataset. This dataset is provided by the China Meteorological Data Service Center & National Meteorological Information Center. 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 Yangtze 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 the 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 as well. In this study, the records of temperature, atmosphere pressure, relative humidity, wind speed, and sunshine duration during the period 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 [26,27], bias correction methods are widely used to post-process the original outputs of GCMs. To include the uncertainty caused by different classes of bias correction methods, four single-site methods and four multi-site methods are used in this study. The single-site methods include two univariate methods {daily bias correction (DBC) [28] and Quantile Delta Mapping (QDM) [29]} and two multivariate methods {multivariate bias correction n-dimensional probability density function transform (MBCn) [30] and two-stage multivariate quantile mapping (TSQM) [31]}. In the four multi-site methods (hereafter DBC-EC, QDM-EC, MBCn-EC, and TSQM-EC), a post-processing approach named by the empirical copular (EC) function [32,33,34,35] is used in resorting the outputs of DBC, QDM, MBCn, and TSQM to restore the correlation among sites and variables. The selection of these methods is mainly due to their good performance, which will be demonstrated in the later section. The overall routes of the four single-site method are as follows:

DBC is a hybrid method combining the Daily Translation method (DT) [36] and LOCal Intensity scaling method (LOCI) [37]. LOCI is first used to correct the precipitation occurrence, ensuring that the frequency of the precipitation occurrence of the corrected data is equal to that of the observed data at the reference period for each month. DT is then applied to adjust the frequency distribution of precipitation amounts and temperature for each month. This hybrid method takes advantage of the LOCI method to correct precipitation occurrence. Then, instead of applying the same factor to each daily precipitation for a specific month, it modulates the change as a function of the daily precipitation amount (or temperature) based on differences in the frequency distribution of precipitation amount (or temperature) for the observed data.

QDM is also a quantile mapping bias correction method, while the climate change signal simulated by climate models is preserved. QDM combines two steps in preserving the climate change signal: (1) future model outputs are detrended in terms of quantile and bias corrected to observations by quantile mapping; and (2) model-projected relative changes in quantiles are superimposed on the bias-corrected model outputs.

MBCn is a multivariate generalization of quantile mapping that transfers all aspects of an observed continuous multivariate distribution to the corresponding multivariate distribution of variables from a climate model. Specifically, QDM is first utilized to correct the bias of each climate variable. The inter-variable correlations are then reconstructed by repeatedly adjusting the random linear combinations of climate variables. When applied to future climate model projections, changes in quantiles of each variable between the historical and projecting periods are also preserved.

TSQM is an extension of the DBC method with the introduction of multivariate correlations by using a distribution-free shuffle approach. Specifically, DBC is first used to correct the marginal distribution of a single variable, and then a distribution-free shuffle approach is used to introduce the target multivariable correlations.

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 simulation and projection data, respectively. Each bias correction method is applied to correct the biases 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 illuminated in the Penman–Monteith formula [6], ETp is closely relevant to temperature, radiation, atmosphere pressure, wind speed, humidity et al. [38]. However, the temperature projections of GCMs are generally much more reliable than those of other climate variables (such as wind speed and humidity) [27,39,40]. Therefore, using formulas based on temperature to project the response of ETp to climate change is relatively reasonable [14,41]. In this study, four empirical formulas based on temperature variables are used to calculate the daily ETp, the Blaney–Criddle formula [42] (Equation (1)), the Hargreaves–Samani formula [8] (Equation (2)), the Makkink formula [9] (Equation (3)), and the Priestley–Taylor formula [5] (Equation (4)). They are optimized from approximately a dozen formulas suggested by Xu and Singh [43].

$$\{\begin{array}{l}\mathrm{ETp}=\{\begin{array}{l}{K}_1\times p\times (0.46\times T+8.13)\hspace{1em}\hspace{1em}\hspace{1em}growingseason\\ K_2\times p\times (0.46\times T+8.13)\hspace{1em}non-growingseason\end{array}\\ 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}$$

(1)

where $\mathrm{ETp}$ is the daily potential evapotranspiration, mm/d; $p$ is the percentage of total daytime hours for the used period, dimensionless; ${\omega }_s$ is 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; 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.

$$\{\begin{array}{l}\mathrm{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}$$

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

$$\{\begin{array}{l}\mathrm{ETp}=K_5\times \mathsf{\Delta }\times R_s/[(\Delta +\gamma )\times \lambda ]+K_6\\ \mathsf{\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}$$

(3)

where $\mathsf{\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).

$$\{\begin{array}{l}\mathrm{ETp}=K_7\times \frac{\mathsf{\Delta }}{\mathsf{\Delta }+\gamma }\times \frac{R_n}{\lambda }+K_8\\ \begin{array}{cc}\hfill R_n=& 7.14\times {10}^{-3}\times R_s+5.26\times {10}^{-6}\times R_s\times {(T+17.8)}^{1.87}\hfill \\ & -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\hfill \end{array}\end{array}$$

(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 [14,44,45,46,47,48]. 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 [43,49,50]. 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 [6] (Equation (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 [7]. 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 calculate the ETp during the calibration and validation periods by using the daily China surface climate dataset as well. The mean R² [Equation (6)] and mean absolute value of relative error [ARE, Equation (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 local calibration is evaluated in terms of R², relative error, inter-annual trend, and intra-annual distribution.

$$\mathrm{ETp}=\frac{0.408\times \mathsf{\Delta }\times (R_n-G)+\gamma \times \frac{900\times u_2\times (e_s-e_a)}{T+273}}{\mathsf{\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 2 m 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, the coefficients are 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 frameworks are 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 [51,52,53,54]. 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 (main effect) and the variance resulting from two or more explanatory variables interactively (interaction effect). 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 ANOVA frameworks are similar).

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

$${\mathsf{\Delta }}_{s,g,b}=Inde{x}_{s,g,b}^{pro}-Inde{x}_{s,g,b}^{ref};s=1,\dots ,s^{\prime };g=1,\dots ,g^{\prime };b=1,\dots ,b^{\prime }$$

(8)

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

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

$$\{\begin{array}{l}\overline{{\mathsf{\Delta }}_{*,g,b}}=\frac{1}{s^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{\mathsf{\Delta }}_{s,g,b}}\\ \overline{{\mathsf{\Delta }}_{s,*,b}}=\frac{1}{g^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{\mathsf{\Delta }}_{s,g,b}}\\ \overline{{\mathsf{\Delta }}_{s,g,*}}=\frac{1}{b^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{\mathsf{\Delta }}_{s,g,b}}\\ \overline{{\mathsf{\Delta }}_{*,*,b}}=\frac{1}{s^{\prime }\times g^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{\mathsf{\Delta }}_{s,g,b}}}\\ \overline{{\mathsf{\Delta }}_{*,g,*}}=\frac{1}{s^{\prime }\times b^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{\mathsf{\Delta }}_{s,g,b}}}\\ \overline{{\mathsf{\Delta }}_{s,*,*}}=\frac{1}{g^{\prime }\times b^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{\mathsf{\Delta }}_{s,g,b}}}\\ \overline{{\mathsf{\Delta }}_{*,*,*}}=\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 }}{\mathsf{\Delta }}_{s,g,b}}}}\end{array}$$

(9)

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

$$\{\begin{array}{l}S=\frac{1}{s^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{(\overline{{\mathsf{\Delta }}_{s,*,*}}-\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}\\ G=\frac{1}{g^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{(\overline{{\mathsf{\Delta }}_{*,g,*}}-\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}\\ B=\frac{1}{b^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{(\overline{{\mathsf{\Delta }}_{*,*,b}}-\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}\\ SG=\frac{1}{s^{\prime }\times g^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{(\overline{{\mathsf{\Delta }}_{s,g,*}}-\overline{{\mathsf{\Delta }}_{s,*,*}}-\overline{{\mathsf{\Delta }}_{*,g,*}}+\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}}\\ SB=\frac{1}{s^{\prime }\times b^{\prime }}{\displaystyle \sum _{s=1}^{s^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{(\overline{{\mathsf{\Delta }}_{s,*,b}}-\overline{{\mathsf{\Delta }}_{s,*,*}}-\overline{{\mathsf{\Delta }}_{*,b,*}}+\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}}\\ GB=\frac{1}{g^{\prime }\times b^{\prime }}{\displaystyle \sum _{g=1}^{g^{\prime }}{\displaystyle \sum _{b=1}^{b^{\prime }}{(\overline{{\mathsf{\Delta }}_{*,g,b}}-\overline{{\mathsf{\Delta }}_{*,g,*}}-\overline{{\mathsf{\Delta }}_{*,*,b}}+\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}}\\ SGB=\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=1}^{b^{\prime }}\begin{array}{l}({\mathsf{\Delta }}_{s,g,b}-\overline{{\mathsf{\Delta }}_{s,g,*}}-\overline{{\mathsf{\Delta }}_{s,*,b}}-\overline{{\mathsf{\Delta }}_{*,g,b}}\\ +\overline{{\mathsf{\Delta }}_{s,*,*}}+\overline{{\mathsf{\Delta }}_{*,g,*}}+\overline{{\mathsf{\Delta }}_{*,*,b}}-\overline{{\mathsf{\Delta }}_{*,*,*}}{)}^2\end{array}}}}\\ T=\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=1}^{b^{\prime }}{({\mathsf{\Delta }}_{s,g,b}-\overline{{\mathsf{\Delta }}_{*,*,*}})}^2}}}\end{array}$$

(10)

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 scenarios and GCM; $SB$ is the uncertainty resulting from the interaction effect between emission scenarios 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 scenarios, GCM, and bias correction method; and $T$ is the total uncertainty.

It is worth noting that the total uncertainty is equal to the biased estimated variance of ${\mathsf{\Delta }}_{s,g,b}$. Therefore, the uncertainty components can be termed as the variance caused by the main effects or interaction effects. In addition, the fraction of each uncertainty component to the overall uncertainty is termed relative contribution (in %).

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 [55] 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:

$$\mathrm{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 projection periods. With the observed TX and TN as reference, the 8 bias correction methods are used to correct the TX and TN outputs of the 24 GCMs under the historical forcing and 3 emission scenarios. Then, the temperature trend in Yangtze 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 Yangtze River Delta is evaluated. The uncertainty of ETp projections 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

As all projections are derived from GCMs and good simulation ability is a prerequisite for model projection, the performance of GCMs is evaluated by comparing the simulated temperature to that in the grid dataset of daily surface temperature in China (V2.0) during the period 1961~2014. Two aspects, the inter-annual or inter-decadal trend and the intra-annual distribution, are taken into account in the comparison, as they are of great importance for projecting the overall trend of ETp and guaranteeing the accuracy of the projections. Since the systematic error in the outputs of GCMs is non-negligible and can be corrected by bias correction methods [10,27,28,29,30,31,32,33,34,35], 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 are calculated for each GCM and the observation, and shown in Figure 2. In addition, the R² between the observed series and the ensemble mean simulated series (for the inter-annual trend, the 10-year moving average of the two series is used) is calculated and displayed in Figure 2 as well. 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. Furthermore, R² is higher than 0.7 in terms of the inter-decadal trend and higher than 0.9 in terms of the intra-annual distribution. Overall, the ensemble of GCMs used in this study is capable of reproducing 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

Removing the systematic biases of GCM outputs via bias correction methods can improve the accuracy of temperature and ETp projections to a large extent. However, 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 generally ranging between 1.34 °C and 2.52 °C. However, the systematic biases of climate model simulations are 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 biases of GCM outputs [56,57]. In addition, for the two temperature variables, the performance of bias correction methods are generally comparable. Furthermore, they are comparable with that in previous studies [31,34,35,56,57,58] as well. Overall, the eight bias correction methods used in this study can effectively reduce the systematic biases of temperature outputs of GCMs.

3.3. Temperature Projections

Temperature over the multi-decadal and regional scale is usually projected by correcting the biases of GCMs outputs [31,34,35,56], which is followed in this study. Hence, the mean increase in 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, 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 these increases are very close to those 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 in TN is generally lower than that of TX, except for the SSP2-4.5 scenario during the period 2021~2050. Also, it seems almost impossible (extremely difficult) to limit the temperature rise over the Yangtze 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). However, ambitious mitigation efforts will reduce the increment of temperature to a large extent, especially in the far future (approximately 2.11 °C). Specifically, the increases inin TX/TN under the emission scenarios of SSP1-2.6, SSP2-4.5, and SSP5-8.5 in Yangtze 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 in 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 City, and Suzhou City 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 [59].

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 framework 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 Yangtze River Delta 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 in 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 framework. Then, the magnitudes and relative contributions of all uncertainty components are calculated, and their mean values in Yangtze River Delta are calculated via the Thiessen polygon method and shown in

Table 2. The mean magnitudes of uncertainty components and their relative contributions to the total uncertainty of the increase projections of daily maximum temperature (TX) and daily minimum temperature (TN) in Yangtze River Delta.

Period	2021~2050
Effect	S	G	B	SG	SB	GB	SGB
TX/Magnitude/°C2	0.01	0.25	0.00	0.01	0.00	0.00	0.00
TN/Magnitude/°C2	0.01	0.21	0.00	0.01	0.00	0.00	0.00
TX/Relative contribution/%	4.93	90.76	0.00	4.31	0.00	0.00	0.00
TN/Relative contribution/%	5.58	89.68	0.00	4.75	0.00	0.00	0.00
Period	2061~2090
Effect	S	G	B	SG	SB	GB	SGB
TX/Magnitude/°C2	0.75	0.66	0.00	0.07	0.00	0.00	0.00
TN/Magnitude/°C2	0.79	0.60	0.00	0.05	0.00	0.00	0.00
TX/Relative contribution/%	50.79	44.31	0.00	4.90	0.00	0.00	0.00
TN/Relative contribution/%	55.06	41.30	0.00	3.64	0.00	0.00	0.00

Note: The meanings of S, G, B, SG, SB, GB, and SGB are consistent with those in Equation (10).

. The results show that, for TX and TN, the uncertainty components resulting from the main effect of emission scenarios, the main effect of GCM, and the interaction effect between emission scenarios and GCM all increase visibly along with the projection lead times. The increase in the uncertainty component resulting from the main effect of emission scenarios is obviously larger than that resulting from the main effect of GCM. The increase in uncertainty component resulting from the interaction effect between emission scenarios and the GCM model is very small. Other uncertainty components are all close to 0 during the two projection periods. This indicates that the increase in total uncertainty is mainly due to the substantial increase in 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 in uncertainty component resulting from the interaction effect between emission scenarios and GCM. The total uncertainty is dominantly contributed by the main effect of GCM (90~91%) during the period 2021~2050, but by the main effects of emission scenarios (51~55%) and GCM (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 Yangtze River Delta are, respectively, 2.03 and 2.19 during the period 2021~2050, and 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 Yangtze 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 higher than 0.82 and the mean ARE is 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, the Hargreaves–Samani formula, the Makkink formula, and the Priestley–Taylor formula are, respectively, 2.59 mm d⁻¹, 2.53 mm d⁻¹, 2.49 mm d⁻¹, and 2.56 mm d⁻¹, which are very close to that by the Penman–Monteith formula (2.54 mm 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 Blaney–Criddle formula.

3.6. ETp Projections

Projecting the trend of ETp is conducive to understanding the impact of climate change on the hydrological, ecological, and social systems as ETp is closely relevant to the changes of water resources [1,8,13,60,61,62,63]. Therefore, 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 in ETp in Yangtze 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, relative to the reference period, the mean ETp shows an obvious increasing trend for each projection period and each city in Yangtze 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 in ETp in Yangtze River Delta during the period 2021~2050 is, respectively, 0.15 mm d⁻¹ (6.2%), 0.14 mm d⁻¹ (5.7%), and 0.17 mm d⁻¹ (6.8%) under the SSP1-2.6, SSP2-4.5, and SSP5-8.5 scenarios. During the period 2061~2090, they rise to 0.21 mm d⁻¹ (8.5%), 0.27 mm d⁻¹ (10.9%), and 0.41 mm d⁻¹ (16.7%), respectively. The rapid increase in ETp indicates greater demand for evapotranspiration, which may lead to higher water requirement for vegetation, animals, and humans and further the constrictive water management. In addition, the increase in ETp shows similar spatial distribution patterns under the three emission scenarios during the two projection periods. Specifically, the increase in 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 framework 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 Yangtze River Delta is 0.004 (mm d⁻¹)² during the period 2021~2050, and 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, Fuyang City, and Huaibei 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 Yangtze River Delta are calculated via the Thiessen polygon method and displayed in

Table 3. The mean magnitudes of uncertainty components and their relative contributions to the total uncertainty of the increase projections of the mean daily ETp in Yangtze River Delta.

Type	Period	S	G	B	E	SG	SB	SE	GB
Magnitude/(mm d−1)2	2021~2050	0.0001	0.0025	0.0000	0.0010	0.0001	0.0000	0.0000	0.0000
	2061~2090	0.0071	0.0077	0.0000	0.0042	0.0011	0.0000	0.0006	0.0000
Relative contribution/%	2021~2050	3.24	62.54	0.05	24.24	2.86	0.00	0.37	0.08
	2061~2090	33.76	35.59	0.04	18.26	4.76	0.01	2.81	0.03
Type	Period	GE	BE	SGB	SGE	SBE	GBE	SGBE
Magnitude/(mm d−1)2	2021~2050	0.0002	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000
	2061~2090	0.0009	0.0000	0.0000	0.0002	0.0000	0.0000	0.0000
Relative contribution/%	2021~2050	5.93	0.05	0.02	0.51	0.00	0.09	0.02
	2061~2090	3.93	0.04	0.01	0.70	0.00	0.05	0.01

Note: The meanings of S, G, B, SG, SB, GB, and SGB are consistent with those in Equation (10). E is the uncertainty component resulting from the main effect of the ETp calculation formula. SE is the uncertainty component resulting from the interaction effect between emission scenarios and the ETp calculation formula. GE is the uncertainty component resulting from the interaction effect between GCM and the ETp calculation formula. SGE is the uncertainty component resulting from the interaction effect among emission scenarios, GCM, and the ETp calculation formula. SBE is the uncertainty component resulting from the interaction effect among emission scenarios, bias correction method, and the ETp calculation formula. GBE is the uncertainty component resulting from the interaction effect among GCM, bias correction method, and the ETp calculation formula. SGBE is the uncertainty component resulting from the interaction effect among emission scenarios, GCM, bias correction method, and the ETp calculation formula.

. The results show that the uncertainty components resulting from the main effect of emission scenarios, the main effect of GCM, and the main effect of the 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 scenarios and higher than that resulting from the main effect of the 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 scenarios, and then by those resulting from the main effects of GCM and the 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 the 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 those resulting from the main effects of emission scenarios (34%) and the 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 scenarios and GCM and the interaction effect between GCM and the 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 Yangtze 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 Yangtze 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 Yangtze 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.17 mm d⁻¹ during the period 2021~2050 and by 0.21~0.41 mm d⁻¹ during the period 2061~2090, respectively. The overall framework and conclusion are similar to those in most previous studies [14,41,64]. For instance, in the study by Shi et al. [14], the ETp trend of 8 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 (4 empirical formulas and 3 models based on random forest). And the results reveal that, the across-station mean ETp increases range from 0.09 mm d⁻¹ (2040s) to 0.35 mm d⁻¹ (2090s) under the RCP4.5 scenario and from 0.12 mm d⁻¹ (2040s) to 0.68 mm d⁻¹ (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 the ANOVA frameworks, which are widely used in previous studies [51,52,53,65,66,67] 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 the ETp calculation formula (24%) during the period 2021~2050, and by the main effects of GCM (36%), emission scenarios (34%), and ETp calculation formulas (18%) during the period 2061~2090. However, the study by Shi et al. [14] 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. [13] claims that ETp model-related uncertainty is equal to or, in some circumstances, greater than GCM-related uncertainty. Therefore, much more effort is still required to confirm the main sources of uncertainty in ETp projections.

4. Conclusions

This study intends to evaluate the ETp trend in Yangtze River Delta under the background of global climate change since it plays a vital role in studying climate change and water resources. To this end, the systematic biases in temperature outputs of 24 GCMs in CMIP6 under 3 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 Yangtze 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)

Temperature

(1.1)

In comparison with 1971~2000, the temperature increases are, respectively, 1.63~1.92 °C and 2.19~4.41 °C during 2021~2050 and 2061~2090.

(1.2)

The total uncertainty of temperature projections is dominantly contributed by the main effect of GCM (90~91%) during 2021~2050, and mainly by the main effect of emission scenarios (51~55%) and followed by the main effect of GCM (41~44%) during 2061~2090.

(1.3)

The temperature projections are robust and reliable for the two projection periods.

(2)

ETp

(2.1)

In comparison with 1971~2000, the ETp will increase by 0.14~0.17 mm d⁻¹ during 2021~2050 and by 0.21~0.41 mm d⁻¹ during 2061~2090, respectively. This will lead to higher demand for evapotranspiration and may result in more constrictive water management.

(2.2)

The total uncertainty of ETp projections is dominantly contributed by the main effect of GCM (63%) and followed by the main effect of the ETp calculation formula (24%) during 2021~2050, and mainly contributed by the main effect of GCM (36%) and then by the main effects of emission scenarios (34%) and the ETp calculation formula (18%) during 2061~2090.

(2.3)

The ETp 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 (e.g., the SNR of ETp projections is obviously lower than that of temperature 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 Yangtze 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.