Amplified Extreme Floods and Shifting Flood Mechanisms in the Delaware River Basin in Future Climates

Historical records in the Delaware River Basin reveal complex and spatially diverse flood generating mechanisms influenced by the region's mountains‐to‐plains gradients. This study focuses on predicting future flood hazards and understanding the underlying drivers of changes across the region. Using a process‐based hydrological model, we analyzed the hydrometeorological condition of each historical and future flood event. For each event, at the subbasin scale, we identified the dominant flood generating mechanism, including snowmelt, rain‐on‐snow, short‐duration rain, and long‐duration rain. The rain‐induced floods are further categorized based on the soil's Antecedent Moisture Condition (AMC) before the event, whether dry, normal, or wet. Our historical analysis suggests that rain‐on‐snow is the primary flood mechanism of the Upper Basin. Although most frequent, the magnitude of rain‐on‐snow floods is often less severe than short rain floods. In contrast, historical floods in the Lower Basin are primarily caused by short rain under normal AMC. Given the uncertainties in climate projections, we used an ensemble of future climate scenarios for flood projections. Despite variations in regional climate projections, coherent perspectives emerge: the region will shift toward a warmer, wetter climate, with a projected intensification of extreme floods. The Upper Basin is projected to experience a marked decrease in rain‐on‐snow floods, but a substantial increase in short rain floods with wet AMC. The largest increase in flood magnitude will be driven by short rains with wet AMC in the Upper Basin and by short rains with normal AMC in the Lower Basin.


Introduction
The Delaware River Basin (DRB) is a complex and critical region for investigating flood hazards.With its mountainous areas, coastal plains, and diverse ecosystems, the DRB is particularly susceptible to severe storms that can trigger devastating floods.This susceptibility is compounded by the fact that the DRB is home to over 14 million people across four states, as well as critical coastal infrastructure and sensitive wetland ecosystems.Historical flood events across the DRB exhibit strong spatial variations in flood seasonality (Berghuijs et al., 2016;Smith et al., 2010;Sun et al., 2021;Woltemade et al., 2020), implying that flood mechanisms vary both spatially and temporally.The main flood mechanisms in this region include rain-on-snow (ROS) and snowmelt that mainly occur in the high latitude areas, as well as heavy rainfall from tropical cyclones, extratropical systems, and convective systems.The severity of flooding from these mechanisms can also be notably influenced by the soil's Antecedent Moisture Conditions (AMC).For instance, in January 1996, Pennsylvania experienced major flooding caused by ROS.This was a result of heavy rainfall on saturated ground, coupled with rapid melting of a deep snowpack accumulated over a succession of major snowfalls.Another notable flood event is Hurricane Irene in August 2011 (Sun et al., 2021;Villarini & Smith, 2010;Xiao et al., 2021).Prior rainfall had nearly saturated the soil in the region before Hurricane Irene's landfall, which further exacerbated the flooding (Deb et al., 2023).
The vulnerability of the DRB to severe floods underscores the importance of understanding the impacts of climate change on flood characteristics in the region.Although previous studies have explored this topic (e.g., Sharma et al., 2021;Woltemade et al., 2020), significant knowledge gaps remain regarding the underlying mechanisms that drive climate-induced changes in flood characteristics within the DRB, as well as their implications for the region's vulnerability to flooding under future climate variability.In a warmer climate, it is generally expected that peak snow accumulation will decrease, snowmelt will start earlier, and the transient snow zone will shift toward a rain-dominated regime (Lee et al., 2020;Mote et al., 2018;Musselman et al., 2017;Yan et al., 2019).However, the implication of these changes on future flood potential can differ across regions and vary based on factors such as local topography and basin characteristics.Some studies suggest that flood risk in snow-influenced areas may decline due to a combination of diminished snowpack and earlier snowmelt.This could lead to a slower water release into rivers, as there's less energy available for the melt process (Musselman et al., 2017;Yan et al., 2019).Conversely, other research suggests that winter flood potential may increase as more precipitation falls as rain, making it immediately available for runoff, along with anticipated more intense rainfalls in a warming climate (Cao et al., 2019;Kunkel et al., 2010;Safeeq et al., 2016).The implications of heavy storms on flooding under future climate conditions is also undefined.Recent research has projected an intensifying trend in tropical cyclones in future climates near the U.S. Atlantic coast, particularly for more extreme storms (Balaguru et al., 2022;Walsh et al., 2016;Wehner et al., 2018).More recently, Weaver and Garner (2023) indicated a positive trend in hurricane genesis points moving northward, thereby leading to more landfalls in the mid-Atlantic region.Additionally, uncertainties related to future climate projections introduce further complexities into flood predictions, a challenge that has yet to be addressed.
This study uses a physics-based modeling approach to predict future flood hazards and investigate the flood mechanisms that drive changes in flood characteristics of the DRB under climate change.While previous studies have sought to understand flood mechanisms through observational data sets, they are limited in their ability to diagnose projected changes in floods at the process level.Our approach involved simulating the hydroclimatological conditions for each flood event, allowing for the assessment of the drivers of individual floods under nonstationary climate conditions.Specifically, we configured a process-based watershed-scale hydrological model separately for each 8-digit hydrologic unit code (HUC8) basin within the DRB, which simulates continuous hydrological states at a fine 90 m resolution for both historical and future climate scenarios.By aggregating the hydrometeorological states over each HUC8 basin, we analyzed at the HUC8 level the driving mechanisms for each simulated flood event, defined using annual maximum flows (AMF) at the basin outlet.We chose the HUC8 scale because it is more likely to be dominated by a single flooding mechanism compared to coarser HUC levels, and it is also more relevant for informing management decisions compared to finer scale HUC catchments that focus on tributary systems.By explicitly accounting for climate uncertainties in our flood projections, we assessed how climate uncertainties manifested in projected flood changes attributed to each investigated driving mechanism.Additionally, we explored whether there were common perspectives of future flood changes despite the discrepancies in climate projections to inform future flood management and adaptation strategies.

Study Domain
The Delaware River Basin (DRB) is a significant freshwater source for approximately 5% of the US population.It is located along the east coast of the U.S. and drains an area of about 35,000 km 2 .The basin generally has a humid continental climate with abundant precipitation throughout the year.Hydroclimate conditions vary greatly across the region, ranging from rain-dominated coastal plains to snow-dominated mountainous terrains, with elevations ranging from over 1,200 m to a near-sea level elevation closer to the Delaware Bay.As a result of its location and topography, the basin is vulnerable to devastating floods associated with various flood-generating mechanisms, as noted earlier.The DRB is composed of 12 HUC8 basins, with drainage areas ranging from 1,000 to 5,000 km 2 (Figure 1).Basins at higher elevations are characterized by extensive forest coverage and a large elevational range, while basins at lower elevations have a higher proportion of urban and agricultural areas, reflecting increased human impact (Table 1).Given this study's focus on riverine flooding, we excluded two tide-influenced HUC8 basins from our model domain, including Broadkill-Smyrna and Cohansey-Maurice, both with mean elevations below 15 m.

Modeling of Flood Processes
To simulate flood processes for historical and future climates, we used the Distributed Hydrology Soil Vegetation Model (DHSVM).DHSVM is a spatially distributed, process-based hydrological model for continuous hydrological simulations at the watershed scale.DHSVM solves the full energy and water balance at the grid level to simulate overland and subsurface hydrological processes.Each grid cell of the DHSVM comprises a two-layer canopy model that simulates the interception, transpiration, and evaporation of water from the vegetation canopy as well as snow-canopy interaction, a two-layer mass and energy balance snow model to simulate snow accumulation and melt, and a multi-layer soil model to simulate the movement of water through the soil column.Surface and subsurface runoff are routed downgradient cell-by-cell based on the ground slope until they reach the nearest channel.Once runoff reaches the channel network, DHSVM routes flow through the network using a cascade of linear reservoirs.For more detail, readers are referred to previous literature that provides comprehensive descriptions of the DHSVM model physics (Perkins et al., 2019;Storck et al., 1998;Sun et al., 2015Sun et al., , 2019;;Wigmosta et al., 1994).
For this study, DHSVM was implemented at the 3-hourly time step, matching the resolution of the meteorological input, and a spatial resolution of 90-m.The meteorological data set used for both the historical and future simulations is described in detail in Section 3.2.Other geospatial inputs include vegetation type from USGS National  Land Cover Database for Year 2016 and soil type from U.S. Department of Agriculture (USDA) National Resources Conservation Service State Soil Geographic (STATSGO), with the same 90-m grid resolution.A DHSVM pre-processing tool was used to generate the channel network, which creates river topology and river segment properties based on DEM (Fullerton et al., 2022;Perkins et al., 2019).We initiated DHSVM simulations two years prior to the analysis period  to allow a spin-up period for the model to reach equilibrium.Simulations from this 2-year spin-up were not used for the analysis.
DHSVM was calibrated before producing outputs for analysis.Due to the model's high computational demands for such an extensive region, we adopted a manual parameter tuning approach for calibration.Given the importance of snow as an upstream factor influencing streamflow and floods, calibration of DHSVM started with evaluating its daily snow water equivalent (SWE) simulations.The snow parameters for these simulations were sourced from Sun et al. (2019Sun et al. ( , 2022)), which calibrated these parameters specifically for DHSVM over the mountain ranges of the Conterminous US After evaluating the snow simulations, we used a 5-year span of daily discharge data from USGS gauges (Figure 1) to calibrate DHSVM.In each HUC8 basin with available long-term records from a USGS gauge, DHSVM was individually calibrated.The calibration focused on adjusting the soil lateral hydraulic conductivity and soil porosity parameters, following literature recommendations (Cuo et al., 2011;Du et al., 2013;Sun et al., 2016).Calibration was considered satisfactory when the daily Kling-Gupta Efficiency (KGE; Gupta et al., 2009) met or exceeded our target benchmark of 0.65.For a comprehensive model evaluation and to prevent over-parameterization, we evaluated DHSVM's streamflows against measurements at eight USGS gages (details provided in Section 3.1) over the entire analysis period .We used commonly accepted performance metrics, including bias, Pearson's correlation coefficient, and KGE.KGE is a widely employed composite metric in hydrology, which measures the goodness-of-fit between times series of observations and simulations by analyzing the correlation, variability, and bias, providing a holistic assessment of the model's accuracy.KGE values range from ∞ to 1, with a value closer to one indicating better agreement between simulated and observed data.Earth's Future where r is the linear correlation, α is the variability error, and β is the bias.
The outputs from the calibrated model used for flood analysis include the time series of streamflows at each HUC8 basin outlet and time series of hydrometeorological variables averaged over each HUC8 basin, including rainfall, snowfall, soil moisture, SWE, and snowmelt, as detailed in Section 2.3 for determining the dominant flood mechanism.The flood analysis presented here is based on unregulated flows, despite the presence of several major reservoirs in the DRB to the north, to avoid confounding the influence of reservoirs with hydrometeorological processes in determining the flood generating mechanism.Additionally, previous research has shown that the effect of reservoirs on downstream flood distributions is limited in the DRB (Smith et al., 2010), making the use of unregulated flows reasonable for this analysis.

Flood Mechanism Determination
We defined flood events as the AMF within each water year (WY), spanning from October 1 through September 30.
To analyze floods at the HUC8 level, hydrological simulations were performed independently for each HUC8 using the DHSVM model.The flood series used in the analysis were composed of the AMF series based on DHSVM long-term simulations of daily flow at the outlet of each individual HUC8 basin.The AMF was chosen over another commonly used peaks-over-threshold (POT) for defining floods because it provides consistent flood sample sizes when comparing across basins and between historical and future periods.Additionally, the AMF method ensures that selected flood peaks are independent by avoiding selecting peaks during the recession of preceding flood peaks.
We based our flood mechanism classifications on previous flood studies (Berghuijs et al., 2016;Chegwidden et al., 2020;Stein et al., 2021), and refined the approach to improve our understanding of flood mechanisms at the process level.We used a flow chart (Figure 2) based on key hydrometeorological variables associated with each flood event to attribute the dominant flood mechanism to each event.The identified flood mechanisms include snowmelt, ROS, short-duration rainfall under varying AMC, whether dry (SR ), normal (SR), or wet (SR+), and long-duration rainfall with dry (LR ), normal (LR), or wet (LR+) AMC, or other (i.e., not identified by the flow chart).A key distinction between short and long rain, which is based on Stein et al. (2020Stein et al. ( , 2021)), is whether the maximum daily rainfall exceeds two-thirds of the total rainfall over a 3-day duration.The hydrometeorological variables analyzed for each flood event include soil moisture, rainfall, snowfall, SWE, and snowmelt.Rainfall was taken from meteorological inputs to DHSVM, while the other variables were derived from simulation outputs of DHSVM.These variables were then aggregated from the grid cell levels to the HUC8 basin level by taking daily averages over each basin.While a flood event may be caused by different upstream runoff generation processes, we focused on the dominant flood mechanism for each flood event at the HUC8 basin scale.For rainrelated flood mechanism classifications, the hydrometeorological variables were analyzed over a 3-day time window before each flood event.We used a 3-day time window based on Smith et al. (2010), which indicate a maximum flood wave travel time of 3 days from subbasin to estuary.Additionally, Sun et al. ( 2021) reported a 3day maximum lag from peak rainfall to peak flow based on flow observations from multiple gages across the DRB.For snow-related classifications (ROS and snowmelt), we analyzed the hydrometeorological variables over a 7-day time window prior to each flood event, following the recommendations of Chegwidden et al. (2020).This extended time window is adopted because snowmelt often takes longer to generate floods compared to rainfall.We also tested a 5-day and 10-day time window, which yielded similar classification outcomes as a 7-day window.
We identified the dominant flood mechanism for each flood event based on thresholds set for examined hydrometeorological variables, as outlined in the flow chart (Figure 2).We first evaluated whether snowmelt occurred prior to the flood event, and if so, we classify the event as ROS driven if rain fell on the snowpack and the sum of rain and snowmelt contained at least 20% of snowmelt (Musselman et al., 2018;Sun et al., 2022).If heavy snowmelt occurs without concurrent rain, the event was classified as snowmelt driven.If no snowmelt occurred, we assessed if the flood event resulted from a short-duration, high-intensity rainfall event or a long-duration, lower-intensity rainfall event.We then categorized the AMC as dry, normal, or wet based on the minimum surface-layer soil moisture within 3 days prior to the flood event.If no mechanism could be identified, the "other" class was assigned, and no analysis was performed for that class.

Analyzing Flood Characteristics
Our analysis of flood characteristics focuses on the long-term changes in the frequency (i.e., the number of occurrences), magnitude, and timing of flood events, categorized by their generating mechanisms.We used two metrics to determine the timing of floods: the Seasonality Index (SI) and Mean Date (MD).SI measures the seasonal distribution of flood events, and MD measures the average timing of floods.We computed these metrics using circular statistics (Berghuijs et al., 2016;Sun et al., 2022), as shown in the following equations: (3) is the day of the AMF occurrence relative to October 1st for the ith water year (i.e., D = 1 if the event occurred on October 1st), and n is the total number of water years used in the analysis.SI ranges from 0 to 1, with smaller values indicating less pronounced flood seasonality.MD ranges from 1 to 365 days for a regular year or 366 days for a leap year.An MD value of 1 corresponds to October 1st, which is the start of the water year.For flood events with small SI or greater interannual variability in flood timing, MD is a less accurate representation of the actual mean timing of flood events.

Model Evaluation Data Sets
The evaluation of our model's performance focused on SWE and streamflow.Due to the lack of long-term snow measurements in the region, we used the University of Arizona (UA) gridded (4-km) daily SWE data to benchmark the model's SWE simulations (Zeng et al., 2018).The UA SWE data, available from WY 1982-2016, was developed by assimilating in situ measurements of SWE from hundreds of Snow Telemetry (SNOTEL) stations, snow depth measurements from thousands of National Weather Service Cooperative Observer (COOP) stations (Broxton et al., 2016), and 4-km gridded Parameter Regression on Independent Slopes Model precipitation and temperature data over the Continental U.S (Daly et al., 2000).The robustness of the UA data set in capturing key snow dynamics aspects has been demonstrated through extensive and quantitative evaluations in previous studies (Broxton et al., 2016;Dawson et al., 2018).To evaluate the streamflow simulations, daily streamflow observations were obtained from eight USGS gages located near HUC8 basin outlets (Figure 1).These gages were selected based on the following criteria: (a) they have more than 10 years of continuous daily flow records, and (b) they are located upstream of major reservoirs so that flows are not influenced by reservoir management.

Climate Data Set
We used a gridded land surface meteorological data set at a spatial resolution of 1/16°latitude-longitude to force DHSVM over the historical period from 1980 to 2018.The data set was originally available until 2013 (Livneh et al., 2015a) and later was extended until 2018 (Su et al., 2021a).The data set consists of 3-hourly time series of precipitation, air temperature, wind speed, downward shortwave and longwave radiation, and relative humidity.Wind speed was linearly interpolated using NCEP-NCAR reanalysis data, while the remaining meteorological variables were gridded by interpolating the daily maximum and minimum near-surface temperature and precipitation collected from approximately 20,000 COOP stations.For more details about the data set and gridding methodology, readers are referred to Livneh et al. (2015a).This data set has also been evaluated and applied widely in a large body of published studies, for example, (Bohn et al., 2013;Cao et al., 2016;Sun et al., 2022;Yan et al., 2021).(GCMs) at a 1/16°resolution, under the representative concentration pathways 8.5 emission scenario (Taylor et al., 2012).These GCM outputs were statistically downscaled using the Multivariate Adaptive Constructed Analogs approach and was bias corrected to the historical meteorological data set used in this study (Abatzoglou & Brown, 2012).For each GCM, we calculated the mean annual and seasonal change in precipitation and temperature between the baseline period (WY 1970(WY -1999) ) and the projected future period (WY 2070-2099) for each HUC8 basin.The changes in temperature (ΔT) were calculated as absolute changes (°C) and changes in precipitation (ΔP) as a ratio.All 20 GCMs consistently suggest an increase in annual P up to a ratio of 1.2, equivalent to a 20% increase.They also suggest a warmer temperature ranging between 3.0 and 7.2°C (Figure 3).At the seasonal scale (Figures S1-S4 in Supporting Information S1), temperatures are consistently warmer across all GCMs and seasons, with similar warming levels ranging between 4.7 and 5.6°C.Seasonal precipitation changes vary across GCMs.Although all GCMs project wetter conditions for spring and winter, there is a divergence in projected precipitation changes for summer and fall.Nine out of 20 GCMs indicate a drier summer, and six GCMs foresee a drier fall.For subsequent analyses, we selected 10 GCMs out of the available 20 GCMs.These were chosen based on their ability to capture the range of variations in precipitation (dry-wet) and temperature (warm-cold) changes in the studied region, both annually and seasonally.A summary of the 10 GCMs used in this study can be found in Table S1 in Supporting Information S1.These GCMs were used to drive DHSVM for future hydrological simulations, which informed our future flood analysis.

Model Evaluation
Considering the importance of understanding snow-driven floods, the model's skill was evaluated at the UA data grid resolution (∼4 km) for peak SWE across the snow-covered areas of the DRB from 1982 to 2016, which is the overlapping period of UA SWE data and DHSVM simulations.A grid cell, as defined by the UA data, was considered snow-covered if the mean annual peak SWE exceeded 50 mm.Most of the snow-covered areas were located above an elevation of 360 m above sea level, with peak SWE generally increasing at higher elevations.Figure 4 demonstrates the close agreement between DHSVM-simulated peak SWE and UA peak SWE, with a KGE value of 0.72, a Pearson's correlation coefficient of 0.86, and a percent bias of 17.8%.Averaged over the snow-covered areas, the mean KGE of DHSVM-simulated daily SWE is 0.6.Although the date of peak SWE spans from late January to late March, peak SWE occurs in February for over 83% of the region.DHSVM daily streamflow simulations were evaluated for eight HUC8 basins with long-term flow observations that met the gage selection criteria.This evaluation covered the entire historical simulation period  for all HUC8 basins, except for Crosswicks-Neshaminy, which had overlapping flow records with DHSVM simulations only from 2000 to 2014.The KGE values of daily discharge range from 0.65 to 0.77 (Figure 5).The percentages of bias of AMFs range from 15.0% to 3.4% across basins (Figure S5 in Supporting Information S1).

Historical Flood Analysis
The classification method based on the flow chart successfully attributed over 97% of flood events to a single mechanism in all basins, except for the two low-elevation basins adjacent to the Delaware Bay (Lower Delaware  and Brandywine-Christina), where up to 39% of flood events were classified as "other."This result implies that flood events in lowland areas, characterized by flatter topography and more diverse land cover, are more likely to be influenced by multiple flood mechanisms.As the flow chart is designed to attribute each flood event to a single mechanism, we are unable to identify events driven by multiple mechanisms.Therefore, we excluded the "other" category in subsequent analysis.Overall, our results suggest that historical flood events across all 10 HUC8 basins are driven by various flood mechanisms, with no single mechanism accounting for over 45% of flood occurrences.ROS emerges as the dominant mechanism for HUC8 basins with a mean elevation >350 m, while SR predominates in the remaining HUC8 basins with DRB.Among the eight mechanisms investigated, neither snowmelt nor LR-are identified as key flood mechanisms for any basin.For clarity in our discussions, we've segmented these 10 HUC8 basins into two groups based on their dominant mechanism of historical floods: the Upper Basin and the Lower Basin (Table 1).The Upper Basin comprises the five HUC8 basins with a mean elevation >350 m, including East Branch, Upper Delaware, Lackawaxen, Mongaup-Brodhead, and Lehigh, where ROS is the dominant flood mechanism.The Lower Basin comprises the other five HUC8 basins that are SR dominated, including Schuylkill, Musconetcong, Crosswicks-Neshaminy, Brandywine-Christina, and Lower Delaware.

Flood Mechanisms and Seasonal Patterns
The analysis of flood mechanisms, as illustrated in Figure 6, identifies the top three mechanisms based on their percentage contribution to flood occurrences, calculated as the ratio of events caused by each mechanism to the total number of flood events.A fourth mechanism is included if its contribution to flood events is equal to the third mechanism.
In the ROS-dominated Upper Basin, ROS accounts for 29%-42% of historical flood events across its five HUC8 basins.ROS floods consistently exhibit strong seasonality, occurring from late February through March (based on MD) with SI values ranging from 0.6 to 0.84.Rain-driven floods in the Upper Basin generally occur under normal or dry AMC (i.e., SR, SR , and LR), except in Lehigh where LR+ is a major flood contributor.For Upper Delaware, Lackawaxen, and Mongaup-Brodhead, SR and LR events alternately rank as the second and third largest flood contributors, accounting for 18%-21% and 16%-21% of floods, respectively.In East Branch, SR and LR each contribute 16% of flood events.In Lehigh, SR, SR , and LR+ each account for 16% of flood events.The seasonality of SR floods in the Upper Basin is consistently weak (SI < 0.5), while LR floods show relatively strong seasonality, with MD between mid-June through August and SI ranging from 0.51 to 0.61 across basins.For both East Branch and Lehigh, where SR is a dominant flood contributor, SR floods exhibit strong seasonality (SI > 0.7) with MD in the first week of September.LR+ flood is a dominant flood contributor only for Lehigh and shows strong seasonality with SI of 0.78 and MD in early August.
In the Lower Basin, SR is the dominant flood mechanism when excluding the "other" category for the two lowland estuarine HUC8 basins.Across the HUC8 basins within the Lower Basin, SR accounts for 24%-37% of historical flood events.In Schuylkill and Brandywine-Christina, ROS remains a key flood mechanism, along with SR and SR+.In Musconetcong, floods are entirely driven by rain events with normal to wet AMC; SR, LR, SR+, and LR+ contribute to 24%, 21%, 16%, and 16% of flood events, respectively.In contrast, floods in Crosswicks-Neshaminy and Lower Delaware are driven entirely by rain events with normal or dry AMC.In Crosswicks-Neshaminy, SR, LR, and SR account for 37%, 18%, and 13% of flood events, respectively.In Lower Delaware, SR and SR account for 24% and 13% of flood events, respectively.Floods in the Lower Basin, when driven by mechanisms other than ROS, generally show weak seasonality (SI < 0.4).In contrast, ROS floods maintain strong seasonality (SI > 0.7).Notably, the mean timing of ROS floods (based on MD) in the Lower Basin is in late January, which is about a month earlier than ROS flood timing in the Upper Basin, likely due to warmer winter temperatures.

Flood Magnitude by Mechanism
For each HUC8 basin, we examined the flood magnitudes (represented by peak discharge) attributed to each flood mechanism.Figure 7 shows the magnitudes of historical floods across these basins, categorized by the top three flood mechanisms based on their contribution to historical flood occurrences (see Figure 6 for reference).Short rain events (SR , SR, or SR+) consistently cause the highest average flood magnitude in all HUC8 basins, except Musconetcong.Although LR+ events are among the top three flood mechanisms in only two HUC8 basins and are infrequent in most, they produce the largest flood (i.e., highest peak discharge) in half of the HUC8 basins.6).A fourth mechanism is included if its contribution equals that of the third mechanism.Abbreviations of flood mechanisms and the method for determining these mechanisms can be referenced to Figure 2. Black dot indicates mean flood magnitude averaged over events driven by the same mechanism.Note that the boxplot displays the full range of the data without using whiskers.
Earth's Future   the same basin.In contrast, in other SR dominated HUC8 basins, rain-driven flood events with wet AMC (i.e., LR+ and SR+) result in more severe flooding.For instance, in Schuylkill and Brandywine-Christina, SR+ floods have the highest mean flood magnitudes, while in Musconetcong, LR+ floods have the highest mean flood magnitudes.

Shifts in Flood Mechanisms and Timing
An ensemble of flood projections was generated based on 20 downscaled GCM climate projections (detailed in Section 3.1).For brevity, we will refer to these projections as "GCM scenarios."We analyzed the changes in flood mechanisms by contrasting the contribution percentages of each flood mechanism to the total flood occurrences between the future and baseline periods.As illustrated in Figure 8, positive value suggests an increased contribution of a flood mechanism to future flood events.While the ensemble flood projections present a rather wide range of variations due to inherent uncertainties in climate projections, a consistent trend emerges: the historically ROS-dominated Upper Basin is expected to experience a marked decline in ROS floods and an increase in SR+ events.
Within the Upper Basin, ROS floods are projected to decrease by 12%-23% across its HUC8 basins, according to the average GCM scenarios.The reduction in ROS events is positively correlated with elevation.For instance, East Branch, with the highest mean elevation, is projected to experience the largest decrease in ROS floods, ranging between 7% and 37% across GCM scenarios.Lehigh, at a lower elevation, is also expected to see a ROS flood decrease, although at a gentler rate of 3%-23%.Drawing from MD and SI metrics, most GCM scenarios predict that future ROS flood events will maintain strong seasonality, although occurring about a month earlier, between January and mid-February.A notable increase in SR+ flood events is also projected by most GCM scenarios, with their contribution increasing by 9%-13% across the HUC8 basins based on average predictions across GCM scenarios.Interestingly, these projections also suggest a weaker seasonality of future rain-driven floods, indicating a more dispersed flood timing.
Conversely, in the historically SR dominated Lower Basin, most GCM scenarios indicate less pronounced shifts in flood mechanisms and seasonality.For any mechanism, the average change in their contribution to flood events is predominantly below 5%.Overall, the impact of climate change on flood mechanisms in the Lower Basin is much less significant compared to the Upper Basin.

Increasing Magnitudes of Extreme Floods
We analyze the extreme floods driven by each flood mechanism.The extreme floods are defined based on the 95th percentile of the flood magnitudes related to each flood mechanism.To allow for a standardized comparison across HUC8 basins, the magnitude of extreme floods, denoted as Q95, is normalized by dividing it by the corresponding basin's drainage area.For each HUC8 basin, we compare Q95 driven by each mechanism for the baseline and future periods as predicted by each GCM scenario.The results (Figure 9) suggest that Q95 for most rain-related flood mechanisms are projected to increase under future climate conditions, except for those related to LR-floods, which are projected to decrease in the Lower Basin.
On average, across all GCMs, Q95 of SR floods has the highest magnitude in all basins, with the exception of Mongaup-Broadhead and Lehigh.Moreover, Q95 associated with short rain mechanisms (i.e., SR , SR, SR+) mostly have greater magnitudes than those linked to long rain mechanisms (i.e., LR , LR, LR+), highlighting the more significant impact of short rain events on extreme floods.Interestingly, Q95 for floods with wetter AMC (i.e., SR+, LR+) does not necessarily have higher magnitudes than those with normal AMC (i.e., SR, LR).This implies that changing storm precipitation may play a more important role than AMC in shaping future extreme floods.Yet, Q95 of LR+ floods typically have a higher magnitude across GCMs than that of LR-floods, indicating a strong influence of AMC in this context.
When comparing baseline and future climates, we observe consistent increases in Q95 associated with short rain floods (primarily driven by SR and SR+) across most GCMs.In the Upper Basin where floods are historically dominated by ROS, Q95 of SR+ floods exhibit the highest increase relative to other flood mechanisms.This trend is particularly evident in three basins: the East Branch, Upper Delaware, and Mongaup-Broadhead, where the average increase in Q95 for SR+ floods ranges from 71% to 114% across GCMs.These findings, along with the Earth's Future 10.1029/2023EF003868 SUN ET AL.
projected increase in SR+ flood occurrences, suggest that future short rain events are more likely to occur under wetter AMC, leading to amplified magnitudes of extreme floods in the Upper Basin.On the other hand, while the occurrences of future ROS flood events are predicted to decrease significantly, changes in Q95 for ROS floods are minor, with a range of 3%-9%.In the Lower Basin, Q95 for SR floods are projected to have the highest increase in future climates, ranging from 63% to 113%.In most cases, Q95 for SR+ floods display the second-largest increase in magnitude, following SR floods, with increases ranging from 61% to 86%.Generally, Q95 for floods driven by long rain mechanisms are also expected to increase, but Q95 for LR-driven floods may decrease by up to 32%.

Discussion
Downscaled GCMs employed in our study predict a future for the Delaware River Basin that is consistently warmer and wetter.These climate models inform our hydrological model to generate insights into potential shifts in flood frequency and magnitude driven by various mechanisms in future climates.Yet, inherent variability among GCMs introduces a degree of uncertainty in quantifying these changes.Notably, a minority of GCM scenarios deviate from the consensus, suggesting an increase in ROS flood occurrences in the Upper Delaware River Basin.Such inconsistencies highlight the need for a more nuanced understanding of the climate factors driving these divergent projections.To address this, we examined the sensitivity of extreme flood magnitudes (Q95) to climate change for each flood mechanism at the HUC8 level.We used local polynomial regression to quantify how mean annual changes in temperature (ΔT) and precipitation (ΔP) impact the changes of Q95 between the baseline and future period.Notably, while short-term, localized variations in precipitation often play a more significant role in flooding events, we used annual scale changes in P and T to characterize each GCM, to offer a reference point that aligns more closely with how climate change is typically depicted, not on the event scale, but on a long-term scale.Ultimately, the data set used for this sensitivity analysis included paired ΔP, ΔT, and ΔQ95 for each flood mechanism and each GCM in every HUC8 basin.To ensure a robust sample size for the analysis, we first aggregated this paired data by HUC8 basins, grouping them into the Upper Basin and Lower Basin.Then they were further aggregated by GCMs.Specifically, we used the locfit package (Loader, 2020) in R to fit a smooth surface to the paired data and assessed the goodness of fit between the smoothed surface and the data using the Nash-Sutcliffe efficiency (NSE) coefficient.NSE values range between ∞ and 1.0, where a value closer to 1.0 indicates a better prediction.Additionally, we derived partial derivatives of the smoothed surface with respect to ΔP and ΔT, which is a measure of flood sensitivity to changes in these climatic variables.
Our findings reveal a stronger correlation between flood sensitivity and ΔP compared to ΔT based on NSE.In simpler terms, precipitation changes have a more pronounced effect on future flooding than changes in temperature.This greater precipitation sensitivity is especially pronounced in ROS and SR driven extreme floods (i.e., Q95) (Figure 10).When GCMs predict a ΔP between 0% and 15%, Q95 of ROS floods decreases despite increasing precipitation across a wide range of warming scenarios (up to 8°C).However, when the predicted ΔP exceeds 15%, Q95 for ROS floods increases despite significant warming.A similar trend is observed in Q95 of SR floods, where Q95 decreases for a ΔP less than 10% and increases when ΔP is 10% and higher.Overall, both ROS and SR driven extreme floods show a higher negative sensitivity to precipitation changes under drier GCM scenarios, and a greater positive sensitivity under wetter GCM scenarios.These varying sensitivities of extreme floods to precipitation changes contribute to the observed discrepancies in flood projections among different GCMs.
Lastly, it is important to acknowledge a few limitations in our study.Our analysis relied on a deterministic model, which inherently does not account for the uncertainties arising from model parameterization and structural considerations.Despite the model demonstrating competent performance in capturing flood patterns, it is important to note that these uncertainties remain unaddressed and can be crucial in a holistic understanding of future flood risk.In addition, our study primarily focused on the impact of climate change, overlooking other influential factors such as changes in land use and land cover, as well as water management practices.While this focus allows us to isolate and scrutinize the effects of climate change on floods and their generating mechanisms, our approach could potentially mask the impacts of non-climatic factors that could amplify or mitigate the impact of climate change on floods.

Concluding Remarks
This study analyzed the future flood characteristics and the underlying mechanisms driving future flood changes across the diverse topography of the Delaware River Basin, using a process-based hydrological model and an ensemble of GCM scenarios.The identified flood mechanisms include snowmelt, Rain-on-Snow (ROS), shortduration rainfall under varying AMC, whether dry (SR ), normal (SR), or wet (SR+), and long-duration rainfall with dry (LR ), normal (LR), or wet (LR+) AMC.Our historical flood analysis suggested that ROS was the dominant flood generating mechanism in the Upper Delaware River Basin, characterized by higher elevations.In contrast, SR dominated floods in the Lower Delaware River Basin characterized by lower elevations.In the context of future climate changes, GCMs consistently predict a warmer and wetter future.Most GCM scenarios predict a notable shift in flood mechanisms in the Upper Basin.Specifically, when compared to the baseline period, future floods caused by ROS are projected to decrease by 12%-23%.Meanwhile, flood occurrences driven by SR+ are expected to increase by 9%-13%.Importantly, an increasing trend toward extreme floods is anticipated across the Delaware River Basin, particularly those driven by short rain events (SR, SR , SR+).In the Upper Basin, SR+ floods are projected to have the largest increase in extreme flood magnitude, with an average increase across GCM scenarios ranging from 71% to 114%.In the Lower Basin, SR floods are projected to have the largest increase in magnitude, with an average increase across GCM scenarios ranging from 63% to 113%.Despite these consistent flooding perspectives projected by the majority of GCM scenarios, there are notable discrepancies in projected changes in flood frequency and magnitude among GCM scenarios.This highlights the crucial importance of accounting for uncertainties inherent in climate projections when developing flood management strategies and planning infrastructure to accommodate future climatic conditions.

Figure 1 .
Figure 1.Map of the hydrologic unit code basins in the Delaware River Basin (DRB) grouped into the Upper and Lower Basins alongside the USGS gages used for evaluation of model simulated streamflow and floods.
To account for uncertainties in climate projections and their propagation to flood responses, we analyzed 20 downscaled climate projections from the Coupled Model Intercomparison Project 5 global climate models

Figure 3 .
Figure 3. Mean annual differences in precipitation (P) shown as a ratio, and temperature (T) presented as an absolute difference for the Delaware River Basin (DRB), based on 20 downscaled global climate models (GCMs).GCMs highlighted with circles are selected for future flood analysis to represent the range of variations in P and T changes at annual and seasonal scales.

Figure 4 .
Figure 4. (a) Map showing the mean peak Snow Water Equivalent (SWE) from 1982 to 2016, simulated by DHSVM and smoothed via inverse distance weighted interpolation; (b) Comparison of annual peak SWE between the University of Arizona (UA) gridded SWE data and DHSVM simulations: (c) Comparison of cumulative distribution function (CDF) of daily SWE between the UA data set and DHSVM simulations.

Figure 5 .
Figure 5.Comparison of mean daily flow between DHSVM simulations and USGS flow records for each hydrologic unit code basin with available USGS records.Daily flows were averaged for each calendar day from 1980 to 2018 for all basins, except for Crosswicks-Neshaminy, which was averaged from 2000 to 2014 due to available records.

Figure 6 .
Figure 6.Flood mechanisms in the hydrologic unit code basins of the Delaware River Basin and their percentage contribution to historical flood occurrences.Only the top three mechanisms with the highest contributions are displayed, with a fourth added if its contribution equals that of the third mechanism.In the Upper Basin, rain-on-snow (ROS) is dominant, while in the Lower basin, it is short, intense rain under normal soil antecedent moisture condition (SR).

Figure 7 .
Figure 7. Magnitudes of historical floods across hydrologic unit code basins, categorized by the top three flood mechanisms based on their contribution to flood occurrences (refer to Figure6).A fourth mechanism is included if its contribution equals that of the third mechanism.Abbreviations of flood mechanisms and the method for determining these mechanisms can be referenced to Figure2.Black dot indicates mean flood magnitude averaged over events driven by the same mechanism.Note that the boxplot displays the full range of the data without using whiskers.
dominated HUC8 basins, SR floods generally have the highest mean flood magnitude when compared to the other flood mechanisms within the same basin.While ROS floods occur most frequently, they often have lower magnitudes relative to floods caused by other mechanisms.Specifically, only 9%-19% of ROS floods have a magnitude greater than the 75th percentile (75p) of all flood events, whereas 25%-43% of SR floods have magnitudes greater than the 75p of floods.Transitioning to the SR dominated Lower a variation in the dominance of flood magnitudes is observed.In HUC8 basins such as Crosswicks-Neshaminy and Lower Delaware, SR floods again show the highest mean flood magnitudes when compared to the other mechanisms in

Figure 8 .
Figure 8. Projected changes in flood mechanism within each hydrologic unit code basin, by comparing the percentage contribution of each mechanism to the total number of flood events between the future and baseline periods.Abbreviations of flood mechanisms and the method for determining these mechanisms can be referenced to Figure 2.Each boxplot represents the interquartile range of changes across all global climate model scenarios, with the central white line marking the median, the black dots indicating the mean, and the whiskers extended to the maximum and minimum values.

Figure 9 .
Figure 9.Comparison of extreme flood magnitudes, denoted as Q95 for each flood mechanisms between the baseline and future periods.Here, Q95 represents the normalized 95th percentile of flood magnitudes related to each flood mechanism and is measured in cubic meters per second per square kilometer (cms/km 2 ).Abbreviations of flood mechanisms and the method for determining these mechanisms can be referenced to Figure 2.Each box shows the interquartile range of the Q95 values from global climate model scenarios, with the central line marking the median, black dots marking the mean, and the whiskers extending to the maximum and minimum values.The flood mechanism predicted to experience the largest increase in Q95 is highlighted for each model's projections.

Figure 10 .
Figure 10.Sensitivity of extreme flood magnitude to mean annual precipitation change for (a) ROS floods and (b) SR floods.Here, extreme floods are characterized by Q95, which is calculated as the normalized 95th percentile of flood magnitudes related to each flood mechanism.The values on the contour lines represent the percentage change in Q95 for every 1% change of mean annual precipitation.

Table 1 Characteristics of HUC8 Basins in the Delaware River Basin
Definition and determination of dominant flood mechanisms are described in detail in Section 2.3.In the Upper Basin, flood events are primarily caused by rain-on-snow (ROS), and in the Lower Basin, floods are mainly driven by short-duration, high-intensity rain with normal soil antecedent moisture conditions, terms as SR.Where the USGS gage ID is marked by "x," it indicates that no gage with long-term records is available within the basin.For land use land cover (LULC), the table includes only the three primary types (forest, agriculture, and urban).