Ensemble Assessment of Extreme Precipitation Risk under 1.5 and 2.0°C Warming Targets in the Yangtze River Basin

1.   Introduction

The global mean surface temperature has increased since the late 19th century, with a warming of 1.2°C (IPCC, 2021). In response to changes in global temperature, intensified hydrological cycles and precipitation extremes have been observed in recent decades (Trenberth, 2011; Lehmann et al., 2015; Donat et al., 2016). Faced with threats of ongoing climate change, the “Paris Agreement 2015” set an ambitious goal to “keep global warming well below 2.0°C above preindustrial levels and pursue efforts to limit it to 1.5°C” (UNFCCC, 2015). However, until now, there has been no sign of global warming slowing down. According to the World Meteorological Organization (WMO) report (World Meteorological Organization, 2020), 2020 was the second warmest year on record, and the global mean temperature reached 1.2 ± 0.1°C above the preindustrial average (1850–1900). As the 1.5°C warming target is approaching in the near future, it is imperative to investigate changes in climate-related impacts, as described in the “Paris Agreement,” and take action, as described in COP28-2023.

As global warming continues, the capacity of the atmosphere to hold water vapor is also enhanced, and consequently, more frequent and intense extreme precipitation events may occur, which may trigger disastrous floods, particularly in regions with monsoon climates (Fischer and Knutti, 2015; Guo et al., 2016; Zhang et al., 2018; Ye and Qian, 2021; Qian et al., 2022). In recent years, extreme precipitation events, such as the great storms on 20 July 2021 in Zhengzhou and on 31 July 2023 in Beijing, have occurred increasingly frequently in North China. In the Yangtze River basin (YRB), several disastrous floods have occurred in recent decades, such as the severe flood events in 1954, 1998, and 2020, which caused serious damages to society and economy. The 1998 flood disaster caused direct economic losses of $12,815 million (Zong and Chen, 2000). In 2020, a precipitation event with record-broken intensity in the plum rain season occurred in the Yangtze–Huai River valley and upstream of the YRB, which caused severe flooding with $11.75 billion in economic loss (Ding et al., 2021; Qian et al., 2022). Given unaffordable economic losses and great stresses on both societies and natural ecosystems imposed by these extreme events, understanding extreme precipitation responses under warming scenarios is critical for flood risk assessment and disaster mitigation.

The Coupled Model Intercomparison Project (CMIP) is designed to facilitate the availability and comparability of global climate model (GCM) projections. Although models from the fifth phase of the CMIP (CMIP5) are generally capable of simulating precipitation extremes (Sillmann et al., 2013; Jiang et al., 2015), the projections still underestimate or overestimate the intensities and trends of extreme events compared to climate station-based measurements (Janssen et al., 2014; Freychet et al., 2015; Wu S. Y. et al., 2019). Its successive version, CMIP6, shows some improvements in the simulation of precipitation extremes (Yang et al., 2021, 2023; Zhu et al., 2022), which may be partly attributed to better parameterization schemes of physical processes and the appended indirect effects of aerosols on the formation of clouds and precipitation in the models (Voldoire et al., 2019; Wu T. W. et al., 2019; Xu et al., 2022). However, some studies have reported that precipitation indices are still generally overestimated (e.g., Zhu et al., 2020).

Bias correction and ensemble are so far two ways to raise the accuracies of GCMs. Theoretically, the best bias correction is via modifying physical processes to reach a better match between the simulation results from the GCM and observations. Limited by insufficiency to the processes understanding, this correction with improved physical dynamical processes has been very difficult. Another bias correction is via statistics by constructing the statistical relationship between model simulation and observation, which is usually called transfer function. There are at least four statistical bias-correction methods, i.e., the linear scaling (LS), quantile mapping (QM), distribution mapping (DM), and cumulative distribution function transform (CDFt).

QM needs to assume that the future distribution of a variable of interest will remain similar to that in the reference period. Given this assumption may not hold in realistic cases, DM methods explicitly consider the change of distribution in the future, including Equi-Distant Cumulative Distribution Functions matching method (EDCDFm) (Li et al., 2010) and Equi-Ratio Cumulative Distribution Functions matching method (ERCDFm) (Wang and Chen, 2014). EDCDFm incorporates the change in distribution with a fundamental assumption that the differences between modeled and observed values over the reference period (namely “distant”) will be preserved in a future period. Wang and Chen (2014) found that this equidistant may result in negative precipitation and established ERCDFm to preserve a consistent “ratio” of the observed data to the simulated data quantiles in the reference and projection periods. For DM method, a distribution function is needed. Gamma function is one of the popular functions used. Zhu et al. (2022) compared the performances of these four statistical bias-correction methods in correcting the simulated daily precipitation in summer (June, July, and August) covering 1951–2012 in Yangtze–Huai River basin in China from BCC-CSM1.1-m model. They found that the DM method, based on the gamma distribution, is less effective in characterizing the probability distribution of model-simulated precipitation, especially for the heavy precipitation. It needs more work to check if a gamma distribution is efficient. Another way to raise the accuracy is to mitigate the uncertainty stemming from model initialization, dynamic frameworks, and parameters. CMIP models vary noticeably in their ability to predict precipitation extremes (Wehner, 2013; Sun et al., 2022). The multimodel ensemble (MME) method, linearly combining a number of competing GCMs, has been widely applied to reduce uncertainties from GCMs. The arithmetic mean (AM) and Bayesian model averaging (BMA) methods are two of the most popular MME algorithms. AM is one of the simplest ensemble methods, in which the absolute error of the deterministic simulation increases with the ensemble range, and inconsistency among the models cancels out the climatic fluctuations (Raftery et al., 2005). An improved AM method was reported for selecting or weighting models according to the quantification of their historical performance based on metrics such as relative error, root-mean-square error (RMSE), and correlation coefficient (R²), as well as the ranked probability skill score (RPSS) (Parrish et al., 2012; Yin et al., 2016; Yang et al., 2023). However, the selection of metrics is usually subjective, and the capability of a specific model cannot be comprehensively assessed. In contrast, the BMA method overcomes this limitation by weighting a model directly based on its likelihood of predicting performance. It was first introduced to conduct meteorological model ensembles by Raftery et al. (2005) and then widely employed for the assessment of climate extremes in GCM projections (Yang et al., 2011; Demirel and Moradkhani, 2015; Basher et al., 2020; Zhao et al., 2020). The BMA method assumes a prior distribution to each model representing its uncertainty; then, the prior parameters are updated based on the simulation accuracy by using optimization techniques such as the expectation maximization (EM) algorithm, L-moments, and Markov chain Monte Carlo (MCMC) sampling (Duan et al., 2007). Thereafter, the posterior probability is derived to assign a model weight that reflects the performance of an individual model over a training period.

Given many models submitted to the CMIP are either from the same institution or share part of codes with each other, the independence of the model projection has been considered as a source of uncertainty in multimodel ensemble (Knutti et al., 2017; Sanderson et al., 2017). For a group of model projections to be assessed, the model skill, namely the distance between the model projection and observation, is not a unique index to weigh the model as in BMA. The distances between the model projections are also counted. Originally based on the concepts outlined by Sanderson et al. (2015) using the rules that the models being maximal independent and skillful in reproducing past climate, Knutti et al. (2017) then adapted this approach for constraining a specific future change of sea ice area. Finally, Sanderson et al. (2017) produced a single set of model weights, which was used to combine projections for a range of quantities into a weighted mean result, with significance estimates which also treat the weighting appropriately in conterminous United States and Canada based on CMIP5. Li T. et al. (2021) applied the method in assessing the performance and independence of GCMs in the projection of precipitation under the 1.5 and 2.0°C warming targets over China, also based on CMIP5. Zhao et al. (2022) used this scheme to analyze the change in precipitation over the Tibetan Plateau projected by weighted CMIP6. Zhao et al. (2022a, b) used this scheme in constraining CMIP6 projections of an ice-free Arctic and analyzed historical and future runoff changes in the middle and lower reaches of the YRB from CMIP6 models. The concepts are worthy of being used in testing the model replication before applying BMA for ensemble with more case studies.

In the risk analysis of climate change scenarios, there are two approaches for diagnosing extreme precipitation from climate model outputs, namely, nonparametric (index analysis) and parametric (frequency analysis) methods (Mo et al., 2019). The nonparametric method uses extreme precipitation indices to detect historical and future changes in precipitation extremes. These indices were proposed by the Expert Team on Climate Change Detection and Indices (ETCCDI) at global, national, and regional scales (Alexander et al., 2006; Chi et al., 2016; Scherrer et al., 2016; Pedron et al., 2017). Unlike the nonparametric methods, which lacks the ability to extrapolate extreme events with longer return periods (Raggad, 2018), the parametric methods use theoretical distribution functions [such as the generalized extreme value (GEV) distribution] to describe the tail of precipitation distribution and to investigate extreme events with specific return periods.

Documented studies have shown that the spatial and temporal characteristics of extreme precipitation events have been modified under global change. Much has been found in the East Asian monsoon region (Su et al., 2006; Chen et al., 2014; Gao and Xie, 2016; Zhou et al., 2018; Li X. et al., 2021; Yuan et al., 2021), where increasing extreme precipitation intensity is detected in southern China (Chen, 2013; Wu S. Y. et al., 2019). Studies have reported that the maximum 1-day precipitation (Rx1day), maximum 5-day precipitation (Rx5day), and other extreme indices of heavy and light precipitation have significantly increased in the YRB in recent decades (Gao and Xie, 2016; Yuan et al., 2021). According to CMIP3 projections under greenhouse gas emission scenarios, the increase in extreme precipitation indices will be more remarkable in the upper reaches of the YRB (Xu et al., 2009; Huang et al., 2012; Guo et al., 2013). Based on the CMIP5 projections, 20- and 50-yr return extreme precipitation in the middle–lower reaches of the YRB will increase (Pan et al., 2016; Wang et al., 2019). Studies based on CMIP6 model projections reported that precipitation extremes [total wet day precipitation (PRCPTOT), Rx5day, and very heavy precipitation (R20mm)] may increase across China under the shared socioeconomic pathway (SSP2-4.5 and SSP5-8.5; Xu et al., 2022). Specifically, the changes in precipitation extremes under the 1.5 and 2.0°C warming targets were evaluated with CMIP6 model projections (Wang et al., 2020; Guo et al., 2023), revealing that extreme precipitation in the headstream and lower reaches is more sensitive to the warming climate, while no significant change in the middle reach of the YRB is detected.

In this study, multimodel projections from bias-corrected CMIP6 projections using the BMA approach and weighting scheme are employed (1) to systematically weigh the GCMs and (2) to evaluate seasonal changes in the magnitude and risk of extreme precipitation in the YRB under the 1.5 and 2.0°C global warming targets. As the intensity-based index of the Rx5day is mostly related to the flooding hazard (Alexander et al., 2006; Zhang et al., 2018), changes in Rx5day are specifically analyzed hereafter.

2.   Study area, data, and methods

2.1   Study area

The Yangtze River, the longest river in China and the third longest river in the world, originates from the Dangla Mountain Range in the Qinghai–Tibetan Plateau, runs eastward into the flat and low-lying plain, and eventually flows into the East China Sea. The whole basin area is approximately 1.8 million km², ranging from 24°30′N to 35°45′N in latitude and 90°33′E to 112°25′E in longitude, most of which is dominated by the East Asian summer monsoon (EASM) and Indian summer monsoon (ISM) (Freychet et al., 2015). Furthermore, the basin has been identified as a region highly sensitive and vulnerable to global change (Su et al., 2006; Zhang et al., 2018). As illustrated in Fig. 1, the YRB is usually partitioned into the upper, middle, and lower reaches with 11 principal subcatchments.

Download

Larger image

PowerPoint slide

Fig  1.  Yangtze River basin and its subbasins (1: Jinsha River, 2: Min-Tuo River, 3: Jialing River, 4: Wu River, 5: Upper reach, 6: Hanjiang River, 7: Lake Dongting, 8: Lake Poyang, 9: Middle reach, 10: Lower reach, 11: Lake Tai).

2.2   Climatic data of precipitation

The observed daily precipitation data from 1961 to 2014 are from China’s Ground Precipitation 0.5° × 0.5° Gridded Dataset (V2.0) (http://data.cma.cn/), which is produced by interpolating records at 2472 climatic stations with a thin plate spline approach. The 1961–2014 dataset is used as a reference for bias correction and evaluation of GCM projections, and the 1986–2014 dataset is used to represent the present-day climate.

The projected daily precipitation data from 10 CMIP6 models were retrieved through the data portals of Earth System Grid Federation (ESGF; https://esgf-node.llnl.gov/search/cmip6/). The models used, along with their details, are listed in Table 1. The performances of these models are outstanding in East Asia. Generally, the number of models used for ensemble analysis, which is usually at least three, is not certain (Chen et al., 2015; Wang et al., 2017). The projection data are interpolated to a 0.5° × 0.5° grid by using a bilinear scheme for comparison with the observations. Three shared socioeconomic pathways (SSPs) are selected for future projections in CMIP6, including a world of continuing historical trends (SSP2-4.5), a fragmentation world (SSP3-7.0), and a growth-oriented world (SSP5-8.5) (O’Neill et al., 2016).

Table  1.  List of the 10 selected CMIP6 GCMs

Download as CSV

	GCM	Institution and country	Horizontal resolution (lat × lon)	Nominal resolution atmosphere* (km)
1	ACCESS-CM2	CSIRO-ARCCSS, Australia	144 × 192	250
2	BCC-CSM2-MR	BCC-CMA, China	160 × 320	100
3	CAMS-CSM1-0	CAMS, China	160 × 320	100
4	CanESM5	CCCMA, Canada	64 × 128	500
5	CESM2-WACCM	NCAR, USA	192 × 288	100
6	CMCC-CM2-SR5	CMCC, Italy	192 × 288	100
7	INM-CM5-0	INM, Russia	120 × 180	100
8	MPI-ESM1-2-HR	MPI-M, Germany	192 × 384	100
9	NorESM2-MM	NCC, Norway	192 × 288	100
10	UKESM1-0-LL	MOHC, UK	144 × 192	250
*Data from https://wcrp-cmip.org/cmip-model-and-experiment-documentation/

 | Show Table

DownLoad: CSV

2.3   Time to reach the 1.5 and 2.0°C warming targets

Generally, the time series of the global mean temperature were first smoothed by a 20-yr moving average, and then the time to reach a specific warming target was set as the first 20-yr period when the temperature increase (relative to 1850–1900) exceeded the target temperature. The respective central years of a period for each warming level, model, and SSP scenario were derived using the procedure of Hauser et al. (2019) and then averaged for the 10 selected CMIP6 models, which were extended to 20-yr time slices for statistical stability. Under SSP2-4.5, SSP3-7.0, and SSP5-8.5, the time intervals at which the 1.5°C target was reached were 2023–2042, 2022–2041, and 2018–2037, respectively, and those for the 2.0°C target were 2044–2063, 2036–2055, and 2032–2051, respectively.

2.4   ERCDFm bias correction

To address the systematic bias of the ensemble members, we conducted bias correction with the ERCDFm method for each season and each grid for both historical and scenario periods (Li et al., 2010; Pierce et al., 2015). Notably, prior to applying the ERCDFm correction, a cutoff threshold was set to modify the wet-day frequency to cope with the drizzle effects in the GCM projections (Pierce et al., 2015).

The gamma distribution is generally used to fit the CDFs and the quantile functions (inverse CDFs) of precipitation series. To effectively adjust the upper tail of precipitation, a piecewise bias correction is applied, which uses the 95th percentile of the precipitation distribution as a threshold, and the CDFs below and above the threshold are fitted by a gamma distribution separately. In view of the differences between the CDFs under historical and future climate conditions, the corresponding transfer functions between the model outputs and the observations are also different for the two periods, expressed as

where $x_{\mathrm{m}-\mathrm{h}}$ and $x_{\mathrm{m}-\mathrm{p}}$ are the climate variables of the raw model m in the historical (h) and projection (p) periods, respectively; $x_{\mathrm{m}-\mathrm{h}\_\mathrm{a}\mathrm{d}\mathrm{j}\mathrm{u}\mathrm{s}\mathrm{t}}$ and $x_{\mathrm{m}-\mathrm{p}\_\mathrm{a}\mathrm{d}\mathrm{j}\mathrm{u}\mathrm{s}\mathrm{t}}$ are the corresponding bias-corrected values of the two periods; $F_{\mathrm{o}\mathrm{h}}^{-1}$ and $F_{\mathrm{m}\mathrm{h}}^{-1}$ are the quantile functions, respectively, for the observations (o) and the model projection (m) in the historical period; and $F_{\mathrm{m}\mathrm{h}}$ and $F_{\mathrm{m}\mathrm{p}}$ are the CDFs of model m in the historical and projection periods. The right term in Eq. (2) is the multiplicative factor between the observations and projections in the historical period.

2.5   BMA approach

The BMA method has been used in multimodel ensembles to reduce uncertainties among models and to attain superior performance in simulating precipitation extremes. BMA weights are equal to posterior probabilities of the models generating the simulations, reflecting each model’s relative contribution to predictive skill in the historical period. The BMA output produces an improved probability density function (PDF), which is the weighted average of PDFs centered on the individual bias-corrected simulations. Given the observations (training data) $y_{\mathrm{o}\mathrm{b}\mathrm{s}}$, and M models $y_{1},\mathrm{ }\dots ,y_{M}$, the combined simulation $P\left(y|y_{\mathrm{o}\mathrm{b}\mathrm{s}},y_{m}\right)$ is

where $P\left(y|y_{m}\right)$ is the conditional PDF based on the bias-corrected model $\mathrm{\mathit{y}}_{m}$ alone; ${w}_{m}$ is the weight of model ${y}_{m}$, which is obtained from the normalization of posterior probability $P\left(y_{m}|y_{\mathrm{o}\mathrm{b}\mathrm{s}}\right)$ to ensure that the weights add up to 1; and $P\left(y_{m}|y_{\mathrm{o}\mathrm{b}\mathrm{s}}\right)$ represents the probabilities of the conditional PDF from model ${y}_{m}$ recreating the realistic climatic series and is set statistically in proportion to the product of the likelihood function and the prior probability according to Bayes’ theorem, $P\left(y_{m}|y_{\mathrm{o}\mathrm{b}\mathrm{s}}\right)\propto P \left(y_{\mathrm{o}\mathrm{b}\mathrm{s}}|y_{m}\right) P\left(y_{m}\right)$.

The conditional PDF is generally assumed to be normally distributed, centered at $a_{{m}}+b_{{m}}{{y}}_{{m}}$ instead of the original model outputs ${y}_{m}$, to play a role similar to that of bias correction (Raftery et al., 2005; Zhao et al., 2020). However, model weighting is not a standard procedure in climate modeling, and there is no best scheme thus far (Parrish et al., 2012). Here, by setting a normal distribution parameter, b_m is removed to minimize the interference on the original trends of the model data series, then

where parameters ${a}_{m}$ and ${\sigma }_{m}$ are estimated from Markov chain Monte Carlo (MCMC) sampling. In the process of MCMC iteration, the likelihood function $P\left(y_{\mathrm{o}\mathrm{b}\mathrm{s}}|y_{m}\right)$, defined as the probability of training data with given parameters, can be taken as a function of ${a}_{m}$ and ${\sigma }_{m}$, which obeys a uniform distribution with $a_m \sim U (-\mathrm{10,10})$ and ${\mathrm{\sigma }}_{m}\sim U \left(\mathrm{0,50}\right)$, respectively.

Thus, a deterministic prediction can be obtained from

The BMA is applied to estimate the maximum consecutive 5-day precipitation (Rx5day), an index representing the intensity of extreme precipitation.

The weights are determined by three kinds of factors, namely model system bias, differences of temporal variations, and variability in tendency. After bias correction, the effects of the first two are almost removed. Then, the weights after correction will be more determined by the trends of each projection.

The BMA weights were determined grid by grid in each season to explore the spatial characteristics of model performance and to determine whether the “best” model changes with season. First, we derived the seasonal values of Rx5day from the 54-yr (1961–2014) daily series of GCM projections. Then, the whole BMA calculation was implemented in RStan [R interface to the Stan language that provides Bayesian inference through MCMC methods (Carpenter et al., 2017)]. The output of BMA is a PDF, which is a weighted average of PDFs centered on the bias-corrected projections of precipitation and extreme indices, in which the overall weight was derived as a product of the BMA skill weights and the independence weights.

2.6   Independency test among models

In spite of the advantage of BMA, it does not provide a mechanism for reducing the effect of model replication. The strategies to avoid model replication in this study are in two sides:

First is to choose as different as possible candidate models; second is to use the distance matrix (Sanderson et al., 2017) to judge the model replication. The time series of the observations and model simulations from 1961 to 2014 were averaged to form a seasonal climatology. The Euclidean distances were calculated as the area-weighted RMSE between the two model projections (or model and observation). Then, a matrix corresponding to one variable is normalized by a mean pairwise intermodel distance. The independencies of the GCM projections expressed as the normalized Euclidean distance are presented in Fig. 2.

Download

Larger image

PowerPoint slide

Fig  2.  Diagram of the intermodel Euclidean distance matrix for the 10 CMIP6 GCM projections and observations. Each color grid represents a normalized Euclidean distance between two model projections. All distances are aggregated over the four seasons.

Generally speaking, the replication among the 10 chosen models is low. The two smallest intermodel pairwise distances are 0.47 (CESM2-WACCM and NorESM2-MM) and 0.51 (CESM2-WACCM and CMCC-CM2-SR5). The CESM2-WACCM and NorESM2-MM models present the lowest uniqueness weights. Relatively, ACCESS-CM2 and CanESM5 are more independent.

2.7   Model evaluation

To evaluate the skill of the BMA method to simulate Rx5day over the historical period, four statistical indices were considered, namely, relative bias for the deviation from the mean, RMSE for the overall accuracy, interannual variability skill score (IVS) for the long-term variation, and Taylor skill score (TS) for the similarity between the trends of observations and simulations. The IVS is calculated by the standard deviation (Chen et al., 2011):

where ${\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{s}\mathrm{i}}$ and ${\mathrm{S}\mathrm{T}\mathrm{D}}_{\mathrm{o}\mathrm{i}}$ denote the interannual standard deviation of the simulated and observed data at the ith grid, respectively, and N is the total number of grids. A small IVS indicates low interannual variability.

As a combined index considering both the spatial correlation and the amplitude of variability (Taylor, 2001), TS is computed as Sen’s slope (Sen, 1968) as follows:

where COR is the spatial correlation coefficient (Gao et al., 2012); $x_i$ and ${{y}}_{{i}}$ are the observed and simulated trends, respectively, at the ith grid; $\bar x$ and $\bar y$ are the spatial means of the observed and simulated trends, respectively; ${\mathrm{C}\mathrm{O}\mathrm{R}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum COR, which is set as 0.999 (Zhu et al., 2020); and ${\mathrm{S}\mathrm{T}\mathrm{D}}_x$ and ${\mathrm{S}\mathrm{T}\mathrm{D}}_y$ are the standard deviations of the observed and simulated spatial patterns, respectively. As defined, TS $\in \left[\mathrm{0,1}\right]$, which is associated with a better ability to simulate trends from 0 to 1.

2.8   Risk analysis

The risk ratio is calculated on the basis of extreme value theory (EVT), which focuses on the exceedance probability rather than changes in intensity at different return periods and can also be referred to as the probability ratio (PR) (Kharin et al., 2018). We applied the generalized extreme value (GEV) distribution from EVT, which was fitted to the seasonal Rx5day generated from the BMA approach. The data from all three scenarios were mixed into a single sequence for each warming target to integrate all the information.

The GEV distribution function $G\left(z\right)$, as introduced by Jenkinson (1955), is calculated as

where $\mu$, $\sigma$, and $\xi$ are the location, scale, and shape parameters, respectively. To reduce the noisy sampling variability, the GEV parameters are spatially smoothed by a simple procedure (Kharin and Zwiers, 2005), which considers spatial dependency and strength from eight surrounding neighbors of the target cell. The return level $z_T$ is associated with the return period T, which is the threshold likely to be exceeded in a year with a probability of 1/T and may then be solved by inverting the GEV distribution function [$G\left(z_T\right)=1-1/T$].

Three different return values (20-, 50-, and 100-yr) from the present-day period (1986–2014) are used as a baseline to represent different levels of risk, and the changes in the occurrence probabilities of those exceeding the baseline under the 1.5 and 2.0°C warming targets are quantified by using the PR [defined as in Fischer and Knutti (2015), $\mathrm{P}\mathrm{R}=P_{1}/P_{0}$, where $P_{0}$ is the occurrence probability of exceeding a certain percentile in the GEV distribution during the preindustrial period and $P_{1}$ is the occurrence probability of exceeding $P_{0}$ under the warming targets]. Here, we estimated $P_{0}$ with a present-day period (1986–2014) instead of a preindustrial period. To guarantee the reference period being steady, the precipitation is detrended in the period of 1986–2014. Then the detrended series is considered as stationary for calculating the return period. Specifically, the PR represents a factor by which the probability of an extreme event changes from the current climate to a future scenario.

3.   Result analysis

3.1   Evaluation of precipitation bias correction

The spatial distribution of multiyear mean observed Rx5day illustrates a decreasing gradient from southeast to northwest in the YRB (Fig. 3a), which is similar to the spatial pattern of annual precipitation (Li X. et al., 2021). The relative biases between the observed and GCM-projected Rx5day are remarkable, ranging from −30% to 200%, with an average of 22.4% (Fig. 3b). Rx5day is overestimated in the western part of the YRB, whereas it is underestimated in the middle and eastern parts.

The wet biases are effectively reduced after correction with the ERCDFm through which the areal mean biases of annual precipitation amount decrease from 45% to −3.2% (not shown), and the biases of Rx5day decrease to 1.3% (Fig. 3c).

Download

Larger image

PowerPoint slide

Fig  3.  Spatial distributions of annual Rx5day during the historical period (1961–2014) for (a) multiyear means of observations, (b) biases of the multiyear means of raw projection with AM ensemble, (c) biases of the multiyear means of the bias-corrected projection with AM ensemble, (d) trends of the observations, (e) trends of the raw projection with AM ensemble, and (f) trends of the bias-corrected projection with AM ensemble; 95% confidence is shown with the stippling.

As reported by Pan et al. (2016), CMIP5 models may reproduce the spatial distribution of precipitation well, but they fail to project spatial patterns of precipitation trends. As shown in Figs. 3d–f, the area-average trend of Rx5day is 0.82 mm (10 yr)⁻¹, while the average trends of Rx5day for the raw projections and bias-corrected AM ensemble are −0.02 and 0.42 mm (10 yr)⁻¹, respectively. Across the whole basin, the observed annual precipitation trends are positive in 58% (392 out of 674) of the grids. The observed Rx5day trends in 16% of the grids are significant (p < 0.05). Only 7% of the grids exhibited significant trends (p < 0.05) for simulated Rx5day, among which the trends of Rx5day in 19% of the grids are significantly positive, and only 6% of the grids exhibited significant decreases (p < 0.05) after bias correction. This indicates that in our study region, the trend has not yet been improved in CMIP6 projections. ERCDFm correction can partially improve the trend projection, but the bias of trend is still large.

3.2   Seasonal weights of the BMA

The BMA weights can reflect the model performance. As shown in Fig. 4, the pointwise optimal models illustrate distinctive spatial and temporal heterogeneity. Optimal models show an obvious “clustered” spatial pattern, indicating that the model projections may be sensitive to variations in certain terrains and other land surface physical conditions. However, these “clustered” characteristics are not consistent in the four seasons, confirming that performances of the GCM models vary not only spatially but also temporally in the study basin.

Download

Larger image

PowerPoint slide

Fig  4.  Spatial distributions of the pointwise optimal models based on BMA weights calculated from Rx5day for each season, with the areal mean weights of the superior models (weights > 0.1) shown in the bottom-left corner.

The performances of individual CMIP6 GCMs and MME using the AM and BMA methods are shown in Fig. 5. The models with higher BMA weights are likely to produce lower RMSEs or relative biases, such that the INM-CM5-0 model which receives the highest BMA weight in spring (Fig. 4), has the lowest RMSE (27.6 mm) and relative bias (0.6%) in spring (Fig. 5). MME can minimize the weakness of each single model and produce better performance than all the single participating models, with generally lower relative bias and smaller RMSE. The AM ensemble shows positive biases in winter (18.3%) and spring (0.4%) but negative biases in autumn (−3.0%) and summer (−1.1%). The BMA ensemble presents similar or lower biases than the AM ensemble mostly, especially in winter. For interannual variability, the IVS of the AM ensemble is 7.0 on average for all the four seasons, which is much greater than 0, indicating less ensemble efficiency. The reason is that the AM takes into account every data point in averaging and thus is conducive to underdispersion in the ensemble forecasts (Raftery et al., 2005), especially for data series with large variability, such as extreme precipitation events. However, the BMA approach shows its superiority in adjusting underdispersion by reducing the weights of poorly performing models with low IVS (~1.2), but at the expense of a slightly higher RMSE. On average, across the four seasons, the BMA ensemble outperforms the AM ensemble in reducing the relative bias (58%) and IVS [(7 − 1.2)/7 = 83%].

Download

Larger image

PowerPoint slide

Fig  5.  Relative bias, RMSE, and IVS between bias-corrected (each model, AM, and BMA) and observed Rx5day during 1961–2014 for each season.

As shown in Fig. 6, the uptrends of the observed Rx5day are much greater in summer and winter than in the other two seasons, in which the prominent tendencies are mainly distributed in the middle and lower reaches of the basin. Moreover, slightly decreasing trends are found in autumn and spring, some of which are significant (p < 0.05) in the middle reach. There were some marked differences between the bias-corrected and observed trends in both the areal mean and spatial distribution. Represented by the spatial similarity index between the observed and simulated trends, the TS scores of the AM ensemble ranged from 0.03 to 0.16, while those of the BMA ensemble ranged from 0.12 to 0.43 in the four seasons. The seasonal average TS scores of the BMA ensemble are 2.7 times greater than those of the AM ensemble, demonstrating a stronger spatial correlation between the observations and BMA simulations. Despite the inherent deficiencies in GCMs for projecting extreme precipitation trends, the BMA approach can effectively improve the projection of Rx5day trends.

Download

Larger image

PowerPoint slide

Fig  6.  Observed and simulated trends of the seasonal bias-corrected Rx5day during 1961–2014 by the AM and BMA approaches.

3.3   Projected changes in the risk of extreme precipitation

As shown in Fig. 7, the seasonal average change in Rx5day is 1.9% [i.e., (7.5% + 12.8% + 18.9% + 6.8%)/4 − (6.2% + 12.2% + 15.4% + 4.6%)/4] greater than that in the raw GCM projections under 2.0°C warming under the SSP2-4.5 scenario by using AM. By BMA, this change becomes 0.8% [i.e., (9.7% + 12.5% + 11.0% + 8.4%)/4 − (6.2% + 12.2% + 15.4% + 4.6%)/4], which is 1.1% lower than that by the bias-corrected AM ensemble. Both tiny values of 1.9% and 0.8% indicate that the ERCDFm bias-correction process preserves the changing signal of areal mean extreme precipitation well. Relatively, BMA preserves it slightly better (the difference being 0.8%) than AM (1.9%).

Download

Larger image

PowerPoint slide

Fig  7.  Comparison of the raw AM (first column), bias-corrected AM (second column), and bias-corrected BMA (third column) approaches in terms of changes in Rx5day under the 2.0°C warming target under the SSP2-4.5 scenario relative to the present-day climate.

It is seen that more notable spatial heterogeneity is found in the spatial patterns of the BMA results than of AM results. This can be explained that in BMA the weights of the models are different when the model “meets” the data from different underlying geographical unions with different terrain and climate, as shown in Fig. 4, a feature also found by Huang et al. (2024). As a giant basin, Yangtze River basin contains different geographical units. This superiority of the BMA in merging ensemble members to capture local climatic behaviors poses its higher potential to be applied in the assessment of regional climate change. So below we will just use BMA for the scenario projections.

As illustrated in Fig. 8, consistent increases in the annual Rx5day are projected under all SSP scenarios across the YRB, with 9.8% and 15.5% average increases under the 1.5 and 2.0°C warming targets, respectively. From a seasonal perspective, the average changes are more pronounced in summer than other seasons, which are followed by autumn, spring, and winter. The corresponding increases in Rx5day due to additional 0.5°C warming are 8.3% K⁻¹ in spring, 8.5% K⁻¹ in summer, 12.4% K⁻¹ in autumn, and 19.6% K⁻¹ in winter. Although the uncertainties of projections are associated with the magnitudes of changes under the 1.5 and 2.0°C warming targets, the positive responses of Rx5day represent strong signals of precipitation extremes under climate change scenarios in all seasons.

Download

Larger image

PowerPoint slide

Fig  8.  Responses of Rx5day to different SSP scenarios simulated by the BMA ensemble under (a) 1.5°C and (b) 2.0°C global warming targets and (c) their differences over the YRB. The response is represented with relative change, which is equal to [(each of the scenario − present_day)/present_day] $\times$ 100%.

As estimated from the GCM projections with BMA, the areal mean probabilities of 20-yr return events may increase by 2.15 and 3.05 under the 1.5 and 2.0°C warming targets, respectively, averaged over the four seasons, with the greatest increase occurring in summer (Fig. 9). By using a value of 2.5 as the PR threshold, the subcatchments facing high risks with PRs exceeding the threshold under both the 1.5 and 2.0°C warming targets were identified, as shown in the third column of Fig. 9. The risks of extreme precipitation are increasing in summer and autumn over the northwestern part of the YRB (catchment 1, 2, 3, ref. Fig. 1), while increasing risks are projected in the spring and winter seasons in the piedmont of the Qinghai‒Tibetan Plateau (Sichuan Basin, catchment 3).

Download

Larger image

PowerPoint slide

Fig  9.  PRs for seasonal Rx5day from the present-day climate to the 1.5°C (first column) and 2.0°C (second column) warmer climates and areal means of Rx5day and PR over each of the 11 subcatchments in the YRB (third column). PR is calculated for a 20-yr return period with uncertainty ranges marked.

The subcatchment-mean probability ratios of those extreme events are almost all beyond 1 for the return periods of 20, 50, and 100 yr, indicating that additional risks will be brought up by the additional 0.5°C warming (Fig. 10). The projections by the BMA demonstrate that the areal mean risk of Rx5day will increase by factors of 1.8 in spring, 1.5 in summer, 1.7 in autumn, and 2.4 in winter for the 20-yr return period. Spatially, the distributions of high-risk zones are quite scattered, except for the catchments in the middle to lower reaches in winter (Figs. 10j–l). Although the spatial patterns of high-risk regions are similar for different return periods, the magnitudes of the increases in risk are more noticeable for extreme events with longer return periods.

Download

Larger image

PowerPoint slide

Fig  10.  Changes in the PR from 1.5 to 2.0°C warmer climates for seasonal Rx5day that are expected occurrence once every 20, 50, and 100 yr in the present-day climate.

We further compared the areas that experienced extreme disastrous precipitation events (defined as the current 20-yr return level) under the two warming targets. The areas were spatially aggregated to derive the total area exposed to extreme precipitation events and the fractions to the whole basin for assessment of the impacts exerted by the additional 0.5°C of warming from 1.5 to 2.0°C. As shown in Table 2, the exposed areas increase consistently with warming for all the seasons under the three SSP scenarios except for SSP5-8.5 in autumn. From the perspective of seasonal distributions, there are 10.9% (10.5%–11.7%), 8.7% (4.0%–12.7%), 8.3% (−6.5% to 16.5%), and 20.3% (14.5%–27.4%) expansions of the exposure areas for the 20-yr return events in the four seasons, respectively. Thus, if global warming is limited to 1.5°C instead of 2.0°C, there will be much less exposure to disastrous extreme events, especially in the winter and spring seasons.

Table  2.  Fraction increments of the exposed area (%) by additional 0.5°C warming (from 1.5 to 2°C) for Rx5day events beyond the present-day 20-yr return level

Download as CSV

	Spring	Summer	Autumn	Winter
SSP2-4.5	10.5	4.0	16.5	27.4
SSP3-7.0	11.7	12.7	15.0	14.5
SSP5-8.5	10.5	9.3	–6.5	19.0
Average	10.9	8.7	8.3	20.3

 | Show Table

DownLoad: CSV

4.   Discussion

4.1   Comparison of empirical distribution-based v.s. gamma distribution-based bias correction

We noticed from documented researches (e.g., Zhu et al., 2022) that empirical distribution-based bias correction outperforms gamma distribution fitting. To evaluate the performance of empirical PDF and gamma function, we applied the empirical cumulative PDF and gamma function respectively to correct the historical daily precipitation of BCC-CSM2-MR model (1961–2014) on two grids respectively as an example. Since we focused on extreme precipitation, the correction was carried out by using two different daily precipitation levels, namely the daily precipitations above and below 95% percentile being divided into two groups. The results showed that the relative biases of the corrected precipitation with empirical function were reduced from 30%–87% to −1.9%–0.7% for daily precipitation, and from 33%–99% to −4.7%–8.9% for Rx5day. The relative biases of precipitation with gamma function were reduced to −2.7%–0% for daily precipitation and to −6.1%–7.9% for Rx5day. The determinant coefficient between the corrected daily precipitations of these two methods is above 0.98 and the differences of relative biases are less than 2%.

From above comparison, it is found that for heavy rainfall, empirical distribution-based bias correction also slightly outperforms gamma distribution fitting in our case. However, given that the difference is not obvious, considering the convenience of computation, we remain using gamma function for Rx5day indices. It is worthy to consider this carefully in the future research.

4.2   The uncertainty of spatial resolution

The spatial resolution of the chosen GCMs varies as shown in Table 1 in a different combination of the number of grids in latitudinal and longitudinal directions. The normal spatial resolution or weighted mean diagonal grid box distance (Klaver et al., 2020), for variables in atmosphere is from 100 to 500 km. In this article all the GCM model simulation outputs for variables in atmosphere were interpolated into 0.5°, to be matched with the resolution of the observation. That is, the goal resolution of statistical downscaling is 0.5°.

There are other options of statistical downscaling. For example, if the goal resolution is set to 1°, all the GCM model simulation outputs for variables in atmosphere will be interpolated into 1° and the observation dataset will be upscaled to 1°. It is expected that there will be uncertainty by using different spatial resolutions.

However, for some variables, such as temperature, documented research showed that it is not sensitive to spatial resolution, in which for regridded daily mean temperatures, the PDFs were shown very similar to that from the original data under the spatial resolutions of 0.5°, 1°, and 2°, in different regions (Sun et al., 2015). By comparing some CMIP6 models in high and low resolutions with the observations, Dong and Dong (2021) found that improving resolution does not necessarily improve the model’s simulation of extreme precipitation, showing that the accurate representations of physical processes are still a challenge for current global climate models.

Thus in this article, by assuming the model biases caused by downscaling smaller than the errors in the models (e.g., Li T. et al., 2021), no further analysis by using different spatial resolution has been done. The target resolution interpolated is determined only with the consideration to the resolution of observation dataset. The bias is expected to be adjusted by using the method of ERCDFm.

4.3   Why shall we still need BMA given that the bias has been corrected?

The uncertainty regarding long-term trends stems from the projections of GCMs that integrate anthropogenic influences and internal climate fluctuations (Deser, 2020). Bias correction method can well correct the systematic bias of GCM projections. However, it can partially improve the trend projection. From Fig. 3, it is seen that ERCDFm correction can perfectly correct the bias of mean value of Rx5day. However, it can just partially improve the trend projection. The bias of trend after the correction is still large.

Interannual fluctuations and decadal trends in extreme precipitation are intrinsic characteristics of GCMs. To preserve the intrinsic characteristics of the projections in the BMA procedure and to minimize BMA interference on the original trends of the Rx5day series, the regression coefficients of the GCM output items are omitted by setting a conditional PDF (e.g., Zhao et al., 2020). Therefore, for GCMs with more efficient performance in simulating fluctuations and trends, it is relatively easier to obtain more important weights (Knutti et al., 2017). From Fig. 5, BMA shows its superiority in reflecting interannual variability with the IVS as low as ~1.2, compared to 7.0 by AM ensemble. This supports us to use BMA for better projection in addition to bias correction.

Furthermore, as shown in Fig. 7, more notable spatial heterogeneity is found in the spatial patterns of the BMA-with-bias-correction results than that of the AM-with-bias-correction results. This can explain that in BMA the weights of the models are different when the model “meets” the data from different underlying geographical unions with different terrain and climate, as shown in Fig. 4. The Yangtze River basin is a giant basin containing different geographical units. The superiorities of the BMA in merging ensemble members to capture local climatic behaviors pose its higher potential to be applied in the assessment of regional climate change. Basher et al. (2020) and Zhao et al. (2020) also showed that BMA outperforms the AM approach in simulating the mean state and interannual variability in extreme precipitation.

4.4   More spatial heterogeneous patterns of Rx5day after using BMA

We noticed that there is more spatial heterogeneity in BMA projections as shown in Figs. 6 and 7. The more significant spatial inconsistency of Rx5day with BMA method is related with Rx5day projections by the 10 GCMs. Since the BMA weights of each GCM are different over the sub-basins, Rx5day projected by each model is different. The output of BMA has reasonably illustrated remarkable spatial heterogeneity in the projections after combining the weights and Rx5day for each projection. This phenomenon is acceptable as the Yangtze River basin is characterized with broad climate zones, terrain, as well as land use/land cover. The precipitation pattern is highly spatial heterogeneous for each GCM projection, especially for the extreme events.

4.5   Increasing Rx5day in the scenarios

Consistent increasing trends in precipitation extremes are projected by the CMIP6 ensemble, which are associated with the intensified thermodynamic and dynamic processes of climate systems under global warming (Xu et al., 2009; Pan et al., 2016; Vittal et al., 2016). In the thermodynamic process, this increase is due to the improved water-holding capacity of the warmer air (Christensen and Christensen, 2003). However, the dynamic processes, representing changes in large-scale oceanic and atmospheric circulation oscillations, are the primary dominant factors in shifts in regional extreme precipitation characteristics. In YRB, the varying transport routes of moisture caused by changes in monsoon activities, especially the EASM, have played critical roles in determining the spatial and temporal patterns of summer precipitation (Chen et al., 2014; Li X. et al., 2021).

In general, the GCM projections illustrate an intensified land–ocean contrast in the scenarios, which contributes to stronger EASM circulation patterns and more unstable local atmospheric stratification and leads to greater risks of extreme precipitation in the midlatitude than in the low latitude (Chen, 2013; Wang et al., 2018). The risks of extreme precipitation in the middle and lower reaches of the YRB will increase at ratios of 2.2 (2.8) and 2.9 (3.9) under the 1.5°C (2.0°C) warming target in summer, respectively, which may threaten the local infrastructure for flood inundation protection.

Further exploration will be the critical factors controlling monsoon variability and enhance the understanding of its impact on regional extreme precipitation in large river basins.

5.   Conclusions

A reliable analysis of the changes in the magnitudes and risks of seasonal extreme precipitation in the YRB under the 1.5 and 2.0°C warming scenarios was conducted. For this purpose, the BMA and AM approaches were applied to 10 bias-corrected CMIP6 GCM projections. To effectively adjust the wet biases and preserve the long-term tendency signals for extreme precipitation, piecewise ERCDFm bias correction was applied to the GCM projections. To further reduce the uncertainties of the projections, a BMA-based approach was designed to minimize interference with the original trend of extreme precipitation when optimizing the weights of each GCM. Generally, the BMA outperforms the AM ensemble by an overall 58% reduction in relative bias, an 83% reduction in IVS, and a competitive RMSE for precipitation extremes. The tendencies of extreme precipitation from the BMA approach are considered more reliable.

The BMA projections under the 1.5 and 2.0°C warming scenarios revealed that the maximum area-mean risk of Rx5day may occur in summer, followed by autumn, spring, and winter. Moreover, some regions at high risk of extreme precipitation under warming scenarios were identified spatially, with the northwestern parts of the YRB experiencing the greatest risk in summer and autumn and the Sichuan Basin experiencing greater risk in spring and winter. For the additional 0.5°C warming from 1.5 to 2.0°C, the seasonal risk ratios of Rx5day increase for extreme events with 20-, 50-, and 100-yr return periods. It is hoped that the shifts in spatial and temporal patterns of extreme precipitation will propel local administrations to improve flood protection infrastructure and social awareness education.