Climate change increased the compound extreme precipitation-flood events in a representative watershed of the Yangtze River Delta, China

A compound perspective on hydrological extreme events is of paramount significance as it may lead to damages with larger losses. In this study, an integrated framework, based on downscaled climate variables and hydrological model, i.e. the Soil and Water Assessment Tool, was applied to generate extreme precipitation (Rx1day) and extreme streamflow (Sx1day) series under historical and future climate conditions. Then the potential impacts of climate change for univariate and bivariate joint frequency of extreme precipitation and flood in Xitiaoxi River Basin (XRB), a representative watershed of the Yangtze River Delta, are detected. The compound risk of extreme precipitation and flood under different levels of joint return period for historical and projected periods is estimated by copula‐based two-dimensional approaches. The Rx1day and Sx1day under future scenarios changed by − 0.4% to 11.7% and 0.7% to 20.4%, respectively, compared to historical period based on univariate frequency analysis, indicating the increasing magnitude of the flood in the future. Climate change with different emission scenarios all have a driving effect on the rising coactivity of extreme precipitation and flood under compound flooding frequency analysis. In addition, the enhancement of climate change to extreme events is more apparent for extremes with higher return period and under the periods of 2080s. Moreover, the flood frequency designs are deduced by bivariate joint distribution are safer than that by univariate distribution. This study may provide actionable insights to formulate the planning scheme of flood control and disaster reduction under the changing environment.


Introduction
Global precipitation characteristics and flood situations are expected to undergo significant changes as a consequence of global climate change including increasing temperatures and varying precipitation patterns projected for the future climate (IPCC 2014). The incessant increase in magnitude and frequency of weather extremes (Yang et al. 2019;Blöschl et al. 2019;Vormoor et al. 2015), and uncertainties under climatic change scenarios may cause extremely severe disasters (Stott 2016;Chen et al. 2018;Swain et al. 2020;Hao and Singh 2020). For instance, extreme hydrological events are becoming one of the most serious natural disasters in densely populated regions, such as the Yangtze River Delta, the most developed region of China (Han et al. 2015;Jiang et al. 2020;Wang et al. 2016). Many studies have therefore explored the impacts of climate change on the occurrence of extreme events, and demonstrated that the climate was projected to be warmer and wetter with time and the magnitude of the changes in extreme hydrological events would experience vital changes (Dong et al. 2020;Niu et al. 2021;Tavakol et al. 2020;Zhu et al. 2021).
It is effective to utilize Global Circulation Model (GCM) projections to drive hydrological model, such as the Soil and Water Assessment Tool (SWAT), Variable Infiltration Capacity (VIC), MIKE-SHE and so on, to explore how climate conditions work on the hydrological process (Pastén-Zapata et al. 2020;Li and Fang 2021). Given that GCMs lack accuracy in spatial and temporal resolution, researchers use regionally more relevant scale climate models to downscale future climate scenarios (Chen et al. 2011;Fowler and Kilsby 2007). Particularly, the Statistical Downscaling Method (SDSM) model possesses the competence to apply climate information with greater physical significance and more accurate simulation in GCMs output to the statistical model (Wilby et al. 2002;Shen et al. 2018;Wilby and Dawson 2007). In hydrological and climatological studies, future climate scenarios are a reasonable description of the time and space distribution of the future climate state, which is based on certain driving forces and scientific assumptions. Representative Concentration Pathways (RCPs), as the most widely used future climate scenario, are advantageous for understanding the risks connected with different emission scenarios (Bhatta et al. 2019).
Over the past decades, many studies, highlighting the projected changes to weather and climate extremes, have been conducted, and most of them mainly focused on individual extremes or variables (Tofiq and Guven 2015;Meaurio et al. 2017;Bulti et al. 2020). Conventionally, each hydrological variable interacts intricately instead of existing independently. Considering the complex interactions exist between climate variables and ecosystem, a study by Jha et al. (2019a, b) utilize a probabilistic approach with copula framework to model the response of the dynamic ecosystem to extreme climate events, which needn't to rely on the assumption of stationarity or scale invariance. The statistical behaviors of high and low flow in two hydrological stations are thoroughly investigated by the joint probability function (Zhang et al. 2011a). Zscheischler and Seneviratne (2017) unraveled that the dependence structure between warm season temperature and precipitation lead to a much higher occurrence frequency of multivariate extremes.
The copula function (Skarl 1995), which without restriction on the specification of marginal distributions and their uniformity, is conducive to deduce the joint probability distribution of random variables (Favre et al. 2004;Joe 2014). A number of univariate and copula-based multivariate approaches were developed for co-occurrence frequency analysis of extremes, such as drought (Ayantobo et al. 2019;Sun et al. 2019;Xu et al. 2020), and flooding hazard (Muñoz et al. 2020;Balistrocchi et al. 2017;Yu et al. 2019;Yin et al. 2018;Liu et al. 2020;Vinnarasi and Dhanya 2019), thereby providing a joint frequency perspective of the associative structure between the depicted factors. For instance, Sun et al (2019) utilized a multidimensional copula model to combine the two kinds of disaster in terms of investigating the potential two future scenarios impacts on the simultaneous occurrence and joint return period of droughts and hot extremes in the Loess Plateau of China, and found that the compounding occurrence of drought with long-term hot extremes will be more severe and frequent. Muñoz et al (2020) evaluated compounding effects of terrestrial and coastal flood drivers and wetland elevation accuracy on maximum floodwater height and velocity with a bivariate statistical analysis framework. Yin et al (2018) explored the variations of flood peak and volume in Ganjiang River basin under different climate scenarios by fitting univariate and copulabased bivariate distribution, which indicated that the impacts of climate change on the future bivariate flood quantiles are considerable.
Centered on the Yangtze River Delta region, the joint frequency analysis of the multivariate characteristics based on the coupled extreme precipitation-streamflow under changing climate is rarely inquired about in literature. Furthermore, flood is one of the most destructive and widespread natural disasters over China, and climate change has anticipated to exacerbate the frequency of extreme flood events in the future. Herein, this research aims to detect the joint responses of extreme precipitation and floods under the climate changing environment by establishing suitable marginal and two-dimensional Copula function of the extreme precipitation and streamflow, which expands on our previous work . The results would improve the acknowledge of the disaster-causing process in watershed flood; thereafter, actionable insights are provided for decision makers to formulate the planning scheme of flood control and disaster reduction.
Clearly, the overarching structure of the research is organized as follows. In Sect. 2, the study area and data involved in this paper are displayed. Section 3 is devoted to the methodologies used in the paper, including the prediction technique of extreme events and the copula method. In Sect. 4, the best fitted marginal and joint distributions are established through Copula functions, in order to evaluate the impact of climate change on the frequency of univariate and bivariate joint design events. At last, conclusions are presented.

Study area
The Xitiaoxi River Basin (XRB), located between 119°14 0 E-119°45 0 E and 30°23 0 N-31°11 0 N within the Yangtze River Delta, is one of the major tributaries of the upper reaches of Taihu Lake (Fig. 1). The south and northwest part of XRB distributes mountains, and flat alluvial plain lies in the northeastern part (Zhang et al. 2011b;Zhao et al. 2011). The basin dominates by subtropical monsoon climate with an average annual temperature of 15.5°C and a mean annual rainfall of 1584 mm. Specially, the flood season (from April to October) contains over 70% of the rainfall. As one of the three largest water volume supplies (approximately 27.7%) of the Taihu Lake, the XRB is acknowledged to be an area experienced rapid development of urbanization and turbulent climate change that may affect the occurrence mechanism of local extreme hydrological events. Thus, based on the representative region of the Yangtze River Delta-the XRB, this research aimed to demonstrate the inter-response between extreme hydroclimatic events under future climate scenarios.

Datasets
The datasets used mainly includes the reanalysis data desired by the SDSM model, and the data required for driving the SWAT model. Detailed presentations of data usage are shown in Table 1. The National Centers for Environment Prediction (NCEP) reanalysis data, GCMs future scenario output factors and historical meteorological data were applied to conduct statistically downscaling. The NCEP reanalysis data comprised the daily series of 26 atmospheric circulation factors for the time period 1961 to 2015 with a 2.5°9 2.5°horizontal resolution. Three GCMs, namely the BCC-CSM 1.1 (Beijing Climate Center, Climate System Model, version1.1), the NorESM1-M (Norwegian Earth System Model, version 1), and the CanESM2 (Second Generation Canadian Earth System Model), under RCP2.6, RCP4.5 and RCP8.5 emission scenarios, were utilized for exporting the historical and future climate time series data. All GCM datasets were uniformed to the same resolution of NCEP to avoid any biases caused by differences in scale. The observed weather data included daily precipitation (mm), minimum and maximum temperature (°C) covering the period between 1961 and 2015. Two representative periods of climate projection scenarios  are selected, i.e., the mid-century scenario (2036-2065, the 2050s) and the late-century scenario (2070-2099, 2080s), compared with the baseline period (1961-1990, 1970s) of the three GCMs.
To minimize climate change projections uncertainty caused by the differences of parameterization schemes in various GCMs, the ensemble means of the output in the three GCMs were then implemented in this study. Geographic data and observation data are required to set up the hydrological SWAT model. A spatial database including the digital elevation model (DEM), land use and land cover, and soil map datasets. Observation data included hydrological and meteorological data. The observed daily streamflow data during 1994-2015 was obtained from the Laoshikan Reservoir, Fushi Reservoir, and Hengtangcun Station. Observed meteorological data contained daily precipitation data, taken from 1972 to 2015, were collected from 16 rainfall stations, and five climate variables, namely, daily maximum and minimum temperature, mean wind speed, relative humidity, and solar radiation measurements were provided as input with daily interval for the period from 1961 to 2015. Hydrometeorological data taken from 1994 to 2015 were used for the construction, calibration and validation of the hydrological model.

Methods
For this study, a hybrid method, coupling hydroclimatic models with copula functions, was applied to investigate the univariate and multivariate variation characteristics of hydrological extreme event under future scenarios. Using a statistical downscaling model (SDSM), daily precipitation and other climate variables for the historical and future periods were projected from the outputs of three GCMs under three climate change scenarios. The hydrological processes of XRB were simulated by SWAT model to export streamflow during 1961-2099. Copula models were utilized to investigate the impacts of climate change on the univariate and multivariate extreme precipitation and flood events. A schematic flowchart of the methodology applied in this research is given in Fig. 2.

Hydrological model and extreme indices
The SWAT model, developed by United States Department of Agriculture's Agriculture Research Service (Arnold et al. 1998), can effectively simulate hydrological processes such as runoff and evapotranspiration for a watershed. SWAT is popular in assessing runoff responses in watersheds around the world for its strong portability, suitability for scenario analysis, and friendly interface. The SDSM (Wilby et al. 2002) is a combination of conditional weather generator and multiple regression analysis that can be used to evaluate and predict time-varying meteorological parameters describing local variables. Hence, SDSM establishes mathematical relationship and downscaling the climatic data between regional predictors from GCMs and local predictands by regression and stochastic approaches.
Our previous study  showed that the SWAT and SDSM model exhibits sufficiently capable of simulating the hydrological processes and estimating regional climate scenarios for the XRB. Through the integrated application of SWAT and SDSM model,   and future scenarios (2006-2099) The National Centers for Environment Prediction (NCEP); the website of the Earth System Grid Federation (https:// esgf-node.llnl.gov/) extreme precipitation and streamflow data under future scenarios are simulated and predicted to detect the variation characteristics of the joint frequency. In this study, two extreme flood event indicators are selected: (1) a precipitation index (Annual maximum 1-day precipitation, Rx1day, mm); (2) a river discharge index (Annual highest daily flow, Sx1day, m 3 /s). They are extracted from data series of precipitation and streamflow in different periods under diverse climate scenarios in XRB by Annual maximum sampling method.

Copulas theory
It is important to note that hydrological variables are interrelated, hence, univariate analysis is not competent to provide integrated research of flood ). The Copula function, initially raised by Sklar in 1950s, is an expression that can construct the multidimensional joint distribution of two or more independent variables by marginal distributions and correlation framework (Sklar 1959;Nelsen 2006;Li et al. 2018;Salvadori and De Michele 2007). The copulas have been broadly used to explore the correlation structure between variables regardless of their marginal distribution types in hydrological extremes study. According to the Sklar's theorem, the two-dimensional copula function can be expressed as follows: in which F (a, b) is joint distribution. The function C is a two-dimensional copula in the form [0,1] 9 [0,1] ? [0,1] of the pair (a, b). Where F 1 (a) and F 2 (b), denoted by u 1 and u 2 respectively, represent marginal cumulative distribution functions of random variables a and b. And h indicates parameter of the time-dependent copula. Most researchers concentrate on a handful of copulas (i.e., Frank, Clayton, Gaussian and Gumbel) resulting a negligence of exploring more suitable copulas. These are desirable for hydrological extremes analyses owing to the fact that they can be easily constructed and can describe either positive or negative dependence (Zhang et al. 2011a;Zhang and Singh, 2007). Copula families differ in describing dependence structure and the parameter uncertainty of variables in details. More copulas that have not been completely explored in the hydrological and climatological literature should be considered. The copula functions finally selected in this study are the Clayton copula type, the Joe type, the Tawn type and the Burr type, which are more flexible in modeling various dependence structures between random variables. The Clayton and Joe copula is advantageous when modeling datasets with inconsistent dependence (Jha et al. 2019a,b). The Tawn and the Burr family are able to describe asymmetric skewed dependence structures exhibited in variables (Sadegh et al. 2017).

Establishment of marginal distribution function
Parameters of seventeen probability distribution functions (i.e., Beta, Birnbaum-Saunders, Exponential, Extreme value, Gamma, Generalized extreme value, Generalized Pareto, Inverse Gaussian, Logistic, Log-logistic, Lognormal, Nakagami, Normal, Rayleigh, Rician, t location-scale, and Weibull distributions) are calculated by a maximum likelihood algorithm that minimizes the distance between observations and model simulations. Then the BIC (Bayesian Information Criterion) metric is used to select the most suitable marginal distribution. What's more, the chis-squared (v 2 ) goodness-of-fit test is applied to statistically describe whether or not a probability distribution fits the observed data. The Quantile-Quantile (Q-Q) plotting aims to compare the fitting effect of the function and empirical probability values.
To determine whether there was an inter-dependency between the extreme precipitation and extreme streamflow for different correlation coefficients namely, Kendall rank, Spearman's rank-order, and Pearson product-moment correlation coefficients, are used in this paper to assess the correlation relationship between the two variables.

Selection of fit copula function
Quantitatively identifying the reliability of copula function is required to effectively construct the structure of dependence between extreme precipitation and streamflow. The optimal copula which we have used in this investigation were identified using five model performance indices, the root mean square error (RMSE) (Willmott et al. 1985), the Nash-Sutcliffe efficiency (NSE) (Zhang and Singh 2007), Akaike's information criterion (AIC) (Akaike 1987), the maximum likelihood estimation (MLE) (Sorooshian and Dracap 1980) and BIC (Schwarz 1978). Their equations are formulated as follows: where N is the number of observations; i is the serial number of the observation; x e (i) represents the multivariate empirical probability of observations and x o (i) represents the theoretical probability using fitted probability distribution; k is the unknown number of parameters of probability functions. The RMSE varies in the interval 0 to ? !, and a lower RMSE value associates with a more perfect model fit. In general, a copula model is reliable for hydrological extreme event when the NSE is larger than 0.75, in which the range of changes is from -1 to 1 (Motovilov et al. 1999). Compared with MLE, where a larger value is preferred, the lower the coefficients of AIC and BIC (negative), the better the copula model.

Univariate and bivariate return period
Once the marginal distributions of Rx1day and Sx1day are defined, the univariate return periods of extreme risk are also derived under the historical and future scenarios. The relevant calculation formulas can be expressed as follows: where T denotes the return periods for a Rx1day or Sx1day greater than or equal to certain values of x. F means the cumulative probability distribution functions of Rx1day or Sx1day. To obtain the joint probabilistic characteristics between the two dependent variables, the copula function of joint return periods is determined: where T joint represent the joint return period when R C r or D C d.
The bivariate statistical analysis under uncertain is assessed based on the Multivariate Copula Analysis Toolbox (MvCAT) developed by Sadegh et al. (2017Sadegh et al. ( , 2018, which can calculate the marginal distribution for each extreme indicator and select the most suitable copula function for pair-wise variables. Firstly, to find an appropriate model that optimally fits the available data, the parameters of the 17 different marginal probability distribution functions were estimated in MvCAT. And then, a total of 26 copula distribution families for bivariate modeling were selected to develop the one that best describes the joint dependence structure of extreme precipitation and extreme streamflow based on the univariate marginal function. The mathematical expression of these copulas is available in Table 1 of Sadegh et al. (2017).

Marginal and joint probability distribution of Rx1day and Sx1day
The Kendall, Spearman and Pearson correlation coefficients were figured to evaluate the dependence between (Rx1day, Sx1day) pairs under historical and future climate scenarios, and their caculated values are listed in Table 2. The coefficients exhibit a statistically significant relation at the 0.05 significance level, indicating that the factors under each scenario present strong positive correlation. The results also show that the future correlation between Rx1day and Sx1day under low emission scenario (RCP2.6) increased firstly and then decreased compared to the history. While in RCP4.5 and RCP 8.5 scenarios they exhibit an opposite trend, that is decreased firstly and then increased. This implies a probability that the impact of Rx1day and Sx1day can be altered by emission scenarios.
Hence, the two variables can be used to construct the joint distribution function through the Copula function. Seventeen univariate distributions, out of which the most suitable one is chosen based on the BIC metric, were selected as candidates to fit the simulative and predicted data in terms of Rx1day and Sx1day from the XRB. The parameters and goodness-of-fit of each marginal distribution was estimated using Chi-square tests. All of the fitted univariate distributions passed the Chi-square test, that is, at 5% level. A summary of best marginal distributions fitted to Rx1day and Sx1day is presented in Table 2. According to the optimal marginal distribution function mentioned in the above analysis, empirical and theoretical data for each variable is established. Figures 3 and 4 illustrate the corresponding relations between the fitted and empirical probability of Rx1day and Sx1day. Q-Q plots visually present goodness-of-fit of fitted data to empirical data, showing that they are suitable for building joint dependence. The joint probability distributions between Rx1day and Sx1day under different scenarios were examined using RMSE, NSE, MLE, AIC and BIC. The evaluation metrics suggest that these corresponding best models exhibit an acceptable level of fit, indicating that the selected Copula Functions should be used in this investigation (Table 3).

Impacts of climate change on individual
extreme hydrological events using univariate frequency analysis Figure 5 demonstrates the marginal return period for Rx1day and Sx1day for the history and future periods under RCP 2.6, RCP 4.5 and RCP 8.5, respectively, with the corresponding values of Rx1day and Sx1day for different return periods (T = 2a, 5a, 10a, 25a, 50a, 100a) as illustrated in Table 4. The results showed that the corresponding values of Rx1day and Sx1day for a certain return period in future climate conditions (i.e., 2050s and 2080s) are higher than those in historic period, i.e., 1970s, especially under the RCP 8.5 scenarios. In addition, climate change under the low emission scenario (RCP 2.6) have relative greater impacts on precipitation events with lower return periods, such as 2a, 5a, and 10a, during the periods of 2050s, while the impacts of climate changes under higher emission scenarios (RCP 4.5 and RCP 8.5) on rainfall events with higher return periods, such as 25a, 50a, and 100a, are greater at the same return periods (Fig. 5a and Table 4). With the respects of the periods of 2080s, climate change has relative greater impacts on rainfall events with high return periods, such as 25a, 50a, and 100a, except for RCP4.5 scenarios ( Fig. 5b and Table 4). Similarly, climate change impacts on Sx1day with high return   Table 4). It can be seen that the changes of Rx1day and Sx1day in future period, compared to those values in historical scenarios, increase with the return period, and Sx1day show a larger increment than Rx1day. The Rx1day and Sx1day under future scenarios increased by -0.4 to 11.7% and 0.7 to 20.4%, respectively, compared to historical periods. Above results imply that the changing environment has severer impacts on Sx1day than Rx1day, especially under the RCP 8.5 scenarios. This is consistent with the findings in the literatures that the future projected streamflow appears to increase in general and shows fairly dramatic growths under RCP 8.5 scenario (Nilawar and Waikar 2019; Zhang et al. 2016). These findings suggest that the XRB may suffer higher flood risk in the future. In addition, the impacts of climate change on hydrological events differ under different emission scenarios, due to different warming effects caused by different emission scenarios. High emission scenarios (such as RCP 4.5 and RCP 8.5) would lead to an increase in the frequency and intensity of extreme rainfall, and resulting in a corresponding flood response. There has a greater impacts of climate change on small flood under low emission scenario (RCP 2.6), while on large flood under higher emission scenario (RCP 4.5 and RCP 8.5). High emission scenarios would lead to intensified climate warming and more frequent occurrence of extreme rainfall, which leads to more remarkable response of corresponding higher magnitude floods. We have added the discussions about the results in the pape. Since precipitation is recognized to be the direct influence factor of the flood, these changes in streamflow were closely related to simultaneous changes in precipitation (Li and Fang 2021). The preliminary correlation analysis of the data also suggested that extreme streamflow is significantly and highly influenced by the extreme precipitation. Variation of floods will also magnify the impact of climate change on precipitation in reverse, therefore, further analysis is needed through joint frequency.

Impacts of climate change on compound extreme precipitation-flood events using bivariate frequency analysis
The curves for joint return period of Rx1day and Sx1day under historical and future scenarios are demonstrated in Fig. 6. The colored lines represent the combined probability of extreme precipitation and extreme flood events, and the possibility of simultaneous occurrence of the two events is reflected in the value of probability densities. Acorrding to the summery report from MvCAT, except for a few distribution types that are obviously not suitable for the fitting of extremum sequence, most distribution types can fit the low-frequency events well. The difference of these distribution types is mainly reflected in the fitness of tail data, which also shows the importance of marginal distribution optimization in developing multi-dimensional joint distribution from the side. The most suitable tail fit is chosen from seventeen univariate distributions based on the BIC metric. Finally, the optimal copula functions of Rx1day and Sx1day in the 1970s, in the 2050s under the RCP 2.6 scenario, in the 2080s under the RCP 2.6 scenario, in the 2050s under the RCP 4.5 scenario, in the 2080s under the RCP 4.5 scenario, in the 2050s under the RCP 8.5 scenario, in the 2080s under the RCP 8.5 scenario are Tawn, Burr, Tawn, Clayton, Joe, Burr and Joe respectively. For the sake of explanation, the interval where the combined probability of Rx1day and Sx1day is greater than 50% is expressed as the combination interval with high probability. For the historical and RCP 2.6 scenarios, the intervals of Rx1day are narrow and which of Sx1day are wider under each return period, that is, Sx1day is greatly affected by Rx1day in these situations (Fig. 6a-c). In this sense, the probability that the XRB suffer from the extreme precipitation and floods simultaneously is very high. In the time periods of 2050s under the RCP 4.5 scenario, it is evident that the corresponding relationship between the Rx1day and Sx1day with joint return periods higher than 5a is not obvious (Fig. 6d). As to the periods of 2080s, the combination interval with a high probability of Rx1day and Sx1day is relatively concentrated, which indicates that the variation range of Rx1day, as well as Sx1day, is slight at this stage (Fig. 6e). Figure 6f reflects that the higher the return periods are, the larger the values of Rx1day are in high probability combination interval in the 2050s under the RCP 8.5 scenario. Whereas the variation range of Sx1day does not change much with the regular change of return period, which revealed that the influence of Rx1day Values shown in parenthesis () are variations (%) of future scenarios compared to historical period on Sx1day gradually weakens with the increase of the joint return period. During the 2080s, the fitted results represent that Rx1day and Sx1day have a good correlation because of the little variation ranges of Rx1day and Sx1day with return period. Table 5 shows the variation range of Rx1day and Sx1day in high probability combination interval. Here we mainly focus on events under historical, RCP 2.6 and RCP 8.5 scenarios, since the changing trend of the two variables is not obvious under the RCP 4.5 scenario. From the perspective of different return periods, the joint return periods for compound events higher than 5a in the listed combination interval under the two future scenarios are visibly greater than the those in historical scenario, particularly for the RCP 8.5 scenario. And the relevant Rx1day and Sx1day values cover a broader range, for example, under RCP 8.5 scenario, the 50a combination interval covering a range with marginal Rx1day values in the 1970s, 2050s, and 2080s ranging from 31.7 to 31.9 mm, 34.2 mm to 35.1 mm, and 35.2 mm to 36.6 mm, respectively. This suggests that smaller combinations of Rx1day and Sx1day that can induce a 50a event in the 1970s are less likely to prompt such dangerous events in future scenarios and the flood risk is comparatively higher in the future. It can be qualitatively judged that emission scenarios have a driving effect on the rising coactivity of bivariate combination of Rx1day and Sx1day.
In the 2050s, the climate change during RCP 8.5 scenario in larger return periods have a more significant influence on the rise of combination values than in small return periods, compared to the RCP 2.6 scenario. The design values of Sx1day in combination interval for 2a joint return period under the RCP 2.6 and RCP 8.5 scenarios range from 158.97 to 176.98 m 3 /s and from 156.55 to 168.45 m 3 /s, respectively. While for 100a joint return period, the values range from 218.63 to 228.25 m 3 /s under the 2.6 scenario and from 224.27 to 235.27 m 3 /s under the 8.5 scenario. However, the climate under the RCP 8.5 scenario always keeps a stronger driving effect on the increase of Rx1day and Sx1day under each joint return period, and the corresponding values of Rx1day and Sx1day in combination interval are also covering a wider range in the 2080s. This indicates that future extreme flood events are more uncertain and more difficult to predict. Climate change, characterized by global warming and precipitation extremes, as well as anthropogenic activities including rapid urban growth, are involved in the furious turbulence of hydroclimatic extremes (Clark et al. 2016;Nathan et al. 2019). Variation making in persistence, predictability, and probability of atmospheric patterns and the selection of copula function would exert an influence on the uncertainty range of extreme events (Faranda et al. 2020;Stottl 2016;Liu et al. 2020). Overall, the driving effect of high emission scenario on extreme precipitation and floods is more obvious for events with higher return period, especially in the period of 2080s.

Comparison of univariate and bivariate joint analysis
The extreme precipitation and flood under certain return periods generated by bivariate joint analysis during historical and future periods are showed in Fig. 7. Both the Rx1day and Sx1day show an upward trend, and the design values of bivariate joint analysis are larger than that of univariate analysis for most return periods. It can be noticed from Fig. 7 that the change degree of Rx1day and Sx1day during the historical period and future scenarios show little difference for larger return periods. For example, the design values of Rx1day under 2a return period increases from 28.8 mm in univariate analysis to 28.5 mm in bivariate joint analysis by 2.4% for the RCP2.6 in the 2050s (Fig. 7b). While the values under 100a return period are 33.7 mm in univariate analysis and 33.9 mm in bivariate joint analysis (increased by 0.6%). The changes of other design values of Rx1day and Sx1day under two climate scenarios during 2050s and 2080s are similar. These results imply that the joint combination contribute to a safer design criteria, and the flood frequency design deduced by bivariate joint distribution is securer than that by univariate distribution. These are consistent with the results of Lin et al. (2021) in Ganjiang River basin and Duan et al. (2016) in Huai River Basin, indicating that the results are reasonable and universal. Considering the security of water engineering, flood control design based on the joint distribution of two variables can reduce the risk of flood disaster (Lin et al. 2021).

Summary and conclusions
Flood is a main natural disaster occurring on the Yangtze River Delta, and climate change has anticipated to exacerbate extreme flood events. Thus it is important to take the change into account when exploring the hydrological response for future flood characteristics. To identify the univariate and bivariate joint frequency of hydrological extreme events at different climate scenarios during historical and future period, we applied an integrated framework with the combination of the climatic simulation under three GCMs, streamflow simulation and projection based on the SWAT model, and univariate and bivariate frequency analysis according to marginal distribution functions and copula-based models. Compared with the historical period, the impacts of climate changes on extreme precipitation and streamflow in the XRB were quantified. The main conclusions are summarized as follows: Based on univariate frequency analysis, the corresponding threshold of Rx1day and Sx1day in future period are larger than those in the historical period under the same return period. Compared to historical period, the Rx1day and Sx1day under future scenarios changed by -0.4% to 11.7% and 0.7% to 20.4%, respectively. Most of them exhibit an increasing trend, indicating that the magnitude of flood would be higher under future climate conditions. During bivariate analysis the significant nonlinear correlation between climate variables were captured. There is a considerable increase of Sx1day with the increase of Rx1day, and flood is greatly affected by precipitation in these situations during the 1970s under RCP 2.6 scenario. However, the effect of extreme precipitation to flood gradually weakens with the increase of joint return period in the 2050s under the RCP 8.5 scenario. Rx1day and Sx1day with high probability combination interval in the future are dramatically larger than those of historical scenario, together with a broader variations range, especially during 2080s under RCP 8.5 scenario. This demonstrates that future extreme flood events are more uncertain and more difficult to predict. The effects of climate change on the enhancement of Rx1day and Sx1day are more obvious for higher return period events, especially during 2080s. By the comparation of univariate and bivariate frequency analysis, the results reveal that the flood frequency design deduced by bivariate joint distribution is more safer than that by univariate distribution. Affected by the changing environment, the rainstorm and flood disasters in XRB are likely to occur more frequently, which indicating the urgent needs for the safer flood frequency design.
In general, the intensity and frequency of extreme hydrological events are increasing under the changeful climatic environment, and this study also indicates that the different enhancements extreme events under different emission scenario of carbon dioxide, which should cause the attention of relevant departments. It should be noticed that this study only focuses on the impact of different climate scenarios on regional hydrological processes. However, the influence of climate change on floods is a comprehensive result; the responsiveness of the different driving factors as well as the influence degree on the frequency of changes in regional hydrological processes, should be quantitatively analyzed in further research.