Observation‐Constrained Projection of Flood Risks and Socioeconomic Exposure in China

As the planet warms, the atmosphere's water vapor holding capacity rises, leading to more intense precipitation extremes. River floods with high peak discharge or long duration can increase the likelihood of infrastructure failure and enhance ecosystem vulnerability. However, changes in the peak and duration of floods and corresponding socioeconomic exposure under climate change are still poorly understood. This study employs a bivariate framework to quantify changes in flood risks and their socioeconomic impacts in China between the past (1985–2014) and future (2071–2100) in 204 catchments. Future daily river streamflow is projected by using a cascade modeling chain based on the outputs of five bias‐corrected global climate models (GCMs) under three shared socioeconomic CMIP6 pathways (SSP1‐26, SSP3‐70, and SSP5‐85), a machine learning model and four hydrological models. We also utilize the copula function to build the joint distribution of flood peak and duration, and calculate the joint return periods of the bivariate flood hazard. Finally, the exposure of population and regional gross domestic product to floods are investigated at the national scale. Our results indicate that flood peak and duration are likely to increase in the majority of catchments by 25%–100% by the late 21st century depending on the shared socioeconomic pathway. China is projected to experience a significant increase in bivariate flood risks even under the lowest emission pathway, with 24.0 million dollars/km2 and 608 people/km2 exposed under a moderate emissions scenario (SSP3‐70). These findings have direct implications for hazard mitigation and climate adaptation policies in China.

HMs. More recently, post-processing techniques such as machine learning (ML) have emerged as a popular alternative method for enhancing the performance of hydrological simulation . ML post-processing contributes to the reduction of systematic errors and improves the accuracy of climate or hydrological simulations and projections (Baran et al., 2019;Gong et al., 2022;Hemri et al., 2015). ML can enhance the modeling performance of physical processes (Bonavita et al., 2021;Ma et al., 2020), enabling researchers to acquire new insights on climate change impacts. In the last two decades, support vector machines (SVM; Lin et al., 2009), artificial neural networks (ANN; M. Cheng et al., 2020) and Long Short-Term Memory (LSTM; Lee et al., 2020;Lees et al., 2022) have been utilized to strengthen hydrological predictions by using large-sample (multi-site) data. LSTMs in particular can achieve promising results, provide insights into hydrological processes (Lees et al., 2022), and reduce uncertainty in projections of future flood dynamics (D. Y. Li et al., 2021;W. B. Liu et al., 2021).
As floods are multifaceted hazards in terms of their timing, peak discharge, volume, duration, and hydrograph shape (e.g., Xu et al., 2019), deciding how to realistically represent flood characteristics is an important issue. Over the past decade, numerous studies have employed climate models to evaluate climate change impact on flood dynamics, but most studies have used a univariate hydrological analysis framework (e.g., Kim et al., 2019). These univariate frameworks can only provide limited assessment of the flood occurrence likelihood. As flood peak and duration are both important characteristics for flood risk assessments and water resources management in China, this study focuses on evaluating the joint occurrence of flood peak and duration-referred to as "bivariate floods" hereafter-by using copulas within a bivariate framework. Yin et al. (2020) found that flood peak and volume can show different patterns of variation in a warming climate, and suggested that flood risks should be explored in a multivariate framework. It is well established that river floods with either large peak discharge or a long persistence can result in severe effects, including infrastructure failure and environmental vulnerability (Barth et al., 2017;Brunner et al., 2019;Gaal et al., 2012). Some studies have pointed out that the global climate crisis is linked to socio-economic factors (D. Jiang et al., 2022), such as gross domestic product (GDP) and population. As human adaptability to climate change also depends on increases in income per capita (Leimbach et al., 2017;Yin, Gentine, et al., 2023;Yin, Guo, et al., 2023), environmental degradation can increase exposure or vulnerability to climate hazards. Flood events are among one of the costliest climate hazards (UNDRR, 2020;UNISDR, 2015), thus it is necessary to assess their impacts on socioeconomic (e.g., GDP) exposure under climate change. As floods may cause widespread societal damages, research is increasingly considering socioeconomic exposure to floods. However, most previous evaluations (e.g., W. Liu & Sun, 2019;Park et al., 2018;Peters, 2016) have assumed a fixed socioeconomic state, likely due to data constraints, while the population and assets exposed to flooding are actually expected to change dynamically over time in the real world. Without considering such dynamic changes, the results might be misleading for climate policy planning. Recently, T. Jiang et al. (2017Jiang et al. ( , 2018 released population and GDP data under different shared socioeconomic pathways (SSPs), which consider the "universal two-child policy" in China and realistically capture China's socioeconomic conditions. However, few studies have assessed the impacts of flooding on socioeconomic systems by using this SSP data set under both a univariate and bivariate framework.
Here, we investigate how bivariate flood risks (i.e., flood peak and duration) and socioeconomic exposure might change in China, a flood-prone country, under three CMIP6 SSPs (i.e., . To explore the physical mechanisms behind runoff generation, we first construct a random forest model to evaluate the roles of meteorological variables driving runoff. Then, we drive the HM-LSTM hybrid model with the outputs of five bias-corrected CMIP6 GCMs under three SSPs to project future daily river streamflow over 204 catchments in China. Subsequently, the joint distributions of flood peak and duration are established by using copulas under both past and future climates. The shifts of bivariate flood characteristics are evaluated by using the most likely realization, and the population and GDP exposed to increased flood risks are quantified. Finally, we assess the uncertainty of flood projections sourced from different procedures.

In Situ Observation Data Set
A high-resolution (0.5° × 0.5°) gridded meteorological data set in China from 1961 to 2016 is employed, which includes daily temperature (maximum, minimum and average, °C), and daily precipitation (mm) from the National Meteorological Information Centre. This data set is considered the most up-to-date gridded meteorological data 10.1029/2022EF003308 3 of 21 in China and has been extensively employed in hydrological research (Wu et al., 2018;. To represent flood dynamics, we also collect daily streamflow data from 463 hydrological stations, covering different periods from 1961 to 2016. After strict data quality control and screening for robustness of the analysis, we retain 204 catchments covering eight major river basins in China. Further details about the catchment screening process are in Text S1 of Supporting Information S1.

ERA5-Land Data
ERA5-Land is a replay of the ERA5 climate reanalysis' land component, which generates a total of 50 variables characterizing the water and energy cycles over global land at a high resolution (0.1° × 0.1°) from the European Centre for Medium-Range Weather Forecasts (ECMWF, Malardel et al., 2016). We employ hourly data of six climate variables: precipitation (pr), air pressure (ps), 2 m temperature (T 2m ), dew-point temperature (T dew ), surface downwelling shortwave radiation (rsds), and surface downwelling longwave radiation (rlds). All the hourly variables are translated into a daily scale before implementing data analysis.

GCM Outputs and Socioeconomic Data
For the future climate scenarios, we use a multi-model ensemble consisting of five GCMs listed in Table S1 of Supporting Information S1 under three SSPs (i.e., SSP1-26, SSP3-70, and SSP5-85) from the latest CMIP6. Numerous studies have found that the raw outputs of precipitation and temperature from CMIP6 climate models are overestimated in Asia and have non-negligible uncertainties (Chai et al., 2022). To reduce the systematic biases of climate models, we use the bias-corrected daily outputs from the Intersectoral Impact Model Intercomparison Project 3b (ISIMIP3b). Our selected variables include pr, ps, temperature (daily-average, maximum and minimum), rsds, rlds, relative humidity (RH), and specific humidity (SH). The ISIMIP3b data has been systematically bias-corrected and downscaled to 0.5° × 0.5° spatial resolution by using climate observations from 1850 to 2100. In the bias correction process of the ISIMIP3b, Lange (2019) employed the parametric trend-preserving method and considered the interdependence between different variables, thus showing substantial advantages over its predecessor (i.e., the ISIMIP2). Many studies have employed this data set and confirmed its robustness in assessing climate change (H. Yin, Gentine, et al., 2023;Yin, Guo, et al., 2023). We choose 1985-2014 as the historical period and choose 2071-2100 to depict the end-of-century scenario, taking into account the long-term effects of climate change.
To assess the potential socioeconomic exposure to floods, we use the population and GDP data from the SSPs to represent various climate change mitigation and adaptation options. The population and GDP data from the Intergovernmental Panel on Climate Change's (IPCC's) three selected SSPs are combined with the corresponding greenhouse gas emission scenarios, providing a total of three matrix frameworks (i.e., . In this study, we use the open-access data set of "the universal two-child policy" prediction released by T. Jiang et al. (2017Jiang et al. ( , 2018, which considers the results of China's successive population and economic censuses as well as year-by-year statistical yearbooks. The model for producing this data set has different parameters (e.g., fertility, mortality, migration rate, education level) under different SSPs and considers parameter modification according to China's population and policies. The data set also provides an extensive survey of the fertility intentions ("expected number of children") of Chinese citizens, which has been reported as highly reliable (Xing et al., 2022). The economic data is projected by using the Cobb-Douglas model and the Population-Environment-Development (PED) model. The Douglas model and the PED Analysis model are used to estimate China's socioeconomic index from 2010 to 2100 (T. Jiang et al., 2017). This data set has been utilized extensively for evaluating the socioeconomic effects of extreme hydrologic events Su et al., 2018;.

Methodology
The workflow of this study is illustrated in Figure 1, which contains six modules. First, we perform a sensitivity analysis of the in situ streamflow using meteorological factors from the ERA5-Land data set. Then, the HM-LSTM model is built by combining physical model and data-driven model, with consideration of lag time. The ISIMIP3b output is used to drive the HM-LSTM model for projecting future streamflow scenarios, and then the flood peak and duration are extracted. Under the bivariate framework, the most likely realization method is 4 of 21 employed to estimate changes in flood risks under future climates. Subsequently, we systematically evaluate the exposure of future GDP and population to bivariate floods. Finally, we clarify the physical mechanisms and apply the analysis of variance (ANOVA) method to complete the uncertainty analysis.

Deriving Near-Surface Specific Humidity and Relative Humidity
The C-C relationship describes the relationship between saturated water vapor pressure e sat and air temperature T (Koutsoyiannis, 2012): where T 0 and e s0 are the integration constants, taken as 273.16 K and 611 Pa, respectively; L v denotes the latent heat of vaporization, taken as 2.5 × 10 6 J kg −1 ; R v denotes the water vapor gas constant, taken as 461 J kg −1 K −1 .
T dew characterizes the temperature at which the air cools to saturation condition under constant water vapor and pressure. Substituting T 2m and T dew from the ERA5-Land data into Equation 1, we can deduce the near-ground RH: RH = e sat (T dew )/e sat (T 2m ).
SH is the ratio of the mass of water vapor to the total mass of the air mass, which is deduced by using ERA5-Land near-ground ps and T dew from the following equation (Simmons et al., 1999):

Sensitivity Analysis of Runoff to Meteorological Variables
The random forest is an ensemble learning method for classification or regression that works by constructing a large number of decision trees during the training. In our study, the random forest model is first constructed using meteorological variables (i.e., pr, ps, SH, RH, rlds, rsds, T 2m ), and then used to quantify each meteorological variable's contribution to streamflow. Assuming that perturbing one variable does not change the value of the other variables, for each of the 100 trained trees in the random forest (Genuer et al., 2010), we perturb one of the meteorological variables by 1 standard deviation (SD), and streamflow is again predicted using the same random forest model with the same meteorological variables, but including the perturbed variables. This method has been widely used to quantify the contribution of a factor (e.g., pr, ps, SH, RH, rlds, rsds, T 2m ) to the target variable (Green et al., 2020). The percent of variability in streamflow is calculated as follows: where S M denotes the sensitivity of streamflow by one of the meteorological variables M, which can either be pr, ps, SH, RH, rlds, rsds, or T 2m ; R obs is the observation of streamflow (in m 3 /s); R (RF M+1SD) is the random forest simulated streamflow by +1 SD stirring M and other variables; R (RF all var) is the random forest simulated streamflow with original data; stdev (R obs ) represents the SD of R obs .

Hydrological Models
We use four conceptual lumped HMs as candidates in this study, namely the GR4J, HBV, HMETS, and Xinanjiang (XAJ) models. The GR4J model is an upgraded version of the GR3J initially developed by Edijatno et al. (1999) and subsequently enhanced by Perrin et al. (2003). Estimates of precipitation and potential evapotranspiration serve as inputs for the GR4J model, but only streamflow is employed to optimize parameters. Here, we employed a temperature-based approach proposed by Oudin et al. (2005) to estimate potential evapotranspiration in calibrating the model. The HBV model was first applied to hydropower plant flood forecasts in the 1970s by the Hydrologiska Byråns Vattenbalansavdelning (Bergström & Forsman, 1973). The HBV model is a simple lumped HM that stands up in comparisons of models used to simulate snowmelt runoff. HMETS is a basic but efficient ensemble conceptual HM developed by École de technologie supérieure. The HMETS records essential hydrological activities such as flow direction, evapotranspiration, infiltration, snow accumulation, melting and refreezing (Troin et al., 2016). The XAJ model has transitioned from a reservoir flood forecasting model to a multifunctional HM by Zhao (1992). The XAJ model has been used extensively for streamflow simulation in China (Zhang et al., 2012), and most parameters have a clear physical meaning.
The four HMs' parameters are optimized using the Shuffled Complex Evolution (SCE-UA) method developed by Duan et al. (1992), with the objective function of maximizing the Kling-Gupta efficiency (KGE). The SCE-UA algorithm is a commonly used stochastic optimization method for parameter calibration and model optimization. To achieve global optimization of the parameters of the four HMs, multiple simplexes are employed to search the solution space in parallel. Because the daily streamflow series cover different periods across the selected catchments, we choose the 20-year observation period with the best continuity and completeness. We employ a cross-validation approach suggested by Arsenault et al. (2017), in which the calibration (validation) period utilizes data in odd (even), and the first 2 years are used for warm-up. Numerous studies have calibrated the four HMs following the above methods in a large sample of catchments, verifying their promising capacity in streamflow simulation in China (Yin, Guo, Gentine, et al., 2021;.

HM-LSTM Model
In recent years, ML has made remarkable strides in geospatial data modeling (Reichstein et al., 2019). ML excels at fitting high-dimensional data, allowing it to construct trustworthy high precision models, although the "blackbox" feature makes physical processes more challenging to interpret. In order to make the physical processes more interpretable, this study additionally employs the output streamflow of the HM as a constraint of the LSTM to form a HM-LSTM model (i.e., the HM output is one of the input variables of LSTM).
The LSTM is a specific sort of RNN developed to handle the challenges posed by diminishing gradients and carry information from earlier time steps, which is achieved by utilizing backpropagation through time (Hochreiter & Schmidhuber, 1997). The LSTM processes inputs while remembering the states of hidden neurons, passing on information as it propagates forward. One LSTM unit consists of an forget gate, an input gate, an output gate and a memory cell. The forget gate determines what information is discarded and what is maintained in the cell state, whereas the input gate applies a sigmoid function on the prior hidden state and current input. The output gate calculates the next hidden state and outputs the hidden state of the current time step. The following equations describe a typical memory block of an LSTM structure: where x t , i t , f t , and o t are input variables, input gate, forget gate, and output gate at time t, respectively. c t and h t denote the cell state and the hidden unit at the current time t respectively, while c t−1 and h t−1 are from the previous time t − 1. At the initial time, both the cell and hidden states are set to all-zero arrays. b i , b f , b o , and b c are bias items; σ (⋅)and tanh (⋅) denote the sigmoid function and the hyperbolic tangent function, ranging from 0 to 1. Note that we connect a fully-connected layer after the LSTM layer, that is, remap h t to the new output.
To take advantage of the LSTM's ability to handle long time series, the lag times of meteorological variables and simulated streamflow are also considered as inputs of the LSTM. As a result, we obtain four HMs and four HM-LSTM models, which we refer to as eight hybrid terrestrial models (TMs). The outputs of five GCMs under three SSPs are aggregated over the 204 catchments to provide basin mean time series during 1985-2100, which are utilized to drive the calibrated HM to estimate daily streamflow series. We aim to employ the HM-LSTM model to replace traditional HMs that do not perform well in some catchments, thus improving the simulation and explaining the physical mechanisms involved. This multi-model hybrid approach can also increase the number of catchments available and thus enhance the robustness of the study, which is important for projecting future flood risk and the associated socioeconomic exposure.

Calculation of Catchment Lag Time
The implementation of a lag time is prompted by the fact that it takes some time for changes in some meteorological factors to alter the streamflow at the outlet of the catchment; for example, a heavy rainstorm event always requires a certain response time to be translated into the flood peak. This implies that a period of lag time is necessary for training the HM-LSTM model. We use the catchment response time (lag), which measures the amount of time it takes for precipitation to travel from the catchment to the streamflow gauge downstream (Ganguli & Merz, 2019): where A lag (km 2 ) represents the catchment area. For each catchment, we use the meteorological variables during day T and day T-lag to force the TMs. As a summary of Section 3.3, we provide a summarized Table S2 in Supporting Information S1 on the parametrization used to generalize our developed model.

Most Likely Realization and Socioeconomic Exposure Assessment
This study employs flood peak (Q) and duration (D) to explore the joint probability distribution of floods using the copula function, which is a powerful instrument for interdependent variables with known marginal distributions. The flood process is extracted by sampling the annual maximum peak, and the flood duration is identified by using the curve of base flow (Figure 2a). The tenth percentile of daily flow in the study period is determined as the threshold for representing base flow. The marginal distributions of Q and D are fitted with seven candidate distributions (i.e., Gamma, GEV, inverse Gaussian, Log-normal, Normal, Pearson type-III and Weibull). Several types of joint return periods (JRPs, e.g., OR, AND, Kendall dynamic) have been tested within the copula-based framework (Salvadori et al., 2016;Shiau, 2003). Among these, the OR case (T or , i.e., the flood event occurs when either the flood peak or the duration reaches the threshold value) is typically used in flood occurrence evaluation (e.g., Yin et al., 2020): where ; F Q (q) and F D (d) denote the marginal distribution functions of Q and D, respectively.
Under the bivariate copula function, the selection of a suitable T or results in numerous combinations of flood peak and duration (Figure 2b). The points along the T or -level curve usually have different impacts from an engineering perspective (Chebana & Ouarda, 2011), so the probability of each event must be taken into account while determining the appropriate joint quantiles . In this work, the most likely realization method suggested by Salvadori et al. (2011) is employed to select the most probable flood scenario along the T or -level isoline. For a given T or , the following formula can be used to derive the most likely combination: where, f (q,d) represents the joint probability density function of flood peak and duration, f Q (q) and f D (d) are probability density functions of Q and D; μ denotes the mean inter-arrival time between two consecutive floods (in this study μ is equal to 1), and c[F Q (q),F D (d)] = dC(F Q (q),F D (d))/d(F Q (q))d(F D (d)) denotes the density function of the copula, respectively. As a result of the difficulty in generating analytical solutions in Equation 11, a numerical algorithm such as the harmonic mean modified Newton's method can be used to solve this issue .
Previous studies usually define the socioeconomic exposure in the future period as 0 or 100% (e.g., Gu et al., 2020); however, this definition cannot fully represent shifts in future climate hazards. To solve this problem, we define the socioeconomic exposure with consideration the ratio of shifts in JRP: where E POP (E GDP ) characterizes the exposure of population (GDP) affected by increased risk of bivariate flooding; T h and T f are the JRP during historical and future period, respectively; I(·) is the indicator function, when T h − T f > 0 is recorded as 1, and vice versa (i.e., I = 0 when T h − T f ≤ 0); POP and GDP indicate the values of population and GDP in this catchment at the end of the century, respectively.

Uncertainty Contribution Analysis
As China is a flood-prone country, national future socioeconomic exposure has received a lot of attention (H. J. F. Liu et al., 2013). Since the international political conditions and global economic development paths vary greatly under different SSPs, this can directly affect the distribution and size of GDP and population. We partition the overall uncertainty into contributions from the different modules of the integrated model, which consists of three SSPs, five GCMs and eight TMs under ISIMIP3b. In this study, the overall uncertainty is represented by the variance of the updated JRP of the historical  50-year flood during the future (2071-2100) period, and then it is decomposed into the contribution from different sources and their interactions using the multivariate ANOVA method. The change in the climatic indicator is assumed to follow the following model: where μ represents the mean change of the model ensemble of the climatic indicator; R i , G j , and H k represent the effects on the climatic indicator of the ith SSP, jth GCM, and kth TM, respectively; and I i,j,k represents the sum of the effects due to the interaction between different sources.
where VR, VG, and VH present the variance contributed by the effects of SSPs, GCMs, and TMs, respectively; VI RG , VI RH , and VI RGH represent the variance from interaction effects between SSP-GCM, SSP-TM, and SSP-GCM-TM, respectively. By dividing the variance from each source by the total variance, the fractional contributions of individual sources to overall uncertainty are determined.

Machine Learning-Constrained Hydrological Simulations
We first use the random forest model to measure the contribution of each meteorological variable to streamflow, as represented by the sensitivity of streamflow to +1 SD in each predictor variable. Figure 3 shows considerable spatial heterogeneity of the sensitivity results. It is essential to emphasize that the density of stations varies in different regions, and thus the sensitivity gradients vary spatially, with fewer stations in the western half of the country. Precipitation is usually the most dominant factor in streamflow generation in the Yangtze and Pearl River basins, and especially in coastal areas. Most of the variables are positively correlated with streamflow, except for shortwave radiation. The negative relationship between shortwave radiation and streamflow is most pronounced in the upper Yangtze River basin, with a sensitivity of −5% to −15% per SD. The positive feedback of temperature and precipitation on streamflow also implies that streamflow extremes may increase under a warming climate, and intensity of precipitation is a primary determinant of flood risk (Yin et al., 2019).
Among the eight TMs, the XAJ-LSTM model performs the best on average, with the highest KGE over most catchments, while the GR4J and GR4J-LSTM are the best-performing in the majority of catchments models with a total number of 85 selected in the 204 catchments (Figures 4a and 4b). Simultaneously, HBV and HBV-LSTM show best performance only in 8 catchments, probably because they suffered from over-parametrization during the calibration period (Orth et al., 2015). The KGE values of the optimal models are presented in Figures 4c and 4d during calibration and validation periods, respectively. The KGE values during the validation period are larger than 0.7 over 70% of catchments, indicating that the selected model generally demonstrates a satisfactory simulation performance. Even though we observe a few catchments with relatively lower KGE during the validation period, the KGE values of all observed catchments are still larger than 0.55. Therefore, we are confident that the selected the best one in each catchment can model streamflow accurately for all 204 catchments. We also find that the KGE values during the validation period are significantly lower in the northern catchments than in the southern catchments. These findings are consistent with prior research (e.g., , because the precipitation-streamflow relationship is stronger in the humid southern regions and thus facilitates hydrological simulation.

Projected Changes in Univariate Flood Quantiles
After projecting the daily streamflow by driving the best-performing TMs with bias-corrected CMIP6 outputs, we extract the annual maximum peak discharge and corresponding flood duration, and then fit all the distributions mentioned in Section 3.4 with the maximum Bayesian Information Criterion, respectively. Figures 5  and 6 illustrate the multi-model ensemble (five GCMs) mean flood peaks and duration under the 50-year historical return period (RP), with substantial spatial variation throughout China. Due to different catchment sizes, the flood peaks under the 50-year RP range from 2,000 m 3 /s to 20,000 m 3 /s in the historical period. We also depict the relative changes in univariate flood quantiles under four given RPs, revealing that both flood peak and duration are anticipated to increase in the future due to a warmer climate, with substantial regional heterogeneity in magnitude. We detect a strong growing trend in the Yellow River basin, Yangtze River basin, and Songhua River basin for the 50-year RP, which is widely recognized as a vital indication of flood occurrence characteristics in China, suggesting a possible threat to the current infrastructure in those regions under future climate warming. Additionally, a large number of catchments would witness an increase in the 50-year RP by more than 50% in those regions ( Figures 5 and 6), whereas only a few catchments, such as the Jinsha River basin, would experience a decrease in flood risk.
To systematically analyze the shifts in univariate flood quantiles, we provide the relative changes in design flood peaks and durations for the five selected GCMs and RPs under the three SSPs in Figures S5-S34 of Supporting Information S1. Although the GCMs have considerable variation in quantiles, all five GCMs predict an overall upward shift in flood quantiles, implying a strong increase in flood risk. The future climate projections suggest that only a few catchments are likely to experience a minor decrease in flood peak and duration. Overall, the flood conditions in 9 out of 10 studied catchments are projected to worsen under all SSPs, with relative change rates of peak approaching +40%. A few of catchments exhibit opposite tendencies under different GCMs; for instance, catchments near the Liaohe River would experience decreasing flood quantiles under the MPI-ESM1-2-HR GCM (Figures S11-S13 in Supporting Information S1), whereas the remaining four GCMs project a substantial increase. The relative change results for all 204 catchments, for the five GCMs under the 50-year RP are also provided explicitly by the violin plots in Figure 7. The violin depicts the kernel density method's probability distribution of relative change. It is shown that the violin shapes are not sensitive to SSPs, but change in both flood peak and duration across different GCMs, which display significant uncertainty. The relative changes in flood quantiles for various RPs range from −50% to 300%, and only a few catchments exhibit a decreasing tendency under the five GCMs. In addition, the relative change rates are concentrated between 20% and 100% in most catchments, indicating that the flood quantiles in China would substantially intensify. Overall, even with the significant GCM uncertainty and regional heterogeneity, it is anticipated that the flood quantiles assessed for different RPs are likely to increase rapidly as the climate warms.

Bivariate Flood Risk Changes and Socioeconomic Exposure
Meteorological disasters frequently lead to vast socioeconomic losses due to their various characteristic properties (Ayantobo et al., 2017), and it is necessary to synthesize the bivariate flood risk changes. Figure 8 illustrates 9-year. This indicates that it is beneficial to limit carbon emissions to reduce the societal vulnerability to bivariate river floods. Figure 9 further presents the future GDP and population exposure to increased bivariate flood risk under the three SSPs. We find that the GDP exposure is higher in the southeastern coastal region and Central China, while the population exposure is higher in the Pearl River Delta region, the Yangtze River basin, and the Huaihe River basin. It is informative to note that the GDP exposure is lowest under SSP3-70 (24.0 million dollars/km 2 vs. 29.8 million dollars/km 2 under SSP1-26 and 54.7 million dollars/km 2 under SSP5-85), probably due to the stagnation of China's economic growth after 2050 in the SSP3 scenario (T. Jiang et al., 2017). These areas have a large number of industrial and transportation hubs, and frequent flooding could damage their functionality (L. Liu et al., 2020).
In the Yangtze River basin, which has many chemical industries, the large exposure of GDP is also associated with environmental pollution and damage to production (Xia et al., 2021). In the future, this impact may diminish with upgrades in industrial structure and technology, but disregard for environmental issues could have disastrous consequences in the context of climate change (e.g., SSP3). The highest population exposure is also projected under SSP3-70 (607 people/km 2 vs. 361 people/km 2 under SSP1-26 and 420 people/km 2 under SSP5-85), due to the faster growing population, the large number of newborns, the lower life expectancy, and the "pyramidal" shape of the population compared to other paths (T. Jiang et al., 2018). Although different SSPs produce large differences in the magnitude of exposure, the relative spatial distribution of exposure is approximately the same under future warming conditions, which has positive implications for the government to formulate corresponding policies.

Underlying Physical Mechanism and Uncertainty Analysis
By investigating future changes in the meteorological variables, we attempt to comprehensively assess the physical mechanisms behind the increased risk of future flooding in China. Figure 10 ( Figure S35 in Supporting Information S1) illustrates the changes in extreme meteorological variables, defined as changes in the 99th (90th) percentile values from the historical to the future period. As anticipated, near-surface temperatures and precipitation are projected to substantially increase in China under all three scenarios, with the largest increases under SSP 5-85 (Figures 10a-10c and 10g-10i). Given that precipitation is among the most significant factors driving the HMs, increased extreme precipitation would likely explain the increasing floods in most catchments in China. In Section 4.1, we found stronger sensitivity of streamflow to precipitation in Southeast China (Figure 3a), therefore large declines in the JRP (enhancing bivariate flood risk) are projected in these regions. Due to the limits imposed on oceanic water vapor movement, RH decreases in some land regions (Figure 10j-10o), but SH is expected to continue to rise over the majority of the land surface in the future mainly due to the increases in saturation water vapor caused by atmospheric warming , which may intensify future extreme precipitation events. Decreasing streamflow (increased JRP) is projected by the end of this century in the Jinsha River (upstream of the Yangtze River) basin (Figures 8a-8c). This might be explained by the fact that most of the meltwater discharge from glaciers on and around the Tibetan Plateau is expected to peak in the 2030s and 2040s , so the diminishing glaciers might cause declines in snow-generated runoff in downstream areas.
A large collection of hydrological simulations is derived from the three SSPs, five GCMs, and eight TMs, allowing the uncertainty components to be decomposed by ANOVA. As discussed above, the estimates of China's population and GDP also vary considerably under different scenarios, so the SSPs account for a large part of the uncertainty in the exposure of GDP (12.6%) and population (11.3%), as depicted in Figure 11. When evaluating the exposure of GDP to floods, we find that the GCM, GCM-TM, and TM are the top three major uncertainty sources. The population uncertainty results are similar to the GDP results, with the contributions of the GCM, GCM-TM, TM of 23.6%, 22.6%, and 16.4%, respectively. This indicates that the boundary conditions driving the model and the complicated coupling relationship between GCMs and TMs can have a substantial effect on the output. In the uncertainty analysis of both GDP and population exposure, the total influence of the separate factors (i.e., SSP, GCM, TM) accounts for more than 50%. The JRP, on the contrary, is dominated by coupled influences (i.e., SSP-GCM, SSP-TM, GCM-TM, SSP-GCM-TM), which account for 73% of the uncertainty, while the impact of the SSP is minimal (2.8%), especially compared with the significant impact of the SSP-GCM-TM (28.8%). The significant proportion of uncertainty stems from these coupled influences represents that the interaction within the cascade model chain should be paid more attention. For example, combining different models (e.g., GCMs and TMs) may amplify the uncertainty in projecting future hydrological series, which may represent the non-linear uncertainty translation process in a top-down projection framework. Future work should be devoted to quantifying the complicated interaction process from the perspective of physics. Despite the inherent uncertainty, we still project significant increasing tendency of flood risks over China, highlighting that flexible governance measures should be developed to reduce GDP and population exposure.

Discussion
By combining the data-driven and physical HMs in the hybrid HM-LSTM model, we are able to take advantage of the memory structure so that more relevant variables can be taken into consideration. Particularly, the ML can improve the simulation performance in catchments with more intense human activity impacts. Therefore, this hybrid model is able to achieve reliable performance of streamflow simulation and projection, which can also reduce uncertainty by considering the physical factors (e.g., radiation, pressure, and humidity). Moreover, we examine the changes in the JRP of historical 50-year floods in a bivariate framework and find that the JRP decreases in the vast majority of the studied catchments under the three SSPs. The incidence of the historical 50-year flood events is projected to increase 2-10 times on average by 2071-2100 across China under different SSPs. Moreover, we find large geographic differences in the projected GDP and population of flood-prone Chinese regions, providing valuable information for the government to develop appropriate measures. Within the constraints of socioeconomic conditions, future work could employ the WEAP (Water Evaluation And Planning) model to fully investigate the optimal water resources management under climate change. Finally, we analyze the uncertainty of the results. For GDP and POP exposure, the main components of uncertainty are the coupled effects of the SSP, GCM, and TM; for the future JRP, uncertainty is mainly influenced by SSP or GCM or TM change Figure 11. Fractional uncertainty contributions of all sources (shared socioeconomic pathway, global climate model, terrestrial model, and combinations). The bar charts on the right quantify the mean fractional contribution of each source for gross domestic product, POP, and joint return period.
factors alone. Since the uncertainties associated with different GCMs and TMs are substantial, it is reasonable to adopt a multi-model approach. China is actively implementing an ambitious carbon neutral strategy to achieve peak CO 2 emissions by 2030 and carbon neutrality by 2060. Many studies have focused on how to achieve carbon neutrality in the context of climate change, but overlook the impact of climate change on floods and particularly the associated socioeconomic impacts . China's National Disaster Reduction Committee has stated it in its latest Comprehensive National Disaster Prevention and Mitigation Plan (2022) that there are weaknesses in the country's disaster prevention system in the context of global warming. Improving natural disaster preparedness and upgrading infrastructure for agriculture, communication, transport, housing and electricity have become the key to improve environmental and societal resilience in China. In our study, we find that SSP3-70 and SSP5-85 show higher flood risks and larger socioeconomic exposure than SSP1-26, confirming the importance of keeping a low-emission pathway in China. As China is one of the fastest-growing economies, there is capacity to improve investments in disaster mitigation. Policy makers need to balance economic development and climate mitigation, and should pay more attention to areas with high socioeconomic exposure (e.g., Yangtze River basin), improving weather and climate forecasts and infrastructure sustainability under climate change.
Some limitations must be considered. Although we tried to include as many relevant features as possible to train our HM-LSTM model, the uncertainty that arises from human and geophysical factors cannot be neglected. More specifically, the government's climate change and economic policies are a critical driver of future socioeconomic exposure (Ge & Lin., 2021). Currently, it is challenging to find a data set considering future policy factors (e.g., more aggressive fertility policies), as this may require more complex scenarios (J. Cheng et al., 2021) and the existing SSPs are far from adequate. In trying to relate our model to socioeconomics, the assumption of unchanged subsurface conditions is also implicit. In practice, as the population or an economic center move, the environmental conditions (natural or man-made) are also likely to change, and these shifts are not accounted for in our model. Consequently, our variable screening deserves a more comprehensive investigation in future research. Since flood generation mechanisms are also influenced by multiple factors such as changes in underlying surface conditions (Munoz et al., 2018), future research should aim to reveal the coupling between the drivers, and to minimize the noise (e.g., precipitation and RH may contain duplicate information). In addition, the use of a larger number of GCMs may increase the robustness of global and component uncertainty, so mobilizing computational resources to bound the uncertainty of GCMs in future work is also a next priority.
China is facing many contradictions and obstacles, such as faltering economic growth, and an ecological and environmental carrying capacity nearing its upper limit. This work underscores the need to focus on climate change, which poses significant risks to China's population and economic growth. The Sixth Assessment Report from the Intergovernmental Panel on Climate Change (IPCC AR6; IPCC, 2021) shows that each additional increment of global warming will be accompanied by a greater change in extreme events. With atmospheric warming, storms and floods are likely to become more frequent and intense in most regions across the globe. Our study quantifies socioeconomic exposure on this basis and suggests that policy makers should be well prepared to meet the challenges. Future increases in flood risk seems highly likely; nevertheless, adaptable governance solutions can be developed to better prepare for GDP and population exposure. Otherwise, anthropogenic warming-induced flood intensification may place a considerable strain on our infrastructure and ecosystems.

Conclusions
We present the HM-LSTM model and employ a multi-model ensemble under three scenarios (SSP1-26, SSP3-70, and SSP5-85) to project long-term streamflow. Changes in the quantiles of future (2071-2100) flood peaks and durations relative to the historical period  are analyzed and the most likely realization method is employed to obtain a updated JRP for the future period. Future socioeconomic exposures are derived from the changes in bivariate flood risk for uncertainty analysis. The main findings are summarized as follows: 1. Our HM-LSTM model reduces reliance on traditional hydrological observations in contrast to traditional conceptual lumped HMs, and replaces observations with a well-established reanalysis data set. By ensuring that the KGE coefficient of each selected catchment is higher than 0.55 (most catchments have KGE of over 0.7), we reduce the uncertainty in simulated future flood risk projections. 2. We find the magnitudes of flood peaks and durations are projected to increase significantly in most catchments by the end of this century, especially in the Yellow River basin, the Yangtze River basin, and the Songhua River basin. In the bivariate flood framework, the JRP of historical 50-year flood events strongly decreases in most catchments. SSP 5-85 exhibits the most significant changes, with a median JRP of 22.9 years by late of the 21st century. These findings emphasize the urgency of emissions reductions. 3. Population exposure is highly intense under SSP 3-70, and the SSP5-85 projects most severe GDP exposure, suggesting that a significant increase in future socioeconomic effects by flood risk. In view of the severe situation, it is essential to adhere to and implement the Sendai Framework for Disaster Risk Reduction, including strengthening international cooperation and infrastructure development.