Global Projection of Flood Risk With a Bivariate Framework Under 1.5–3.0°C Warming Levels

Global warming increases the atmospheric water‐holding capacity, consequently altering the frequency, and intensity of extreme hydrological events. River floods characterized by large peak flow or prolonged duration can amplify the risk of social disruption and affect ecosystem stability. However, previous studies have mostly focused on univariate flood magnitude characteristics, such as flood peak or volume, and there is still limited understanding of how these joint flood characteristics (i.e., magnitude and duration) might co‐evolve under different warming levels. Here, we develop a systematical bivariate framework to project future flood risk in 11,528 catchments across the globe. By constructing the joint distribution of flood peak and duration with copulas, we examine global flood risk with a bivariate framework under varying levels of global warming (i.e., within a range of 1.5–3.0°C above pre‐industrial levels). The flood projections are produced by driving five calibrated lumped hydrological models (HMs) using the simulations with bias adjustment of five global climate models (GCMs) under three shared socioeconomic pathways (SSP126, SSP370, and SSP585). On average, global warming from 1.5 to 3.0°C tends to amplify flood peak and lengthen flood duration across almost all continents, but changes are not unidirectional and vary regionally around the globe. The joint return period (JRP) of the historical (1985–2014) 50‐year flood event is projected to decrease to a median with 36 years under a medium emission pathway at the 1.5°C warming level. Finally, we evaluate the drivers of these JRP changes in the future climate and quantify the uncertainty arising from the different GCMs, SSPs, and HMs. Our findings highlight the importance of limiting greenhouse gas emission to slow down global warming and developing climate adaptation strategies to address future flood hazards.


Introduction
The Sixth Assessment Report from the Intergovernmental Panel on Climate Change reported that the global mean surface temperature experienced an increase of 1. 1°C (2011-2020), in comparison to the preindustrial levels   (IPCC, 2023).In July 2023 alone, the global average temperature was 1.54°C above preindustrial period (Tollefson, 2023), reaching the Paris Agreement's proposed threshold to control temperature rise well below 2°C or even within 1.5°C (UNFCCC, 2015).In a warming world, climate and weather-related hazards have become an important global issue, affecting human well-being and the sustainable development of socioecosystems worldwide (Khurana et al., 2022;Klomp, 2016;Toya & Skidmore, 2014).
Flooding is a common destructive hazard and has a wide range of impacts on society and ecosystems, posing significant challenges to sustainable development (M.Jibhakate et al., 2023;Moftakhari et al., 2017;B. T. Pham et al., 2021).The number of individuals residing in flood-prone areas, as detected by satellites, has grown from 58 million to 86 million between 2000 and 2015, and the proportion of people affected by flood disasters worldwide continues to grow every year (Tellman et al., 2021).Given the phenomenon of climate warming, the waterholding capacity of the atmosphere shows an increase of about 7%/K, which leads to more intense precipitation extremes and may result in more frequent flood hazards (Gu et al., 2022).Hence, in a warmer climate, it becomes crucial to investigate the factors driving floods and project their future risks.
Flooding is the result of a combination of natural and social factors.Natural factors, including heavy rainfall (Acosta et al., 2016;Dolojan et al., 2023), high antecedent moisture and base flow (W.R. Berghuijs & Slater, 2023; W. R. Berghuijs et al., 2019), rapid snowmelt (Flanagan et al., 2020) and storm surge (Chao et al., 2021;Genovese & Green, 2015), play a dominant role in generating floods and triggering hazards.For example, an extraordinarily severe storm event occurred in the Henan province in China from 17 to 23 July 2021, resulting a large flooded area and substantial economic losses (Ding et al., 2023).Social factors primarily encompass the impact of human activities on the processes of flood occurrence and control measures (Merz et al., 2021).These include changes in land use (Anderson et al., 2022;Brath et al., 2006), vegetation cover (Sheng et al., 2022), and river engineering (Munoz et al., 2018).In recent years, lots of regions across the globe have experienced an escalation in the frequency and intensity of floods, and a few regional studies have projected that future flood risks would increase in most areas (Bilskie et al., 2022;Ganguli et al., 2020).Under such conditions, a systematic projection of future flood risks under global warming is important for designing climate adaptation strategies.
Numerous studies have utilized bias-corrected global climate models (GCMs) as inputs to force hydrological models (HMs) for projecting future floods changes (e.g., Kang et al., 2023;Tofiq & Güven, 2015;Yin et al., 2020).The spatial scales of studying flood risk under climate change range from the catchment scale (Ma et al., 2021;T. Wang et al., 2022;Wu et al., 2023) to the regional scale (Lin et al., 2018;Roudier et al., 2016;Sairam et al., 2021).Nonetheless, there are limited studies that have explicitly examined global-scale flood assessments under various global warming targets.For instance, Asadieh and Krakauer (2017) showed that the magnitude of high runoff extremes would increase in the northern high latitudes, with a projected doubling under the 4.0°C global warming as compared to the 2.0°C global warming.Alfieri et al. (2017) projected global river flood risk under 1.5, 2.0, and 4.0°C warming, and found that over 70% of the population worldwide would face a five-fold surge in flood impacts at 4.0°C warming level.Liu et al. (2021) employed an observation-constrained multi-model method to explore how global flood magnitudes change by 1.5 and 2.0°C warming.Although these studies have analyzed the global flood evolution pattern under climate warming, global hydrological models (GHMs) employed for this purpose typically have missing or overly simplified descriptions of important hydrological processes, simplifying the flow concentration process or neglecting the interaction between surface water and groundwater, and thus tend to be less-well calibrated for flood simulation (Dottori et al., 2018;Gosling et al., 2017).Hirabayashi et al. (2021) used the annual maximum daily discharge to project flood frequency, and found that flood exposure is projected to increase with climate warming.More importantly, floods are complex events which can be described by multiple characteristics such as peak flow, flood volume, and duration (Yin et al., 2020).Numerous studies have mostly analyzed flood risk by using a univariate framework (Arsenault et al., 2017;Huziy et al., 2013), which fails to capture the multiple aspects of flood events simultaneously and has certain limitations (Requena et al., 2016;Yin et al., 2020).High peak discharge or prolonged duration, which respectively reflect the magnitude and temporal characteristics of floods, pose significant threats to human life and socio-economic activities (Brunner et al., 2019;Woods et al., 2023).Therefore, it is particularly important to consider multiple flood characteristics jointly when assessing the impact of climate warming.Copula functions can be used for constructing joint distributions due to their good ability to link different marginal distributions (Fu & Butler, 2014;Guo et al., 2018;L. Zhang & Singh, 2006).However, most studies of multivariate flood characteristics have focused on regional flood risk (Naseri & Hummel, 2022;R. Yang et al., 2021), leaving a poor understanding in the flood's multiplex nature from a global perspective, especially within the context of global warming.This study primarily emphasizes the utilization of copula within a bivariate framework to assess the joint occurrence of flood peak and duration globally.
The main aim of the study is to examine the projected alterations in bivariate flood characteristics (i.e., peak and duration) over 11,528 catchments across the globe.First, we use temperature data from ERA5 and precipitation data from Multi-Source Weighted-Ensemble Precipitation version 2 (MSWEP-V2) along with daily streamflow to calibrate five HMs.Then, we drive five lumped HMs with the outputs of five bias-corrected GCMs from the latest sixth Coupled Model Inter-comparison Project Phase-6 (CMIP6) under three shared socioeconomic pathways (SSP126, SSP370, and SSP585) to project future daily streamflow for each catchment.After extracting different global warming levels based on multi-model ensemble (MME), the bivariate framework of flood peak and duration is established based on copulas during historical and different climate warming scenarios at the 1.5, 2.0, 2.5, and 3.0°C warming levels.Finally, we use the most likely realization to evaluate the shifts in bivariate flood risk, and systematically unravel the factors behind the changes of joint return period (JRP) under future climate as well as quantify the potential uncertainties.

Observational Hydrological Data
The daily streamflow data is sourced from various regions worldwide (Text S1 in Supporting Information S1) and is available from directly accessible links (Table S1 in Supporting Information S1).We conduct a three-step quality assessment to ensure data integrity and reliability: we initially eliminate the stations that have experienced changes in measurement instruments or station datum in order to maintain data consistency.Subsequently, we retain 12,895 stations having a 20-years valid data minimally, ensuring that each year has less than 10% missing data.We then exclude 1,367 catchments that exhibit inadequate performance in hydrological simulation, as characterized by a Kling-Gupta efficiency (KGE) value of less than 0.5 across all five candidate HMs.Overall, these filtering steps leave 11,528 catchments that span diverse climate regions worldwide.

ERA5-Land Data
The European Center for Medium-Range Weather Forecasts (ECMWF) develops the ERA5-Land data, which is an application of the land portion of the ERA5 climate reanalysis (Malardel et al., 2016), with a consistent global coverage and high resolution (0.1°× 0.1°).To calibrate the HMs, we use the hourly 2 m air temperature (T dry ) to derive daily average, maximum and minimum temperatures.The ERA5-Land data has been extensively employed in risk assessments related to climate change and hydrological extremes, demonstrating a high level of data quality and reliability (Sullivan et al., 2020;Yin, Gentine, et al., 2023;Yin, Guo, Gu, et al., 2021).

MSWEP-V2 Data Set
We employ the precipitation data from MSWEP-V2, which is a high-resolution (0.1°× 0.1°) global precipitation data set covering 1979 to present at a 3-hourly temporal resolution.This meteorological data set is created by ideally integrating an extensive variety of meteorological data observation methods, which has been successfully applied both regionally and globally (Beck et al., 2019).The 3-hourly precipitation is aggregated to the daily scale prior to conducting data analysis.

Coupled Model Inter-Comparison Project Phase-6 (CMIP6) Scenarios
We employ a MME, including five GCMs (i.e., IPSL-CM6A-LR, GFDL-ESM4, MPI-ESM1-2-HR, MRI-ESM2-0, and UKESM1-0-LL) from the latest version of CMIP6, to reduce the uncertainty in analysis.Table S2 in Supporting Information S1 presents the details of these selected GCMs.The simulations of these GCMs under three SSPs (i.e., SSP126, SSP370, SSP585) are used, which represent sustainable pathways, regional rivalry, and fossil-fueled development, respectively (O'Neill et al., 2016;Riahi et al., 2017).The daily climate variables are precipitation, maximum, minimum and average temperature, relative humidity and the 2 m specific humidity, which have been methodically bias-corrected and progressively downscaled to 0.5°× 0.5°with high precision under the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3b, Lange, 2019).Vapor pressure deficit (VPD) can be calculated by using ISIMIP3b temperature and relative humidity.The ISIMIP3b protocol aims to quantify the risk linked to climate change at various socioeconomic and climatic levels.This data set has been used in numerous studies, all of which have demonstrated its reliability for evaluating climate changes (Janssens et al., 2020;Yin, Guo, et al., 2023;Yun et al., 2021).

Methodology
The study's flowchart encompasses five components (Figure 1).First, five HMs are calibrated employing temperature data from ERA5, precipitation data from MSWEP-V2 and observational streamflow.Next, the calibrated HMs are driven by the ISIMIP3b outputs to project historical and future streamflow series, after which the flood peak and duration are extracted.We extract different global warming levels according to the time-sampling Earth's Future 10.1029/2023EF004312 HUANG ET AL. method.Subsequently, the copula function and most likely realization method are integrated to calculate flood risk change under various warming scenarios, and we also analyze the drivers of the future JRP changes.Finally, the analysis of variance (ANOVA) method is applied to evaluate the JRP uncertainty estimates stemming from different procedures.

Definition of a Baseline Period and 1.5°C-3.0°C Global Warming Periods
In the light of earlier research (Döll et al., 2018;Ji et al., 2020;Marx et al., 2018), the time-sampling method is employed to decipher the periods of distinct stages of global warming.Since CMIP6 also produces outputs during the historical period defined as spanning from 1850 to 2014, we select the most recent 30-year interval  as baseline period.The HadCRUT4 data set indicates that the surface average temperature during this historical period is recorded as 0.66°C, relative to preindustrial levels (Morice et al., 2012).Consequently, the warming thresholds of 1.5, 2.0, 2.5, and 3.0°C correspond to temperature increases of 0.84, 1.34, 1.84, and 2.34°C above the reference period, respectively.SSP126 represents a sustainable pathway with a low forcing level of 2.6 W/m 2 by 2100, which envisions the world adopting a green development trajectory, thereby encountering the least significant obstacles in terms of climate change mitigation and adaptation efforts.SSP370 is a central trajectory with a medium-to-high anthropogenic radiative forcing level of 7.0 W/m 2 by 2100.Lastly, SSP585 portrays a fossil-fueled development scenario with the highest radiative forcing level in the future.This indicates that the world would adopt an energy-intensive and fossil-based economic model, thereby facing greater challenges in mitigating climate change.We use three scenarios to represent different development models in the world.For a specific SSP, we define distinct warming periods by utilizing the 30-year moving average of global mean temperature derived from MME. Figure 2 illustrates the periods corresponding to various warming levels for each SSP.Notably, it is important to highlight that under the SSP126 scenario, an increased temperature in the GCM simulations exceed no more than 3°C by the end of the 21st century.The 3°C warming level period is therefore defined as the last 30 years of the century under this scenario, with a view to comprising with other scenarios for convenience.

Hydrological Modeling
Despite recent progress in global flood risk assessment, the accuracy and improvement of GHMs-based flood simulation is insufficient for most applications (T.Yang et al., 2019).Using well calibrated watershed HMs is important for simulating floods or extreme flows reliably at the outlet of the watershed (Veldkamp et al., 2018).Considering that the 11,528 catchment areas span different climate zones and underlying surface conditions worldwide, a single HM may not be able to reliably capture the heterogeneity of hydrological processes.The Xinanjiang (XAJ) model (Zhao, 1992) and GR4J-6 model (Perrin et al., 2003) are two common conceptual HMs, which are suitable for humid and semi-humid areas.Considering that snow covered watersheds are mainly dominated by snowmelt runoff, we employ the Hydrologiska Byråns Vattenbalansavdelning (HBV) model (Bergström & Forsman, 1973) and HMETS model (Martel et al., 2017), which both contain a snow module.The SIMHYD model (Chiew et al., 2002) has advantages in semi-humid and semi-arid areas by considering both saturation excess runoff and infiltration excess runoff mechanisms simultaneously.Therefore, we use these five HMs, with different parameters, modules and mechanisms that can effectively identify different land surface hydrological processes in different climates of the globe.

Model Description
The XAJ model employs mathematical equations based on empirical relationships derived from observed data to estimate various components of the water balance equation, consisting of four primary modules.The GR4J-6 model relies on a simplified representation of natural processes, which incorporates two parts of production storage and routing storage.In the GR4J-6 model, estimates of precipitation and potential evapotranspiration are utilized as input variables.In our study, we apply a temperature-based method in model calibration as described by Oudin et al. (2005).The HBV model was initially developed to generalize the hydrological evolution at a catchment scale, including three main modules, that is, soil moisture routine, degree-day snow routine and runoff routine.Due to its structural simplicity and high performance in capturing key hydrological processes at various temporal scales, the HBV model has been widely applied worldwide for purposes such as flood forecasting, water resource management planning, and climate change impact assessments (Shagega et al., 2020;Shiwakoti, 2017).The HMETS model is designed to simulate the complex interactions between the terrestrial surface and the hydrological cycle.The model utilizes a pair of interconnected reservoirs to represent both the unsaturated and saturated zones within the catchment.It effectively simulates essential hydrological phenomena, including evapotranspiration, infiltration, snow accumulation, as well as melting and refreezing processes (Martel et al., 2017).The SIMHYD model is a seven-parameter conceptual model that simulates infiltration excess runoff, saturation excess runoff and base flow.It contains four stores for groundwater, interception loss, routing process and soil moisture (Y.Zhang et al., 2009).

Model Simulation
The outputs of the lumped model are based on the watershed scale.To ensure simulation accuracy, we calibrate the model using the highest-resolution data set available (i.e., 0.1°× 0.1°ERA5-Land temperature data and 0.1°× 0.1°MSWEP-V2 data set), and the average watershed values are computed using Thiessen polygon.For the future projections, it is challenging to acquire data with a high-resolution of 0.1°× 0.1°.Therefore, we utilize the commonly used high reliable ISIMIP3b data set with 0.5°× 0.5°following numerous previous studies (J.Li et al., 2022;Yin, Gentine, et al., 2023).
The optimization of the parameters in five HMs is achieved by the Shuffled Complex Evolution (SCE-UA) method (Duan et al., 1992).The main objective is to maximize the KGE, which is a widely used index for evaluating and comparing model performance in simulating streamflow.
Assuming Q obs represents the observed streamflow and Q sim represents the simulated streamflow, where γ denotes the Person's linear correlation coefficient between Q obs and Q sim , α is the ratio of standard deviations of Q obs and Q sim , and β is the ratio of mean value of Q obs and Q sim .
The SCE-UA effectively utilizes a blend of deterministic search strategies and random elements, enabling simultaneous benefits from the collected topological response surface information throughout the optimization process (Naeini et al., 2019).This method can quickly search for the global optimal solution in HM parameter optimization.When using the SCE-UA method, it is necessary to set relevant parameters (such as the number of vertices in each composite shape, the number of vertices in each sub composite shape, and the population size) based on the number of parameters to be optimized in a specific HM (Duan et al., 1992).When employing this approach to calibrate parameters in HMs, the objective function is computed using streamflow data from odd years, while streamflow data from even years is utilized for model validation (Arsenault et al., 2017;Yin, Guo, Gentine, et al., 2021).

Joint Return Period and Most Likely Realization Method
Following previous research (Brunner et al., 2019;Kang et al., 2023;Tosunoglu et al., 2020;Vittal et al., 2015), in the extraction of flood characteristics, the annual maximum peak is sampled, while the identification of flood duration is achieved through the utilization of the base flow curve (Figure S1 in Supporting Information S1).
Flood magnitudes are often more strongly related to variations in antecedent base flow than short-term (≤3-day) extreme precipitation (W.R. Berghuijs & Slater, 2023).The threshold of base flow is represented by the 80th Earth's Future 10.1029/2023EF004312 HUANG ET AL.
percentile of daily flow during a given period.The copula function, known for its ability to model interdependencies between variables with known marginal distributions, is employed to analyze the relation between flood peak (Q) and duration (D).When identifying the best-fitting function that adequately captures Q and D within each catchment, the selection is based on the utilization of the smallest Akaike information criterion (Akaike, 1974) and the "OR" criterion (T or ).The T or states that a flood event happens when Q&D reaches the specified threshold value, and can be expressed as the following equation: (2) where and F D (d) represent the marginal distribution functions of Q and D, respectively.
In this study, the optimization process focuses on identifying the most likely realization by maximizing the joint probability density, as demonstrated by Salvadori et al. (2016) and Yin et al. (2022).To determine the most likely combination for a given T or , the following formula can be employed: 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; c[F Q (q), F D (d)] denotes the copula probability density function linked to the marginal distributions, respectively.Due to the analytical challenges for obtaining solutions in Equation 2, it is necessary to employ a numerical algorithm to address this computational challenge (Yin et al., 2018).

Uncertainty Contribution Quantification
The projection in implications of global warming on extreme water events may be significantly impacted by multiple uncertainty sources (H.-M.Wang et al., 2020).We employ a comprehensive impact modeling chain, encompassing 75 distinct scenarios that were generated by combining three SSPs, five GCMs, and five HMs.A multivariate ANOVA method (Addor et al., 2014) is employed to decompose this total uncertainty into individual contributions arising from various aspects and their interplays.We assume that change in the climate uncertainty index Δy l,m,n (i.e., updated JRP of the historical 50-year flood in this study) can be modeled using the following equation: where μ denotes the average variation in the model's total climatic uncertainty index; S l , G m , and H n denote the implication on the uncertainty index of the lth SSP, mth GCM, and nth HM, respectively; and I l,m,n signifies the aggregate of the consequence resulting from the association among various sources.
The total variance (also known as overall uncertainty, V T ) can be divided into discrete contributions coming from various aspects using the multivariate ANOVA approach as shown in the following equation: where V S , V G , V H are the variance produced from the potential impacts of SSPs, GCMs and HMs, respectively; V SG , V SH , V GH and V SGH reflect the variance through coupling impacts involving SSP-GCM, SSP-HM, GCM-HM, and SSP-GCM-HM, respectively.

Projected Changes in Flood Peak and Duration Under 50-Year Return Period
Overall, the HMs exhibit satisfactory performance across most catchments (Figures S2 in Supporting Information S1).During the calibration period, over 85% of the catchments obtain a KGE value higher than 0.6, while about 74% achieve this threshold during the validation period.After projecting 75 (combination of five GCMs, five HMs, and three SSPs) scenarios of daily streamflow by driving five well calibrated HMs with GCM outputs from ISIMIP3b, we employ the most likely realization method mentioned in Section 3.3 to estimate the bivariate flood quantiles at various warming levels, with an extraction of the annual maximum peak flow and relating flood duration.Figures 3 and 4 illustrate the MME median relative change of flood peak and duration at the return period (RP) of 50-year from the baseline period  to different future warming periods, presenting large regional heterogeneity at global scale.In most catchments in China, both flood peak and duration are predicted to surge strikingly, which aligns with previous work (Kang et al., 2023).With increasing temperature levels, the growing trend becomes more significant (from approximately 10% of flood peak relative change at 1.5°C warming period to greater than 30% at 3.0°C warming period, Figure 3), implying a potential threat to the existing infrastructure in China.However, for some regions like central Europe, parts of South America, southern Canada, and inland parts of Asia, there is a slight reduction in flood risk within their respective catchments, with most reductions not exceeding 15%.Recently, in both regional and global analyses, a nonmonotonic relationship has been observed between precipitation and runoff, wherein the intensity of extreme events initially increases with rising temperatures, reaching a maximum or peak point, and subsequently decreases (Gao et al., 2020;Yin, Guo, Gentine, et al., 2021).For example, in the Amazon region, the relative change in flood peak under SSP585 scenario does not exhibit a significant increase compared to SSP126 and SSP370 scenarios, which may be attributed to the presence of a nonmonotonic relationship between temperature and runoff.Our analysis reveals that as an outcome of anticipated climatic warming, both the magnitude of flood peak and the duration of flood events are predicted to increase in most global lands, while a few regions (e.g., central Europe, southern Canada, parts of South America) may experience a decreasing flood risk.

Projected Changes of Bivariate Flood Risk Under Different Warming Levels
Weather and climate hazards are defined as major disturbances to the operations of societies caused by natural physical occurrences, leading to widespread damages that exceed the affected society or community's ability to manage the disaster (Freitas et al., 2022), and integrating the alterations in bivariate flood risk is crucial.There are most catchments across globe experiencing a decrease in the JRP by the 1.5°C warming period, to a median JRP of 44.9-year, 36.3-year, and 33.8-year in SSP126, SSP370, and SSP585 in the 50-year flood event, respectively (Figures 5a-5c).In some high latitudes, the increases of JRP (i.e., decrease in floods) may be attributed to prolonged snow-rain transition during precipitation (W.R. Berghuijs et al., 2014).Taking SSP126 scenario as an example, the overall JRP influenced by flood peak shows a decreasing trend (Figure 5d), whereas the JRP influenced by duration exhibits an increasing trend (Figure 5g); this effect is more pronounced in high latitude under a bivariate framework, leading to an overall JRP increase for many catchments in the future.Alterations in the marginal distributions of flood peak and duration, along with their interaction, primarily contribute to the observed changes in JRPs.If only the flood peak were to undergo a shift, a significant increase in flood risk would occur (Figures 5d-5f).For instance, at a global level, the projected average occurrence of the historical 50-year flood event is anticipated to reduce to around 41-years under SSP126 (approximately 38 years under SSP370 and  5i).Under SSP370 and SSP585, dependence between flood peak and duration results in weak increases in the flood event for most catchments, while almost no effect under SSP126 (Figures 5j-5l).Overall, changes in flood peak are the primary factor altering the flood event.We obtained similar results at warming levels of 2.0, 2.5, and 3.0°C (Figures S3-S5 in Supporting Information S1).
Figure 6 provides a more intuitive visualization of the variations on the flood event in the RP of 50-year under different scenarios.As global warming level increases, the JRP in different catchments around the world tends to decrease on average (Figure S6 in Supporting Information S1).The global average JRP ratio, which can intuitively reflect changes in flood risk to a certain extent, represents the ratio of future JRP to 50 years.The average Figure 6.Joint return period (JRP) ratio for the future 50-year return period relative to the reference period for the 1.5 (a-c), 2.0 (d-f), 2.5 (g-i), and 3.0°C (j-l) warming periods.The JRP ratio represents the ratio of future JRP to 50 years under different shared socioeconomic pathways.
ratio decreases from approximately 0.75 at 1.5°C warming to around 0.5 at 3.0°C warming, implying that the flood risk that has almost doubled.This indicates that limiting carbon emissions and mitigating the rate of global warming are beneficial for reducing the vulnerability of society and ecosystem to bivariate flood events.

Projected Changes in Future Climate Variables Behind the Variations of JRP
By exploring the changes in future climatic variables, we attempt to understand the underlying mechanisms affecting the variations of the JRP.Under the SSP126 scenario, future average precipitation changes exhibit substantial geographical variation (Figures 7a,7d,7g,and 7j).Some parts of southeast Asia and North America are anticipated to witness increased rainfall, which leads to an increase in flood frequency in these regions (Figures 6a, 6d, 6g, and 6j), while other regions such as Australia and Brazil exhibit decreases in precipitation.It is worth noting that due to the complexity of rainfall-runoff relationship and growing human stressors in the Amazon region, precipitation decreases with climate warming but the JRP decreases to 35 years at 3°C warming (Figure 6j).In contrast, there is much greater spatial coherence in SH and VPD.Since there is an exponential increase in saturated water vapor with the rise in surface temperature, SH is projected to increase across the majority of the land surface (Figures 7b,7e,7h,and 7k).The increase in temperature leads to an amplified temperature difference between the atmosphere and the surface.Simultaneously, there is a non-linear change in humidity differences from the land to the air (Barkhordarian et al., 2019).Under this coupling effect, VPD is expected to increase, with atmospheric saturation water vapor rising faster than vapor deficit, thus intensifying precipitation extremes (Figures 7c,7f,7i,and 7l).We obtained similar results under the SSP370 and SSP585 scenario (Figures S7 and S8 in Supporting Information S1).Overall, these climatic variables variations are vital in controlling the hydrothermal conditions and extreme weather events in the future.

Variance Decomposition
Hydrological simulations are generated from a comprehensive ensemble consisting of three SSPs, five GCMs, and five HMs.This extensive ensemble enables the decomposition of uncertainty components through the application of ANOVA.At different warming levels, the overall uncertainty of the future JRP estimates is broken down into uncertainty components from various sources.Figure 8 shows the contribution of seven uncertainty sources and their spatial distribution in the globe for 1.5°C warming period.The JRP is predominantly influenced by coupled factors (i.e., SSP-GCM, SSP-HM, GCM-HM, SSP-GCM-HM), which collectively contribute to 69.80% of the overall uncertainty.GCMs are contributed by various institutions and regarded as one of the major uncertainty sources in coupled factors (Abramowitz et al., 2019).These models undergo meticulous parameter adjustments to refine their performance (Mauritsen et al., 2012), and their outputs are simulated without a comprehensive assessment of internal variability (Maher et al., 2021).
GCMs comprise a compilation of models through diverse organizations that frequently collaborate in information and knowledge sharing (Abramowitz et al., 2019).These models are fine-tuned through intricate parameter adjustments and their simulations are presented lacking an assessment of inner fluctuation (Maher et al., 2021;Mauritsen et al., 2012).Global climate model uncertainty is thus a key factor, especially in central South America, northwestern Europe, and southeastern Asia (Figure 8b).The influence of the SSP alone is found to be minimal (4.23%), particularly when compared to the substantial impact of SSP-GCM-HM (24.63%).The spatial distribution of seven uncertainty contributions for the 2.0, 2.5, 3.0°C warming periods is similar to the 1.5°C warming period (Figures S9-S11 in Supporting Information S1).The interplay among SSPs, GCMs, and HMs is found to be significant in contributing to the overall uncertainty, indicating that more emphasis should be placed on understanding the interaction within the cascade model chain.This highlights the importance of utilizing a comprehensive ensemble of models when projecting future changes in flood risk.
Furthermore, we find that the uncertainty contribution resulting from the interplay between different sources undergoes slight variations with increasing temperatures.Regarding the changes in JRP across different global warming, the whole uncertainty can significantly affected by the selected GCM and the interaction between the GCM and other sources of uncertainty (i.e., GCM, SSP-GCM, GCM-HM, SSP-GCM-HM), accounting for approximately 81.63% (81.03%-82.09% from 1.5 to 3.0°C warming) in the whole uncertainty.It is important to note the interaction contribution between various uncertainty sources is substantial.Collectively, these interactions account for approximately 72.91% (ranging from 69.73% to 74.22%) of the overall uncertainty.This uncertainty contribution from interactions is closely related to the choice of GCMs.Specifically, the contributions from the interactions between GCMs and SSPs, as well as between GCMs and HMs, are approximately 15% and 25% respectively, while the interaction between SSPs and HMs contributes only around 5%.Previous studies have also indicated that, under climate change, the GCM selection is a significant aspect of uncertainty within hydrological predictions (Chen et al., 2023;Shen et al., 2018).In contrast, regardless of global warming level, SSPs only contribute a small portion (approximately 5%) of the uncertainty, as the warming periods for different SSP scenarios, such as 1.5, 2.0, 2.5, and 3.0°C, are defined based on their respective warming trajectories.

Discussion
Certain limitations need to be considered.Despite attempts to gather daily streamflow data from numerous global locations, the predominant concentration of measurements is found in the United States, Europe, Australia, and China.Data from India and Africa are relatively scarce, however, these regions are especially sensitive to the water-related extremes (Alfieri et al., 2017;Mishra et al., 2020).field surveys, the use of remote sensing data, the analysis of historical records, and the establishment of collaborations with local authorities and communities.In bivariate flood risk analysis, the study does not account for the potential influence of changes in human socio-economic factors which may modify flood regimes, such as the construction of new hydraulic engineering and alterations in revegetation (Wen et al., 2021;J. Zhang et al., 2022).Moreover, an incorporation of broader range of GCMs has the potential to enhance the assessment of great uncertainties both globally and regionally.In order to lower the uncertainty around GCMs in next research, it is imperative to prioritize the deployment of computing resources.Future research must investigate the interactions between these drivers while reducing the impact of disturbance, given the sensitivity of flood generation mechanisms to variables such as modifications in conditions of the environment (Munoz et al., 2018).For example, precipitation and relative humidity may contain redundant information, which requires us to carefully screen during the analysis process.
Streamflow simulation and attribution are important scientific endeavors in studying the response of the terrestrial water cycle and water resource changes to global climate change.For streamflow simulation, GHMs might better project the response of flood regimes to climate change for a large range of the globe, but their large spatial scale and strong heterogeneity limit the precise description of hydrological processes (Padiyedath Gopalan et al., 2022;Zaherpour et al., 2018).For streamflow attribution, the main factors influencing streamflow include climatic factors, watershed underlying surface, and human activities.In the future work, global flood projections based on GHMs and identification of dominant factors driving streamflow deserve further investigation.Moreover, in some mid to high-latitude and high-altitude regions, such as northeast China, North Europe, northern America and southwestern South America, the winters are long and particularly cold, with deep snow accumulation.The accumulated snow melts when temperatures rise in the following spring and summer, resulting in snowmelt floods.Under climate warming, the occurrence time, frequency, intensity, and impact of snowmelt flood disasters have changed significantly (Vormoor et al., 2016).Snow-dominated floods merit further analysis and discussion in future work.
Climate change has profound effects at the global scale, which manifest in disrupted weather patterns that pose significant risks to food production and an increasing likelihood of catastrophic flooding (Nicholls et al., 2021;Taherkhani et al., 2020).This study highlights the significance on addressing climate change and examines variations in flood risk for bivariate flood variables across various climate warming scenarios.The anticipated escalation in flood vulnerability across various regions worldwide highlights the imperative to mitigate greenhouse gas emissions, amplify carbon sequestration endeavors, and fortify climate resilience.According to the IPCC AR6 report, the intensification of climate warming is expected to lead to increasingly intricate and difficultto-manage consequences and risks associated with climate change (IPCC, 2023).It is predicted that both the frequency and intensity of storms and floods would rise worldwide.Policymakers are faced with the challenge of finding a careful balance between promoting economic growth and addressing climate change, with particular attention given to regions that are more susceptible to socio-economic vulnerabilities.Efforts should be focused on enhancing weather and climate forecasting capabilities, as well as ensuring the resilience of infrastructure to well respond to the challenges under climate change.Addressing climate change and effectively mitigating future flood risks requires comprehensive cooperation among governments, agencies, and civil society across nations.One example that highlights this matter is the International Joint Commission for the Great Lakes.The IJC's approach involves comprehensive cooperation and coordination among various stakeholders.It facilitates information sharing, scientific research, and policy development to assess the impacts of climate change, improve flood forecasting, and enhance adaptive measures.Failure to jointly address global warming in this way may strain our infrastructure and ecosystems as a consequence of intensified flooding.Furthermore, countries of the Global South (e.g.,China and India), with a large population and underdeveloped economy, are likely to be flood hotspot areas in the future and thus require special attention.Only by achieving the carbon peak and carbon neutrality can we fundamentally alleviate climate warming and the upward pressure of flood disasters.

Conclusions
This study provides a global-scale assessment of the impacts of 1.5, 2.0, 2.5, and 3.0°C global warming above preindustrial levels on river floods.Specifically, we employ five well-calibrated lumped HMs (i.e., HBV, GR4J-6, HMETS, SIMHYD, and XAJ) to simulate the daily streamflow using the outputs from 5 bias-corrected GCMs under three SSPs (i.e., SSP126, SSP370, and SSP585) from CMIP6.We examine the alterations in flood peak and duration in the future under different warming levels in comparison to the historical period .We get the updated JRP using the most likely realization method and evaluate its uncertainty by utilizing ANOVA method.The main conclusions are summarized as follows: Earth's Future 10.1029/2023EF004312 1.After calibrating five HMs with good simulation performance, we extract the peak discharge and duration from the simulated streamflow.It is observed that the majority of catchments experienced an increase in both peak intensity and duration of floods, with a more pronounced trend as temperatures rise; however, there are also regions of the world exhibiting decreases in flood intensity and duration.2. The JRP of 50-year flood events declines in the majority of catchments within the context of bivariate flood framework, with the flood peak acting as the main driver of JRP changes.Notably, SSP585 shows the most largest variations, with the median JRP values of 33.8, 30.0, 26.9, and 23.2 years under 1.5, 2.0, 2.5, and 3.0°C warming scenarios, respectively, but substantial discrepancies in various parts of the world.These findings underscore the urgency of enhancing efforts in climate change adaptation and carbon dioxide reduction.3. The uncertainty analysis reveals that the future JRP estimates are predominantly influenced by coupled factors, especially the association between GCM, SSP and HM, explaining approximately 27.87% (24.63%-29.17%from 1.5 to 3.0°C) of the total climatic uncertainty.GCM uncertainty is also a key factor, especially in central South America, northwestern Europe, and southeastern Asia.In contrast, the contribution of the SSPs and HMs to uncertainty relatively small across different warming scenarios, and their spatial distribution is similar.

Figure 1 .
Figure 1.Flowchart of bivariate flood projection methodology under different warming scenarios.

Figure 2 .
Figure 2. Projected global warming periods for 1.5, 2.0, 2.5, and 3.0°C upon multi-model ensemble mean.Dashed vertical lines indicate the middle year at which a given warming threshold is reached and the black vertical solid lines denote the temperature rise uncertainty of five global climate models.The horizontal dashed gray lines indicate the temperature anomaly with respect to the historical era (1985-2014) for each warming period.The table shows the 30-year periods of different warming levels for each shared socioeconomic pathway.

Figure 3 .
Figure 3. Multi-model ensemble median change in flood peak magnitude for the 1.5, 2.0, 2.5, and 3.0°C warming periods (rows) and three shared socioeconomic pathways (columns) under the return period of 50-year relative to the reference period.Insets show the histogram of the flood peak relative change, with the dashed vertical line representing the median value.

Figure 4 .
Figure 4. Multi-model ensemble median change in flood duration for the 1.5, 2.0, 2.5, and 3.0°C warming periods (rows) and three shared socioeconomic pathways (columns) under the return period of 50-year relative to the reference period.Insets show the histogram of the flood duration relative change, with the dashed vertical line representing the median value.

Figure 5 .
Figure 5. Future joint return period (JRP) of the historical 50-year flood event and attribution of the change in JRP under 1.5°C warming.Updated JRP, (a-c); JRP change induced by flood peak (d-f), duration (g-i), or dependence (j-l).

Figure 7 .
Figure 7. Future changes in different meteorological variables for 1.5 (a-c), 2.0 (d-f), 2.5 (g-i), and 3.0°C (j-l) global warming periods in comparison to the historical period under the SSP126 scenario.The results are derived from the average of five global climate models outputs from the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP3b).

Figure 8 .
Figure 8. Contribution ratios of uncertainty from seven different sources in the 1.5°C warming period.(a-c) Uncertainty in the future joint return period estimates arising from different shared socioeconomic pathways (SSPs) (a), global climate models (GCMs) (b) or hydrological models (HMs) (c).(d-g) Uncertainty from the interaction of SSP and GCM (d), SSP and HM (e), GCM and HM (f) as well as SSP, GCM, and HM (g).(h) boxplot of seven uncertainty sources across the globe.Insets in panels (a-g) show the histogram of the contribution ratio, with the dashed vertical line representing the median value.
Thus conducting additional research in order to analyze floods in regions with limited data in future studies is necessary, which may entail the implementation of HUANG ET AL.