A two-stage Bayesian multi-model framework for improving multidimensional drought risk projections over China

Understanding future drought risk is a prerequisite for developing climate change adaptation strategies and for enhancing disaster resilience. In this study, we develop multi-model probabilistic projections of multidimensional drought risks under two representative emission scenarios (RCP4.5 and RCP8.5) through a copula-based Bayesian framework. An ensemble of five regional climate simulations, including four from the CORDEX East Asia experiment and one from the 15 Providing REgional Climate Impacts for Studies (PRECIS) simulation, is used to project future changes in hydroclimatic regimes over China. A new Bayesian copula approach is introduced to uncover underlying interactions of drought characteristics and associated uncertainties over 10 climate divisions of China. The proposed Bayesian framework explicitly addresses the cascade of uncertainty in high-resolution projections of multidimensional drought risks. Our findings reveal that precipitation and potential evapotranspiration (PET) are projected to increase for most areas of China, while increasing 20 radiative forcing is expected to amplify the increase in PET but does not cause significant changes in the precipitation projection. In addition, the drought duration and severity are projected to substantially increase for most areas of China. The estimated drought risks in China are expected to become more than double under both emission scenarios. The extreme droughts are projected to intensify in terms of frequency and associated risks as the radiative forcing increases.

3 between drought characteristics regardless of their marginal distributions (AghaKouchak et al., 2014;Ganguli and Reddy, 2014;Masud et al., 2017;Sadegh et al., 2018;Salvadori and De Michele, 2004;Salvadori et al., 2016). For example, Xu et al. (2015) considered the spatial extent of droughts in the copula-based multivariate drought frequency analysis in Southwest China. Zhang et al. (2019) used copula and the convection-permitting climate simulations to assess climate change impacts on the multivariate drought evolution over South Central Texas. One of the most important variables derived by multivariate 70 drought frequency analysis is the drought return period, which represents the average time between drought episodes and thus quantifies drought risks (Kwon and Lall, 2016;Masud et al., 2015). For example, AghaKouchak et al. (2014) indicated that the 2014 California drought can be a 200-year extreme event when considering the combined effects of low-precipitation and high-temperature conditions. Liu et al. (2016) concluded that the average multivariate return period of extreme droughts in China during 1961−2013 was 42.1 years. Previous studies, however, fail to guarantee the global optimization of copula 75 parameters by using the frequentist approach, leading to a potential bias in drought return periods. Another limitation of previous studies is that they fail to quantify the underlying uncertainties of copula parameters. Such uncertainty is considerably large since the samples of drought episodes are typically limited, and ignoring the uncertainty diminishes the scientific credibility in drought assessments (De Michele et al., 2013;Sadegh et al., 2017). Therefore, it is necessary to explicitly address the uncertainty in copula-based drought risk assessments for advancing our understanding of complex mechanisms and 80 potential impacts of droughts.
The objective of this study is to develop a two-stage Bayesian multi-model framework to improve the multidimensional drought risk assessments in a changing climate. Specifically, an ensemble of five regional climate simulations, including four from the CORDEX East Asia experiment and one from the Providing REgional Climate Impacts for Studies (PRECIS) simulation will be used to improve the performance of climate simulations in China based on the Bayesian model averaging 85 (BMA) approach. Drought episodes will be detected using the Standardized Precipitation-Evapotranspiration Index (SPEI) in 10 climate divisions of China (Vicente-Serrano et al., 2010). Drought risks will also be quantified using the joint return period of duration and severity calculated by a Bayesian copula approach. The uncertainties in BMA and copula parameters will be addressed within a Bayesian framework, leading to probabilistic hydroclimatic projections and drought return periods. The Climatic Research Unit (CRU) dataset will be collected to evaluate the BMA-based simulations of precipitation and PET. 90 This paper is divided into four sections. Section 2 will describe models, algorithms, and datasets used to perform Bayesian multi-model climate simulations and multivariate drought risk projections. Section 3 will systematically evaluate the BMAbased hydroclimate simulations and assess climate change impacts on multidimensional drought risks. Finally, Section 4 will provide a summary and conclusions of this study.

Two-stage Bayesian multi-model framework
To improve multidimensional drought risk projections, we propose a copula-based two-stage Bayesian multi-model framework, as shown in Fig. 1. The framework explicitly uncovers the cascade of uncertainty in multidimensional drought projections. The first step is to perform ensemble climate simulations, including four climate simulations available in the CORDEX East Asia experiment and one PRECIS simulation over China. The second step is to uncover the uncertainty in 100 model weights (Figs. 1a−e) through the MCMC-based BMA climate simulation, leading to probabilistic hydroclimatic projections in each grid cell (Figs. 1f and 1g). Details of ensemble climate projections are given in Section 2.2. The BMA climate simulations will be compared against the AEM simulations. The third step is to propagate the uncertainty in hydroclimatic projections to the uncertainty in drought projections based on drought indices, leading to multiple scenarios of drought characteristics (i.e., duration and severity), as shown in Figs. 1h and 1j. The drought detection is performed using the 105 6-month Standardized Precipitation Evapotranspiration Index (SPEI6) over 10 climate divisions in China (see Fig. 2a). The parameters required in the SPEI6 calculation are estimated based on the historical data, and then used to calculate SPEI6 under future climate. The fourth step is to perform the MCMC simulations for copula parameter inference and uncertainty quantification, leading to the uncertainty in the dependence structure of drought variables (Fig. 1i) and thus the uncertainty in return periods of drought episodes (Fig. 1j). The red "whiskers" in Fig. 1j represent the uncertainty of drought characteristics 110 resulting from the climate projection. Details of the copula-based drought risk assessment are provided in Section 2.3.

Bayesian multi-model climate projection
The PRECIS model developed by the UK Hadley Centre, together with four regional climate simulations from CORDEX available for the East Asia domain, were used to assess the changes in hydroclimatic regimes over China. Specifically, the COnsortium for Small-scale MOdelling in CLimate Mode (CCLM) RCM was used to dynamically downscale four Coupled 115 Model Intercomparison Project Phase 5 (CMIP5) GCMs (CNRM-CM5, EC-EARTH, HadGEM2-ES, and MPI-ESM-LR) in the CORDEX East Asia experiment, while the PRECIS model was driven by the HadGEM2-ES (Huang et al., 2018;Rockel et al., 2008;Zhu et al., 2017). All the five simulations have the same horizontal resolution of about 0.44° × 0.44° (~ 50 km) but differ in the model domain. The computational domain of the PRECIS simulation covers a region with 109 × 88 horizontal grid points with 19 vertical levels in the atmosphere (see Fig. 2a). In comparison, the CCLM model domain is slightly different 120 with 203×167 horizontal grid points (see Fig. 2b). The PRECIS climate simulation covers the historical period  and a future period , while the CCLM climate simulation covers the historical period  and a future period . Future simulations for both PRECIS used in this study and CCLM used in the CORDEX East Asia experiment are forced with two emission scenarios, including RCP4.5 and RCP8.5. The 30-year monthly hydroclimatic variables including precipitation and potential evapotranspiration (PET) for the historical (1975−2004) and future (2069−2098) 125 periods are collected from the five climate projections to assess the impact of climate change on hydrological regimes. The https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.
FAO-56 Penman-Monteith Equation was applied to the calculation of PET, which was suggested to yield more realistic estimates than the temperature-only-based Thornthwaite method (Allen et al., 1998;Dai, 2013;Sheffield et al., 2012).
Bayesian model averaging (BMA), as an effective tool of correcting under dispersion in ensemble climate projections, was used to improve the accuracy of monthly precipitation and PET simulations. Assume that x = x1,…, xK signify the ensemble 130 of all considered climate simulations, and y denotes the climate observations. pk(y|xk) represents the conditional pdf of y given xk. The probabilistic forecast pdf of y for the multi-model ensemble can be expressed as: where wk is the BMA weight of model k in the ensemble. The sum of all wk values is equal to 1 and they are nonnegative, which reflect how well an individual climate simulation matches the observation in the training period. The conditional pdfs, 135 pk(y|xk), are commonly assumed to follow a normal distribution. As a result, the original forecast is usually transformed to the space of normal distribution as: The values for ak and bk are member specific, and they are the linear transformation parameters derived by simple linear regression of observations on each climate simulation in the ensemble, leading to the BMA predictive mean as Eq. (3). 140 BMA has been demonstrated to be a powerful approach to combine an ensemble of climate simulations since it is essentially an "intelligent" weighted average forecast based on the model performance (Raftery et al., 1997;Vrugt, 2016;Vrugt and Robinson, 2007;Wilson et al., 2007;Yang et al., 2011;Zhang et al., 2016). Therefore, BMA was applied to monthly precipitation and PET for each grid cell with CRU's (Climatic Research Unit) gridded monthly precipitation and PET dataset 145 as reference. The CRU dataset is a global gauge-based climate variable product with a 0.5° × 0.5° grid resolution based on thousands of weather stations (Harris et al., 2014).
The BMA weights and the variance σ 2 were estimated in this study using the MCMC simulation instead of the EM algorithm.
The MCMC simulation has been demonstrated to outperform the EM algorithm, which explicitly samples the posterior distribution of the BMA parameters for uncovering the uncertainty associated with model weights and thus improving the 150 reliability of climate projections (Duan and Phillips, 2010;Vrugt, ter Braak, et al., 2008). The MCMC-based BMA simulations were performed with uncertainty ranges, enhancing the reliability and credibility of the multi-model climate projections. The MCMC simulation was implemented using the Differential Evolution Adaptive Metropolis (DREAM) algorithm (Vrugt, 2016).
where p(w, σ 2 ) and p(w, σ 2 |x, y) denote the prior and posterior distributions of BMA weights and variance, respectively. p(x, y|w, σ 2 )  L(w, σ 2 |x, y) denotes the likelihood function; p(x, y) denotes the evidence that acts as a normalization constant, which can be excluded from the Bayesian analysis in practice. Thus, the formulation of Eq. (4) can be simplified as: The likelihood function L(·|·) in the MCMC-based BMA projection is commonly logarithmically transformed to Eq. (6) for numerical stability and simplicity, where n represents the number of observations in the training period.
The prior distribution is set as a uniform prior distribution of w [0, 1] K and σ 2  [0, 3·var(y)]. The MCMC simulation proceeds by running multiple Markov chains simultaneously and proposing a candidate point zp at each step (Vrugt, 2016;Wang and 165 Wang, 2019). The acceptance or rejection of the candidate depends on the Metropolis acceptance probability: where zc represents the current point, and p(‧) represents the probability density. The Markov chain moves to zp or not, depending on whether the candidate point is accepted. The convergence of Markov chains indicates that the MCMC evolution can stop, which is commonly monitored through the multi-chain R̂ diagnostic of Gelman and Rubin (1992). Typically, a R̂-170 statistic value below 1.2 indicates that the posterior distribution converges to the stationary distribution. A more detailed description of the MCMC simulation, together with the DREAM algorithm, is available in Vrugt, ter Braak, et al. (2008) and Vrugt (2016).

Copula-based Bayesian multidimensional drought risk projection
Copulas are multivariate cumulative distribution functions that enable us to link the marginal distributions of multiple random 175 variables together to form the joint distribution (Genest et al., 2007;Zhang et al., 2019). The dependence of drought duration and severity detected by the SPEI6 over each of the 10 climate divisions in China (see Fig. 2a) was thus described using copulas in this study, leading to a bivariate return period of drought episodes. The 10 climate divisions were created based on the long-term mean temperature and precipitation as well as the topography in China. Assume that X = X1,…, Xn denote n random variables, and F1(x1),…, Fn(xn) represent their marginal cumulative distribution functions (CDFs), the joint CDF 180 F(x1,…, xn) can be expressed as Eq. (8) according to Sklar's theorem (Sklar, 1959 where F(x) = P(X ≤ x) and G(y) = P(Y ≤ y) are the marginal CDFs of drought severity and duration, respectively. To identify the marginal CDF of drought characteristics, several types of probability distributions, including Nakagami, exponential, Rayleigh, gamma, inverse Gaussian, t location scale, generalized Pareto, Birnbaum-Saunders, extreme value, logistic, lognormal, Weibull, log-logistic, Rician, generalized extreme value, and normal distributions were included as the CDF candidates. The optimal copula families were chosen from a total of 11 candidates, including Independence, Gaussian,Clayton,190 Frank, Gumbel, Joe, Nelson, Marshal-Olkin, BB1, BB5, and Tawn. Formulas of the copula families are provided in Table 1.
Both the marginal CDF and copula families were selected using the Akaike information criterion (AIC). In addition, a randomization strategy (also known as "Jittering") was used to avoid the potentially adverse impact of repeated drought durations on the bivariate analysis (Chambers et al., 2018;De Michele et al., 2013).
The copula parameters were estimated through the MCMC simulation in a Bayesian framework similar to the BMA 195 parameters, leading to the posterior parameter distribution instead of the deterministic maximum likelihood (ML) estimates.
Here, the Multivariate Copula Analysis Toolbox (MvCAT) was adopted to infer the MCMC-based copula parameters (Sadegh et al., 2017). The log-likelihood function for copula parameter inference in the MvCAT is expressed as: where θ is the copula parameter set; n denotes the total number of observations; σ denotes the standard deviation of 200 measurement error; ỹi denotes the empirical joint probability of observation i calculated using Gringorten plotting position (Gringorten, 1963); yi(θ) is the joint probability of observation i calculated by the parametric copula with the given parameter θ. Different from the BMA parameters, the prior distributions of copula parameters are drawn using Latin Hypercube Sampling (LHS) which is an efficient sampler and has been widely used for implementing robust MCMC simulations (Huang et al., 2018;Stein, 1987;Vrugt, 2016). The Bayesian inference of copula parameter values requires specifying the initial uncertainty 205 ranges, which are provided in Table 1. More details about the MCMC-based inference of copula parameters can be found in Sadegh et al. (2017).
To project the future drought risks, the joint return period of all the episodes in which drought severity (S) and duration (D) exceed their respective threshold is computed using inclusive probability ("OR" and "AND" case) (Salvadori and De Michele, 2004). The drought return period is commonly proportional to the rarity of drought episodes and the relevant losses, 210 and thus climate-induced drought risks can be evaluated by comparing the return periods under past and future climates. The two cases of bivariate return period can be computed using the copula-based approach as: https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.
where μ denotes the average inter-arrival time between the occurrences of drought episodes (Zhang et al. 2017). It should be 215 noted that the return period is not deterministic but probabilistic with uncertainty ranges due to the posterior distribution of BMA weights and copula parameters derived from the MCMC simulation. the corresponding optimal parameters derived using the EM algorithm. In general, an excellent consistency is observed between the MCMC-derived posterior parameter distributions and the EM estimates within the high-density region. More importantly, the MCMC simulation provides uncertainty information on the model weights in the BMA prediction, which 225 provides multiple scenarios that can generate equally good climate simulations. In addition, the PRECIS simulation has the largest contribution (nearly 0.5) in the ensemble simulation to reproducing the temporal pattern of monthly precipitation observation.

Simulation of historical precipitation and PET
Figures 4 and 5 display the spatial distributions of the 30-year annual, winter (December-January-February, DJF), and summer (June-July-August, JJA) mean precipitation and PET, respectively. These spatial distributions are derived from the 230 BMA ensemble simulations and the CRU datasets as well as the absolute model bias generated by the AEM and BMA approaches. In general, the BMA ensemble simulation and the CRU dataset show a similar spatial pattern of the annual, winter, and summer mean precipitation and PET. Compared to the AEM simulations, the BMA ensemble simulations have significantly lower absolute model biases except for the winter mean precipitation over Southeast China. For example, the AEM simulation tends to underestimate the summer precipitation over Southeast China but overestimate over the Tibetan 235 Plateau (Fig. 4k). Such model bias has been largely reduced by the BMA simulation although dry biases remain over Southeast China (Fig. 4l). The improvement of the BMA simulation upon the AEM simulation is more significant for PET than precipitation. For example, a nearly excellent agreement exists between the BMA-simulated summer PET and the CRUderived PET over China (Fig. 5l). In comparison, the AEM-simulated summer PET generally has a positive bias of over 0.8 mm/day over Northwest and Southeast China, as well as a negative bias of more than 0.8 mm/day over the Tibetan Plateau. 240 This indicates that the AEM-based projection of drought risks can be largely overestimated over Southeast China based on the climate simulations currently available in the CORDEX East Asia experiment due to the overestimated evapotranspiration and the underestimated precipitation. Such an overestimated risk can be significantly corrected by using the BMA ensemble simulation.
To further evaluate the performance of AEM and BMA in simulating the annual, summer and winter mean precipitation 245 and PET, the Taylor Diagram is used to quantify the consistency between the patterns from two simulations and the CRU observation (Taylor, 2001). The simulated pattern agrees better with observations if the model has a higher correlation and a more consistent standard deviation (SD) with the observation, as well as it lies nearer the "OBS". Figure 6 presents the relative performance of AEM and BMA in simulating the annual, summer, and winter mean precipitation and PET for 10 climate divisions in China. In general, the BMA simulations have higher CORs than the AEM simulations for both annual and seasonal 250 results. For example, the CORs for the AEM-simulated annual precipitation are below 0.9 for the 10 climate divisions, and nearly all of them are larger than 0.9 for the corresponding BMA simulation (Fig. 6a). The BMA simulation also significantly corrects the variation amplitude of the AEM simulation. For example, the AEM-simulated annual PET for multiple climate divisions has SDs larger than 1.2 while nearly all the BMA-derived SDs lie between 0.8 and 1 (Fig. 6b). In addition, the AEM simulation has a poor performance of reproducing the winter and summer mean precipitation, as well as the summer mean 255 PET since nearly all of their CORs are smaller than 0.6 (see Figs. 6c, 6e, and 6f). But the corresponding CORs for the BMA simulation are increased for all the climate divisions, especially for the summer mean precipitation (Fig. 6e). Therefore, it is evident that the BMA simulation outperforms the AEM simulation in simulating the annual, winter, and summer precipitation and PET for 10 climate divisions in China. It should also be noted that the BMA simulation does not necessarily outperform the AEM simulation in terms of the variation amplitude. For example, although the CORs for the BMA-simulated winter mean 260 PET have a slight increase upon the AEM simulation, the consistency of the SDs with observations is reduced for several climate divisions. However, the AEM simulation generally overestimates the precipitation over the Tibetan Plateau (Divisions 3 and 5) but underestimates the precipitation over the other climate divisions. Such biases can be caused by the convection parameterization schemes used in the RCMs, which misrepresent the diurnal cycle of convective precipitation that is dominant over East China.
BMA can largely reduce the bias of precipitation simulation compared to the AEM simulation, especially for Divisions 3 and 5. For example, the AEM-simulated precipitation has a wet bias of more than 2 mm/day for Division 5 and a dry bias of nearly 270 3 mm/day for Division 8 in July, while the corresponding BMA-simulated biases are reduced to approximately 0.5 and 2.2 mm/day. In addition, the AEM simulations show smaller biases for PET than precipitation, but the biases are significant over several divisions. Specifically, the AEM-simulated PET has a generally negative bias of more than 1 mm/day over the Tibetan Plateau (Divisions 3 and 5) but tends to have a positive bias over Southeast China during the warm season (Divisions 7,8,9,and 10). Such biases can be corrected by the BMA simulation, leading to a nearly excellent agreement between the simulated 275 and observed PET for all the climate divisions. The MCMC-based BMA simulations provide the 95% uncertainty ranges of precipitation and PET due to the sampling of the posterior distributions of climate model weights. Therefore, BMA leads to a better performance than AEM to reproduce the annual cycle of precipitation and PET for 10 climate divisions over China.

Projection of future precipitation and PET
Since the drought occurrence is closely related to precipitation and PET, the BMA-derived projection of future precipitation 280 and PET is vital for assessing the climate-induced changes in drought risks. Figure 9 presents the absolute and relative changes in the 30-year annual, winter (DJF), and summer (JJA) precipitation derived from the BMA projection between past (1975−2004) and future (2069−2098) climates under RCP4.5 and RCP8.5. Note that the presented BMA projections are based on the "best" BMA weights in the MCMC-derived posterior distribution for better visualization. The most significant drying appears over northern Xinjiang, the western Tibetan Plateau, and Sichuan Basin, with a considerable reduction in the amount 285 of annual and seasonal precipitation. The other parts of China are expected to become wet under both scenarios, especially for southern Xinjiang and the northern Tibetan Plateau under RCP8.5 (Fig. 9j). The percentage changes of annual, summer, and winter precipitation will not generally exceed 40%. In addition, there are large discrepancies between summer and winter precipitation changes. For example, the precipitation over southeast coastal areas of China is projected to increase by 1 mm/day in summer but decrease by 0.4 mm/day in winter under both RCP4.5 and RCP8.5 scenarios. It can be seen that the absolute 290 change of precipitation in summer is generally larger than that in winter, but the relative change shows little difference.
Although the absolute change magnitude of summer precipitation generally increases from west to east of China, such patterns are not observed for the relative change. Additionally, the increase of the radiative forcing does not necessarily lead to a significant change in the projected mean precipitation, such as the summer mean precipitation over the southeast coastal areas.
In addition to precipitation, the 30-year annual, winter (DJF), and summer (JJA) mean PET are also compared between 295 past (1975−2004) and future (2069−2098) climates under RCP4.5 and RCP8.5, as shown in Fig. 10. In general, the annual, winter, and summer PET are projected to increase under both scenarios, with an absolute magnitude of less than 0.7 mm/day and a relative magnitude of less than 30%. The absolute increase of PET in summer is remarkably larger than the PET increase in winter, but the relative increase shows no significant difference between summer and winter PET. For example, the projected increase in winter PET does not generally exceed 0.15 mm/day in most areas of China under RCP8.5 (Fig. 10e), but the 300 corresponding increase in summer is generally over 0.3 mm/day and even exceeds 0.6 mm/day in the western Tibetan Plateau, Sichuan Basin, and Northeast China (Fig. 10f). And the relative increase for both summer and winter PET is generally larger than 10% under RCP8.5 for most areas in China (Figs. 10k and 10l). In addition, the degree of absolute and relative changes can also be magnified by the increase in the radiative forcing, especially for the summer mean PET.

Multidimensional drought risk assessment 305
To uncover the interaction of drought characteristics, copulas were used to construct the joint probability distribution between drought severity and duration detected by the SPEI6 for 10 climate divisions in China during 1975−2004. The MCMC simulations were performed for parameter inference and uncertainty quantification of the chosen optimal copula family. Figure   11 presents the marginal posterior distribution of copula parameters derived from the MCMC simulation. The red asterisk in each panel represents the ML estimates derived by the frequentist copula approach. The frequentist copula approach refers to the use of local optimization algorithms (e.g., the L_BFGS-B method) with initial values to estimate parameters (Yan, 2007).
The Marshall-Olkin copula with two parameters was chosen for Divisions 1−3 and 8; the Clayton copula was chosen for Divisions 4−7; the Gumbel copula was chosen for Divisions 9−10. Most of the posterior parameters are well constrained with normal distributions, but some are not, especially for the second parameter θ2 of the Marshall-Olkin copula (e.g., Figs. 11d and 11f), with a nearly uniform marginal distribution. Such unconstrained parameter distributions can be due to the limited 315 samples of drought episodes. In addition, there is generally a plausible consistency between the posterior distribution of copula parameters inferred by the MCMC simulation and the ML estimates from the frequentist approach for most copula families, but divergent parameter estimates exist for several copulas (e.g., Figs. 11c and 11e). Such a divergence does not imply that the frequentist copula approach provides unreliable simulations but indicates that the frequentist approach provides only one plausible estimate. In comparison, the MCMC-derived posterior parameter distribution provides multiple equally good copula 320 simulations, enhancing the reliability of multidimensional drought risk assessments.
To examine the fit quality of copulas, the joint probability derived from the empirical copulas and the parametric copulas are compared against each other, as shown in Fig. 12. The comparisons between the MCMC-based "best" copula and the frequentist copula are distinguished by different colors. The closer the points are to the diagonal in the diagnostic plot, the better the copula fitting is. In general, both the MCMC-based and frequentist approaches provide plausible copula simulations, 325 especially for Divisions 1 and 9. But the frequentist approach tends to underestimate the joint probability compared to the empirical joint probability. Such an underestimation does not necessarily lead to biased copula simulations but can be potentially risky since the frequentist approach fails to guarantee the global optimization for reproducing the joint distribution of observations.
To further uncover the uncertainty in drought risk assessments, Fig. 13 depicts the return period ("AND" and "OR" case) 330 based on drought duration and severity for Division 3 during 1975−2004, derived from the frequentist and Bayesian copulas.
Although Figs. 11e and 12c indicate the difference between the frequentist and Bayesian copula simulations, Fig. 13 generally shows an excellent agreement between the drought return period isolines, indicating that the Bayesian approach uncovers the equifinality of copula simulations. The Bayesian copula also uncovers considerable uncertainty in the drought return period level, especially for those long ones. For example, for a drought episode with a SPEI6 severity of 17.5 and a return period of 335 50 years, the drought duration can reach 8 months at least and 13 months at most considering the uncertainty in copula parameters. Therefore, the Bayesian copula improves drought risk assessments through a robust assessment of multidimensional characteristics of drought episodes and associated uncertainties. Figure 14 presents the comparison of drought severity and duration detected by the 6-month SPEI between the historical 340 (1975−2004) and future (2069−2098) periods over 10 climate divisions in China. In general, both the drought severity and duration are projected to increase over 10 climate divisions. For example, the median drought durations are approximately 2 months over Divisions 4−6 for the historical period (1975−2004) and are projected to increase to 5 months for the future period https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License. (2069−2098. The increase of the radiative forcing leads to no significant increase in the drought duration, but instead causes a slight increase in the drought severity. In addition, the projected drought severity and duration show less variability for most 345 climate divisions in China. For example, the drought duration for Division 2 has the interquartile ranges (IQR) of 7 and 2 months for past and future climates, respectively. A remarkable exception is that there is no significant change in the drought duration and severity for Division 3 (i.e., the west central Tibetan Plateau). This does not indicate that the climate-induced drought risks will not increase since the multiple droughts with relatively long duration and high severity (indicated by the outliers) are projected over the west central Tibetan Plateau. 350 To further quantify the climate-induced change in drought risks, the return periods ("AND" and "OR" cases) of drought episodes based on drought duration and severity are assessed for the historical and future periods, as shown in Fig. 16. The 365 historical drought duration and severity were used to construct the parametric copula, which was then used to calculate the return period for each drought episode under past and future climates, leading to the box-and-whisker plots of return period in Figure 16. Due to the increase in the drought duration and severity as shown in Fig. 14, most climate divisions are expected to more than double the drought return period except Division 3. For example, the median return periods of historical drought episodes over Division 1 are 2.8 and 1.3 years for the "AND" and "OR" cases, respectively. And the corresponding cases of 370 return period are expected to increase by 132% and 183%, respectively, under RCP4.5 for the future period. We also observe that the projected return periods of extreme drought episodes increase from RCP4.5 to RCP8.5 for most climate divisions. For example, quite a few drought episodes with "OR" case return period higher than 10 years are projected for Division 3 under RCP8.5, and such extreme episodes do not occur under RCP4.5. Another typical example is the "AND" case for Division 10, where the occurrence of the drought return periods higher than 10 years significantly increases from RCP4.5 to RCP8.5. This 375 indicates that the extreme droughts are projected to increase in terms of frequency and the associated risks as the radiative forcing increases.

Conclusions
In this study, a probabilistic projection of multidimensional drought risks was developed through a copula-based two-stage Bayesian multi-model framework. An ensemble of five regional climate simulations was used to project future changes in 380 hydroclimatic regimes over China through an MCMC-based BMA approach. A Bayesian copula approach was also introduced to explicitly uncover potential interactions of the SPEI-detected drought characteristics and associated uncertainties, thereby improving the multidimensional drought risk assessment. The proposed Bayesian framework addresses the cascade of uncertainty in the climate-induced multidimensional drought risk projections. We also examined the performance of arithmetic ensemble mean (AEM) and BMA simulations in reproducing the historical climate, as well as Bayesian and frequentist copula 385 approaches used for multidimensional drought simulations.
The BMA climate simulation results indicate that the MCMC simulation provides not only plausible estimates but also the uncertainty information on the relative contribution of individual models in the multi-model climate simulation, improving the reliability of ensemble climate projections. And the PRECIS simulation has a relatively large contribution in the ensemble climate simulations. The AEM climate simulations based on the current CORDEX East Asia experiments show large biases 390 in most areas of China. In comparison, the BMA climate simulation can largely improve the simulation of precipitation and PET, with a higher correlation with observations and smaller biases than the AEM simulation. Both the AEM and BMA climate simulations can well reproduce the annual cycle of precipitation and PET in China, while AEM shows larger errors and BMA successfully corrects the errors. The introduced Bayesian copula approach not only provides equally plausible estimates compared to the frequentist copula approach but also explicitly uncovers the equifinality in the copula simulation. Such an 395 uncovered equifinality can improve the risk assessment of multidimensional droughts by providing multiple scenarios.
The MCMC-based BMA climate projections indicate a general increase in future precipitation and PET under RCP4.5 and RCP8.5 for most areas of China. A considerable decrease in the mean precipitation is also observed in northern Xinjiang, the western Tibetan Plateau, and Sichuan Basin. The increase in the radiative forcing leads to significant amplification of the PET increase but does not cause an obvious difference in the precipitation change. The projected changes in future precipitation 400 and PET cause a general increase in the drought duration and severity with decreasing variability. The drought risks are thus expected to significantly intensify, with more than doubling the multidimensional return period, even though the occurrence of drought episodes shows no significant increase and tends to decrease obviously. These findings reveal that China will experience more frequent extreme droughts, and the associated risks will be elevated due to the increase in the radiative forcing.
The proposed two-stage Bayesian multi-model framework enables a reliable projection of multidimensional drought risks, 405 which is vital for improving mitigation and preparedness of drought hazards and for enhancing sustainable water resources management. This framework can be directly applicable to assessing a variety of multidimensional extreme risks, which is useful for improving the robustness of the risk assessment of extreme events and natural hazards. It should be noted that although the MCMC-based BMA approach significantly improves the ensemble mean climate simulation, the potential errors are not completely corrected. It is thus desired to further improve regional climate simulations using the high-resolution 410 convection-permitting modeling systems in future studies. In addition, the time-invariant BMA weights determined by the historical data in multi-model climate projections may not well represent the nonstationary nature of climate dynamics.
Although the underlying uncertainty in the BMA weights was explicitly addressed in this study and previous studies also yielded plausible results (Olson et al., 2016(Olson et al., , 2018Shin et al., 2019;Terando et al., 2012), it is desired to develop nonstationary frameworks to further improve the credibility of climate projections.                 PET, leading to a probabilistic delineation of drought episodes. (i) represents the uncertainty information of the dependence 680 structure between drought characteristics (drought duration and severity) derived from the MCMC simulation, leading to the uncertainty in the return period of drought episodes as shown in (j). The red "whiskers" in (j) represent the uncertainty of drought characteristics resulting from the climate projection.
https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License.   https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License. Figure 12. Comparison of the empirical and fitted joint probability for 10 climate divisions in China. The fitted joint probability is separately calculated using copulas inferenced by Bayesian and frequentist approaches, as represented by the red and blue dots, respectively. https://doi.org/10.5194/hess-2020-247 Preprint. Discussion started: 20 July 2020 c Author(s) 2020. CC BY 4.0 License. Figure 16. The return periods of all drought episodes for the past and future climates over the 10 climate divisions. The setting of the box-and-whisker plot is the same as Figure 15. The return periods are calculated by the parametric copula constructed 760 for the historical drought duration and severity that are detected by the 6-month SPEI. The uncertainty in the return period results from the cascade of uncertainty as shown in Figure 1.