Drivers of flash flood frequency and intensity in the United States: a quantitative analysis of hydrometeorological interactions

Inconsistent changes in precipitation and flooding have spurred investigations into the underlying mechanisms, yet the quantitative understanding of interactions between precipitation, temperature, and land cover in streamflow dynamics remains limited. We investigate streamflow changes in 294 small and medium-sized catchments across the contiguous United States (CONUS), using over 30 years of hourly data from the USGS river-watching network. We find that 17.3 % of catchments exhibit significant increases in flash flood frequency (p < 0.05), though not all of these trends remain significant after correction for multiple comparisons; 6.5 % show significant increases in flashiness, while the majority experience no substantial changes. Despite significant upward trends in sub-daily heavy precipitation frequency at 67 % of stations, only 23 % show flood frequency increases, indicating complex catchment-specific hydrometeorological interactions. To quantify the contributions of precipitation, temperature, and land cover changes, we employ a novel time-space varying distributed unit-hydrograph (TS-DUH) model integrated with the DRIVE hydrological model and random forest regression. The results reveal that land cover changes across these predominantly upstream, forested catchments have remained stable over the past four decades, with 90.8 % of catchments showing minimal flow change (within ± 3 %) from 1985 to 2015. Precipitation emerges as the primary driver of streamflow changes, but rising temperature and evapotranspiration mitigate the expected increases in flood frequency and intensity associated with heavier precipitation. Model sensitivity experiments show a 3.6 % reduction in flood frequency and an 8.0 % reduction in flood intensity since the 1980s. Additionally, in rapidly urbanizing catchments, a 10 % increase in impervious surfaces could lead to a 20 % peak flow increase, highlighting the importance of urbanization in flood risk. These findings enhance the understanding of spatial-temporal variation in flash flooding. They also provide crucial insights for better flood hazard mitigation strategies.

1 Introduction

Global warming and human-induced land cover changes are widely acknowledged for their effects on enhancing the hydrological cycles (Winkler et al., 2021; Milly et al., 2008). Observational data consistently show increasing frequency and intensity of heavy precipitation events, thereby amplifying flood risks (Allan and Soden, 2008; Westra et al., 2013; Alexander et al., 2006). Flash floods, which are particularly hazardous due to their rapid onset and intensity, pose significant threats to human safety and socio-economic stability (Ahmadalipour and Moradkhani, 2019; Khajehei et al., 2020). These growing concerns underline the need for a better understanding of flash flood causes and characteristics to develop effective disaster management strategies.

While significant progress has been made in analyzing historical flood changes through long-term discharge records (Merz et al., 2021; Slater and Villarini, 2016; Archfield et al., 2016; Liu et al., 2022), the study of flash flood trends is hindered by the lack of high-quality, long-term sub-daily observational data, particularly for small catchments (Gourley et al., 2013; Gaume and Borga, 2008). This data gap, coupled with short time spans, limited spatial coverage, and infrequent observations, complicates the attribution of changes in flash flood behavior. Flash flood dynamics are influenced by the interaction of multiple drivers within the hydrometeorological system (Hapuarachchi et al., 2011; Borga et al., 2014). The mechanisms of precipitation-runoff processes across various landscapes are complex, and the current understanding remains limited (De Niel and Willems, 2019; Merz et al., 2012). Furthermore, in the context of land surface environmental changes and climate warming driven by human activities, the mechanisms behind flood occurrence have become increasingly complex. Therefore, a comprehensive investigation into the trends and causes of flash flood changes not only enhances our understanding of flash flood behavior but also provides a scientific foundation for future risk management and policy development related especially to flash floods.

In the study of attributing changes in flood characteristics, the roles of climate change and land use/land cover (LULC) changes have been increasingly focused (Wagener et al., 2010). In order to comprehensively evaluate these effects, researchers have widely employed methods such as hydrological models (Frans et al., 2013), statistical analysis (Blum et al., 2020), Budyko theory (Wang and Hejazi, 2011) and machine learning techniques (Jiang et al., 2022; Kemter et al., 2023). Climate-related changes, such as increased precipitation frequency and intensity, are widely recognized as the primary drivers of runoff ratio (Wang et al., 2022) and streamflow change (Slater et al., 2021). It is widely acknowledged that increasing impervious surfaces can intensify the hydrological response to precipitation events (De Niel and Willems, 2019); however, the extent of this impact remains uncertain (Blum et al., 2020; Debbage and Shepherd, 2018; Anderson et al., 2022).

Numerous studies utilizing historical observational data have revealed a lack of systematic increase in flood intensity, despite observed increases in extreme precipitation (Sharma et al., 2018; Archfield et al., 2016; Bloeschl et al., 2017; Do et al., 2017; Liu et al., 2022). The low signal-to-noise ratio in hydrological data, coupled with the inherent complexity of rainfall-runoff processes, hinders the derivation of a unified conclusion to explain this inconsistency (Wasko and Nathan, 2019; Hall et al., 2014). Current research indicates that the inconsistency between changes in precipitation and flooding can be attributed to several factors: rising temperature and decreased soil moisture (Sharma et al., 2018; Wasko and Sharma, 2017; Wasko and Nathan, 2019; Ivancic and Shaw, 2015), although the accelerating global transition toward flash droughts may paradoxically amplify flood response through altered soil hydrological properties when intense precipitation follows prolonged dry periods (Yuan et al., 2023; Barendrecht et al., 2024); the complexity arising from diverse flood generation mechanisms, leading to signal mixing (Zhang et al., 2022); and the varying hydrological responses to precipitation events of differing magnitudes (Brunner et al., 2021). While the impact of rising temperatures on flooding is increasingly recognized, quantitative assessments of flood intensity and frequency remain limited, and the comprehensive effects of LULC require further investigation (Sharma et al., 2018; Ivancic and Shaw, 2015). No consensus has been reached regarding the dominant influence of climate change and land cover change. The majority of research on flood changes utilizes daily data because of the availability of relatively longer time series; however, only a few studies have conducted trend attribution analyses of annual maximum flow using sub-daily data (Kemter et al., 2023; Hodgkins et al., 2019). Understanding abrupt flash floods depends on further study of sub-daily event intensity, frequency, and duration in small and medium watersheds (Carmen Llasat et al., 2016). This necessitates a quantitative analysis of changes in sub-daily flash flood characteristics, and the interactive effects of climate and land cover change on streamflow and flash floods.

In this study, we investigate the trends in flash flood intensity and frequency over the past 30 years (1990–2020) across the contiguous United States (CONUS), utilizing extensive observational instantaneous streamflow data from the United States Geological Survey (USGS) in small and medium-sized catchments (drainage area < 3000 km²). Subsequently, the random forest regression method is used to quantify the impact of climate variables and land cover changes, as well as their interactions, on flash floods. Furthermore, we employ the Dominant River Tracing-Routing Integrated with VIC Environment (DRIVE) model (Wu et al., 2011, 2012a, 2014), alongside a newly developed Time-Space varying Distributed Unit Hydrograph (TS-DUH) model (Hu et al., 2024), driven by DRIVE, to quantitatively assess the impacts of temperature and land cover changes on flash floods. We aim to enhance understanding of flash flood trends and their driving factors, providing valuable insights for improving flood risk management and policy formulation across the CONUS.

This paper is organized as follows: Sect. 2 describes the methodology used for analysis. Section 3 introduces the study catchments and datasets. Section 4 presents the results and analysis, including trends in the magnitude and frequency of heavy precipitation, flash floods, and other variables, as well as their impact on flash floods. Finally, Sect. 5 provides discussion and conclusions based on the study results.

2 Methodology

This section outlines the methods used to detect trends in flash flood characteristics and assess their relationship with hydrometeorological variables. Initially, time series datasets for hydroclimatic variables, including streamflow records, precipitation, air temperature, potential evapotranspiration, and soil moisture, were compiled and subsequently analyzed for trends. To quantify the relative contributions of climatic factors (precipitation, temperature), land cover change, and their interactions to annual maximum streamflow, the Random Forest regression model was employed. Additionally, to explore the link between streamflow and temperature, we use both scaling curves and the DRIVE hydrological model for observed and simulated analyses. The DRIVE model, coupled with the TS-DUH model, was used to quantitatively assess the impact of land cover change on the streamflow response. This multi-method approach aims to provide a mechanistic understanding of the key drivers behind changes in flash flood frequency and intensity.

2.1 Extraction of hydrological and climatic characteristics based on observations

In this study, flash floods are operationally identified as rapid, high-intensity runoff responses of short duration, characterized solely from instantaneous streamflow hydrographs without imposing a rainfall duration criterion (Hapuarachchi et al., 2011). This study focuses on small-to-medium catchments (drainage area < 3000 km²). Flash flood streamflow characteristics were extracted using two approaches: the Annual Maximum Flow (AMF) method and the Peaks Over Threshold (POT) method (Lang et al., 1999). The AMF method refers to the practice of identifying the highest peak flow for each year from a given record of streamflow data. Throughout this study, flood intensity denotes peak flow magnitude (m³ s⁻¹), as quantified by AMF and Peaks-over-threshold magnitude (POTM). For the POT method, the discharge threshold was iteratively determined for each station to achieve an average of two flood events per year, yielding time series data for both the number of flash flood events and peak flows. To ensure inter-event independence and appropriately handle multi-peak hydrographs, events were included in the POT time series only if they satisfied both of the following criteria:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }D>\mathrm{5}+\log \left(\mathrm{0.386}A\right)\text{(2)}& {\displaystyle }& {\displaystyle }{Q}_{min}<{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}min(Q_{\mathrm{1}},Q_{\mathrm{2}})\end{array}$$

Here, Eq. (1) prescribes a minimum inter-event separation (D, in days) that scales with catchment area (A, in km²). Equation (2) requires the minimum inter-peak discharge (Q_min) to fall below 75 % of the lesser adjacent peak (Q₁ and Q₂); when this condition is violated, the peaks are treated as a single event and only the maximum is retained. Unlike the AMF method that ensure at least one event per year, the POT method extracts only flood events that exceed a predetermined threshold. This means in dry years, the maximum flow recorded by AMF may not be identified as a flood event.

Table 1 Summary of hydrological and climatic variables used for flood trends analysis.

Catchment boundaries for each basin were delineated using the hierarchical dominant river tracing (DRT) algorithm (Wu et al., 2011, 2012b), employing 90 m resolution Flow Direction Raster (FDR) data. Basin-averaged temperature (°C), precipitation, and potential evapotranspiration (PET, mm d⁻¹) were subsequently calculated. The hydrometeorological variables used for trend detection are listed in Table 1. A threshold of 95th percentile precipitation was employed to define heavy precipitation events (Mallakpour and Villarini, 2015). To capture differences in sub-daily (1, 6, and 12 h) and daily (24 h) precipitation patterns, we analyzed the annual frequency of heavy precipitation (AFHP) at multiple timescales. These time windows were selected to span the plausible range of catchment concentration times across our study domain. Annual maximum 1 d precipitation (RX1D) was calculated by summing hourly precipitation for each calendar day. Additionally, the annual maximum cumulative N-day precipitation (AMP-N) was analyzed, where n=log (0.386A), A is defined in Eq. (1) (Slater and Villarini, 2016). This formulation adapts the accumulation window to each catchment's response time. The mean duration of flood occurrence (MDF) was calculated as the average time span of each POT event above the threshold. The Peaks-over-threshold flashiness (POTFL) is calculated as the difference between the peak discharge and action stage discharge normalized by the flooding rise time and basin area (Saharia et al., 2017). In addition to streamflow and precipitation metrics, we analyzed annual daily maximum air temperature (AMT, °C), annual total potential evapotranspiration (APET, mm d⁻¹), and annual total soil moisture (ASM, kg m⁻²) to characterize the thermal and moisture conditions influencing flood generation.

2.2 Trend detection

Trend detection in hydroclimatic time series was performed using the non-parametric Mann–Kendall (MK) test (Mann, 1945; Kendall, 1948). The MK test's robustness to non-normality and versatility across diverse datasets make it suitable for analyzing streamflow, precipitation, temperature, and potential evapotranspiration time series at individual stations. Additionally, trends in event occurrence data, including zero-event counts, were further assessed using Poisson regression, a method well-suited to analyzing discrete variables such as flood or heavy rainfall frequency (Mallakpour and Villarini, 2015; Slater and Villarini, 2016). Statistical significance was evaluated at the 5 % level.

2.3 The DRIVE and TS-DUH models

The DRIVE model is a distributed hydrological model developed by coupling the Variable Infiltration Capacity (VIC) land surface model (Liang et al., 1996) and the Dominant River Tracing based Routing (DRTR) model (Wu et al., 2012a, b, 2014). The Global Flood Monitoring System (GFMS) has consistently monitored and forecasted flood occurrences and developments for over a decade, utilizing the DRIVE model informed by satellite precipitation data (Hou et al., 2014; Kirschbaum et al., 2017). Within the DRIVE model, the runoff generation process is derived based on a variable infiltration curve that considers the subgrid heterogeneity of infiltration capacity and multiple soil infiltration layers (Cherkauer et al., 2003; Liang et al., 1996). River discharge, surface water storage, water level, and other streamflow variables are derived for each time step by solving the kinematic wave equations for each river path within the hierarchically organized river network defined by DRT algorithm (Wu et al., 2014). The DRIVE model has undergone extensive justification and validation worldwide and has demonstrated reliable performance (Wu et al., 2019, 2014, 2017; Nanding et al., 2021; Huang et al., 2021, 2023; Yan et al., 2020). In this study, the DRIVE model is run at a 0.125° spatial resolution with a 3-hourly time step. The hydrography dataset (e.g., FDR, Strahler order, flow accumulation) is derived using the DRT algorithm. To incorporate land cover information into DRIVE, the high resolution land cover data (details in Sect. 3) were reclassified to match the land cover classification scheme used in DRIVE by aggregating the original fine categories into broader land cover types via a cross-walking table, and subsequently processed to determine the fractional area of each land cover type within each 0.125° model grid cell. These fractional areas were then used to derive grid-averaged vegetation parameters within DRIVE, such as leaf area index, albedo, and roughness length, from the VIC parameter library, reflecting the composite land cover characteristics of each grid cell.

Our newly developed TS-DUH model, integrated with and driven by runoff outputs from the DRIVE model, was employed to refine flood peak estimation under different LUCC scenarios and to contribute to the separation of precipitation and LUCC impacts on streamflow. The TS-DUH model is designed for efficient flash flood analysis by estimating flow velocity and travel time based on topographic and hydroclimatic characteristics, emphasizing the dynamic contributions of runoff linked to variations in rainfall and soil moisture (Hu et al., 2024). The TS-DUH model operates at a significantly finer spatial resolution of 90 m. It leverages high-resolution DEM data for detailed flow path delineation and topographic analysis. In TS-DUH, Manning's roughness coefficients (n) for each grid cell are derived as the area-weighted average of n values assigned to the underlying 30 m land cover pixels. These n values determine spatially distributed flow velocities, thereby capturing the influence of land surface characteristics on flood wave propagation at a finer scale.

2.4 Analysis of potential driving factors

A random forest regression model (Breiman, 2001) was developed using AP, AFHP6h, AFHP12h, RX1D, RX6H, AMP-N, APET, AMT, ASM and TS-DUH simulated peak flows under varying land cover conditions (with values linearly interpolated for intervening years to represent land cover change) for each catchment to identify the relative importance of different drivers of AMF changes. Random Forest is an ensemble method that improves prediction accuracy and robustness by aggregating multiple decision trees, effectively modeling complex non-linear relationships. For each catchment, hyperparameters were independently optimized using the Tree-structured Parzen Estimator (TPE) algorithm (Bergstra et al., 2011), with the data split into 80 % training and 20 % validation sets and the R² on the validation set maximized over 100 evaluations. The optimized hyperparameters include n_estimators (50–500), max_depth (None or 10–40 in steps of 10), min_samples_split (2–20), and min_samples_leaf (1–10). Finally, the mean decrease in impurity (MDI) scores were then used to rank variable contributions. For each catchment, the inflection point of the cumulative curve of ranked MDI values (the rank of maximum absolute second-order difference) was used to guide the choice of a single fixed MDI threshold, applied uniformly across all catchments to separate primary from marginal drivers. These primary drivers were further categorized into three types: precipitation (including AP, AFHP6h, AFHP12h, AMP-N, RX1D and RX6H), land cover (TS-DUH peak flows), and temperature (AMT, ASM and APET).

In order to investigate the drivers of observed trends, this study focuses on precipitation, temperature, and land cover change. Firstly, streamflow (Q) and precipitation (P) scaling curves with temperature $\left(T\right)-Q\sim T$ and P∼T-were developed to characterize the relationships between these variables (Wasko and Sharma, 2017; Zhang et al., 2022; Yin et al., 2018). For each catchment, daily maximum air temperatures were averaged, and standardized peak flows (at catchment outlets) and maximum 3-hourly accumulated precipitation (for each POT event) were grouped into equally spaced temperature bins. The 99th percentile of streamflow and precipitation within each bin was calculated, and these conditional quantiles were connected to form the scaling curves. The scaling curves are generated using 2 °C temperature bins with a 1 °C overlap to reduce the sensitivity to bin boundary selection (Lenderink and Van Meijgaard, 2010). Furthermore, events with temperatures below 5 °C were excluded to minimize snowmelt effects (Wasko and Sharma, 2017). Secondly, we use the DRIVE hydrological model to simulate the impact of temperature changes on streamflow (Wu et al., 2011, 2012a, 2014). The DRIVE model was run separately with dynamic temperature data (1981–2020) (referred to as DRIVE-DT) and static temperature data (repeating only the 1981 3-hourly air temperature for each simulation year, while all other meteorological forcing fields remain dynamic) (referred to as DRIVE-ST) at 3-hourly temporal resolution and a 0.125° spatial resolution. Both simulations shared initial conditions; differences in multi-year average streamflow and POTF/POTM were then compared.

Land cover data from 1985 and 2015, reclassified using the International Geosphere-Biosphere Programme (IGBP) classification (Schaperow et al., 2021), were used to drive the DRIVE model, assessing the impact of land use/land cover change (LUCC) on runoff and multi-year average streamflow. The sensitivity of the DRIVE model parameters to land cover has been validated (Jiang et al., 2020). Subsequently, the TS-DUH was applied to further analyze the influence of LUCC on flood peaks during the routing process. Furthermore, the TS-DUH was used to simulate POT events to separate the impacts of precipitation and LUCC on streamflow. The TS-DUH captures the hydrological response of the catchment by accounting for subsurface characteristics, as well as the spatial and temporal dynamics of precipitation and soil moisture distribution (Hu et al., 2024).

Figure 1 Trends in annual frequency of heavy precipitation (AFHP) for AFHP1h (a), AFHP6h (b), AFHP12h (c), and AFHP24h (d), representing 1, 6, 12, and 24 h accumulation periods, respectively. Red represents upward trends, while blue indicates downward trends. Areas outlined in black represent trends that are significant at the 5 % level (both increases and decreases), with slopes estimated using the Theil–Sen estimator.

3 Study area and data

All instantaneous streamflow data used in this study were resampled to a uniform hourly interval to provide a consistent basis for sub-daily analysis. The selection of gauging stations for this study followed a multi-stage filtering protocol to construct a robust dataset suitable for sub-daily flash flood analysis. An initial screening identified 2492 stations matching the DRT river network (Wu et al., 2014) from the total of 21 344 USGS gauges with an area of less than 3000 km² (Hu et al., 2024). This pool was then subjected to a stringent data completeness assessment. Recognizing the challenges in obtaining complete year-round hourly instantaneous streamflow records, this assessment focused on the May to September period, which typically exhibits high flash flood activity in many parts of CONUS. A seasonal year was deemed complete only if no gap exceeding six consecutive hours of missing streamflow occurred during the entire May-September window. Stations with fewer than 20 complete seasonal years were subsequently excluded. This criterion balances station density with record quality. Once a station met this selection criterion, all available data across its entire period of record were utilized for the trend analysis to maximize the length of the time series. Furthermore, stations identified as experiencing long-term drought conditions or those with an entire year of unobserved data within the potential analysis period were excluded. This rigorous screening process finally resulted in 294 stations, covering both eastern and western regions of CONUS. No interpolation or gap-filling procedures were applied to the retained records. Notably, no pre-screening was conducted to exclude regulated catchments. The focus on small to medium-sized headwater catchments inherently minimizes the influence of major reservoir operations, though minor regulation effects at some stations cannot be entirely ruled out. The drainage areas of these 294 catchments range from 78.8 to 2985.8 km², with a median of 885.6 km² and a mean of 1024.2 km². Most of the resulting data spanned 1990–2020, with the analysis periods aligned with the available streamflow data.

This study utilized the phase 2 of the North American Land Data Assimilation System (NLDAS-2) precipitation product (Xia et al., 2012) for precipitation trends analysis and to force the DRIVE model, at a 0.125° resolution. Soil moisture data in 0–100 cm were also obtained from NLDAS-2. Atmospheric forcing data (air temperature and wind speed) for the DRIVE model and air temperature trend analysis came from NASA's Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) (Gelaro et al., 2017), at a 0.5° × 0.625° resolution. The Global Land-Cover product with Fine Classification system at 30 m (GLC_FC30) (Zhang et al., 2021) provided 30 m resolution land cover data for 1985 and 2015. Potential evapotranspiration data were obtained from the Global Land Evaporation Amsterdam Model (GLEAM, version4.1a) at a 0.25° resolution. Finally, catchment characteristics, including the Digital Elevation Model (DEM) and FDR data consistent with Hu et al. (2024), were used to extract the TS-DUH.

4 Results

4.1 Trends and spatial patterns of heavy precipitation characteristics

Heavy precipitation frequency trends vary across different time scales (Fig. 1). Significant upward trends in 1 h (AFHP1h) and 6 h (AFHP6h) heavy precipitation frequency are observed at approximately 67 % of the stations. However, this proportion decreases to 30 % for 24 h accumulations (AFHP24h). This discrepancy may be attributed to the inherently localized and intense nature of short-duration heavy precipitation events, which can mask consistent trends at longer timescales. Regional analysis reveals increasing trends in heavy precipitation frequency across the eastern CONUS, suggesting a heightened risk of flash flooding, while western regions show decreasing trends.

Figure 2 Trends in annual total precipitation (AP, a), annual maximum cumulative N-day precipitation (AMP-N, b), annual maximum 1 d precipitation (RX-1D, c), and annual maximum 6-hourly precipitation (RX6H, d). The color bar descriptions are consistent with those in Fig. 1.

Regarding intensity, significant increasing trends are observed in RX1D at 12.6 % of stations and in RX6H at 33.7 % (Fig. 2). Spatially, the distribution of these intensity trends is largely consistent with the observed changes in heavy precipitation frequency. Specifically, the increase in RX1D and RX6H near the Gulf of Mexico is considerably greater than that of AMP-N, likely due to increased hurricane influences (Emanuel, 2005; Balaguru et al., 2023).

Figure 3 Trends in annual total soil moisture (ASM, a), and annual total potential evapotranspiration (APET, b). The color bar descriptions are consistent with those in Fig. 1.

4.2 Trends and spatial patterns of soil moisture and PET characteristics

Soil moisture (ASM) trends exhibit varied patterns across the studied catchments (Fig. 3a). Approximately 10.9 % of catchments show a significant downward trend in ASM, while 3.7 % exhibit a significant upward trend. Overall, a decreasing trend in ASM is observed in approximately 65.3 % of catchments, with the remaining 34.7 % showing an increasing trend, irrespective of statistical significance. This widespread drying trend in antecedent soil moisture may significantly impact runoff generation. Consequently, a decreasing trend in ASM may act as a dampening effect on flood magnitudes and frequencies. Concurrently, a significant increase in APET was observed in over half of the catchments (58.2 %), which is typically associated with rising temperatures and a heightened demand for atmospheric water vapor (Dai et al., 2004). While increased precipitation may have an indirect effect on APET, analysis reveals a mean correlation of −0.35 between AP and APET across all catchments, and −0.39 in the eastern region (east of 100° W), suggesting a weak inverse relationship. Increased precipitation generally leads to higher humidity and lower vapor pressure deficit (VPD), reducing atmospheric demand for water. Nevertheless, the overall effect also depends on local conditions, including air temperature, wind speed, and radiation (Li et al., 2022). The combined effect of increasing APET and overall decreasing soil moisture highlights a shifting hydrological regime that could influence flood processes.

Figure 4 Trends in Peaks-over-threshold magnitude (POTM, a), Peaks-over-threshold frequency (POTF, b), mean duration of flood occurrence (MDF, c), annual maximum flow (AMF, d) and Peaks-over-threshold flashiness (POTFL, e). The color bar descriptions are consistent with those in Fig. 1.

Figure 5 Spatial distribution of the drivers influencing annual maximum flow (AMF) trends. Colors indicate the dominant drivers of AMF trends: precipitation (P, as a green circle), land cover (as an orange square), temperature (T, as a pink triangle), and their combinations. Gray points indicate areas where no significant driver has been identified (N/A).

4.3 Trends and spatial patterns of flash flood characteristics

Analysis of flash flood frequency and intensity trends reveals considerable spatial variability among catchments (Fig. 4). 17.3 % of catchments exhibit a significant upward trend in POTF, increasing by an average of 0.5 events per decade; conversely, 5.1 % show a significant downward trend, with an average decline of 0.4 events per decade. AMF displays limited statistically significant changes in most catchments, only 7.1 % exhibit a significant increase (average 2.5 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{decade}}^{-\mathrm{1}}$), and 1.4 % show a significant decrease (average 3.3 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{decade}}^{-\mathrm{1}}$). Additionally, 3.7 % of catchments exhibit a significant increase in POTM, 3.1 % show a significant decrease. Approximately 6.5 % of catchments show a significant upward trend in POTFL, with 3.4 % exhibiting a significant downward trend. Moreover, MDF shows a significant increase in 4.1 % of catchments (average 1.5 h decade⁻¹), while showing a significant decrease in only 0.3 % (average 0.12 h decade⁻¹). The observed MDF is generally on the order of tens of hours (median 21 h). While Sect. 4.1 revealed that approximately 67 % of stations experienced significant upward trends in AFHP6h, only 23 % of catchments exhibited a corresponding significant increase in POTF, suggesting complex catchment-scale interactions among precipitation, temperature, and land cover. Furthermore, statistical trends for catchments with drainage areas less than 1000 km² (169 stations) showed only minor numerical differences compared to the overall trends.

Regionally, the spatial distribution of POTF trends mirrors the distribution of heavy precipitation frequency, consistent with our previous global-scale study (Yan et al., 2020) demonstrating a stronger correlation between flood frequency and daily precipitation frequency, especially at mid- to high-range daily frequencies. While heavy precipitation frequency significantly increased across much of the eastern CONUS (Fig. 1), significant POTF increases were concentrated near the Great Lakes and in the Appalachian Mountains (Fig. 4). Conversely, a limited number of catchments exhibit a significant decrease in POTF in the western CONUS, consistent with the reduced heavy precipitation trends observed in Sect. 4.1. The changes in flash flood intensity (AMF and POTM) exhibit a more mixed pattern, with POTM showing a higher number of declining catchments compared to AMF. The upward trend in MDF is most pronounced in the Great Lakes and northwestern CONUS. This is likely linked to warmer temperatures accelerating snowmelt, leading to more frequent and prolonged spring floods. The combination of increased precipitation, winter snowpack, and spring snowmelt in the Great Lakes region extends flood duration. Similarly, in the mountainous Northwest, combined snow and rain prolong spring flood events.

Figure 6 Scaling curves of standardized peak flow (Q) and precipitation (P) with temperature for Peaks-over-threshold (POT) events below (a1, b1, c1) and above (a2, b2, c2) the 90th percentile across CONUS, the non-arid regions of CONUS, and the arid regions of CONUS. Circular markers highlight the 99th percentile, and solid lines depict the moving average.

4.4 Identifying major drivers of AMF trends using random forest regression

Across the 294 catchments, the single-variable MDI importance at the inflection point has a median of 0.171 (IQR 0.117–0.232). A fixed threshold of 0.15, which falls within this range, was therefore adopted and applied uniformly across all catchments. Random Forest regression analysis indicates that precipitation is the dominant driver of AMF trends across the studied catchments. Approximately 63.9 % of the catchments show that AMF trends are primarily driven by precipitation. However, land cover change and potential evapotranspiration/temperature can also play a significant role in certain catchments, with about 26.2 % of catchments influenced by both precipitation and temperature (P&T). A complete breakdown of all driver categories is provided in Table S1 in the Supplement. Spatially, P&T-dominated catchments tend to be more prevalent in the western CONUS, whereas P-dominated catchments are proportionally more concentrated in the humid eastern CONUS. To further distinguish these two groups, we compare their baseline climate and land cover characteristics. Compared to the 188 P-dominated catchments, the 77 P&T-dominated catchments receive approximately 17 % less mean annual precipitation (median of 982 mm versus 1184 mm) and are cooler on average (14.5 °C versus 15.9 °C). Despite comparable warming rates (0.269 and 0.291 °C decade⁻¹), similar rates of soil moisture decline, and little difference in land cover composition between the two groups, the lower precipitation in P&T-dominated catchments means that temperature-driven increases in evaporative demand consume a proportionally larger share of the water budget, amplifying the sensitivity of runoff generation to warming. Rising air temperature elevates the atmospheric vapor pressure deficit, enhancing evapotranspiration and depleting soil moisture more rapidly during inter-storm periods. The drier antecedent conditions, consistent with the widespread decline in ASM documented in Sect. 4.3, increase available soil storage at storm onset so that a larger fraction of rainfall infiltrates rather than generating direct runoff, thereby attenuating flood peaks.

Figure 7 (a) Spatial distribution of the percentage difference in multi-year average streamflow (1985–2015) between dynamic (DRIVE-DT; using observed temperature from 1981–2020) and static temperature (DRIVE-ST; recycling the complete 1981 3-hourly air temperature series for each simulation year) scenarios (calculated as ((DRIVE-DT − DRIVE-ST) $/$ DRIVE-ST) × 100 %); (b) histogram showing the frequency distribution of the percentage differences in multi-year average flow.

Figure 8 (a) Spatial distribution of percentage differences in multi-year average flow between 1985 and 2015 land cover scenarios, modelled by the DRIVE model; (b) Frequency distribution of percentage differences in multi-year average flow; (c) Spatial distribution of TS-DUH peak change percentage between 1985 and 2015; (d) Frequency distribution of TS-DUH peak change percentage.

4.5 The impact of air temperature changes on streamflow

Figure 6 illustrates the relationship between streamflow (Q), precipitation (P), and temperature (T). The scaling curves for moderate- and high-intensity flash floods, defined by the 90th percentile threshold, are presented for three regions: the entire CONUS, the non-arid regions of CONUS, and the arid regions of CONUS (Beck et al., 2018). Limitations in sample size preclude a detailed analysis across hydroclimatic zones, and events with temperature exceeding 30 °C are excluded. P generally increases between 5 and 30 °C during both moderate- and high-intensity flood events, except in arid regions. For moderate events (below the 90th percentile), while P and Q increase with T, Q declines between approximately 25 and 30 °C across the CONUS, despite minimal precipitation change. This decline accounts for approximately 11.6 % of the total flood events (12852) across the CONUS. This reduction might stem from enhanced evapotranspiration and increased soil water retention. Conversely, for high-intensity events (above the 90th percentile), Q and P generally increase with T. Between 25 and 30 °C, high peak flows persist despite stable or slightly reduced precipitation, likely due to concentrated extreme rainfall and efficient runoff.

In order to further quantify the impact of temperature changes, we compare multi-year average streamflows (1985–2015) from dynamic (DRIVE-DT) and static (DRIVE-ST) temperature scenarios (Fig. 7). Under this experimental design where only temperature varies while precipitation remains identical, this comparison reveals that 63.6 % of catchments experience a flow change between −5 % and 0 %, while 16.7 % experience a decrease exceeding −8 % (Fig. 7b). DRIVE model simulations indicate that from 1981 to 2015–2020, rising temperatures lead to a 3.6 % reduction in catchments with significantly increased POTF and an 8.0 % reduction in those with significantly increased POTM in annual averages.

4.6 Impact of land cover changes on streamflow between 1985 and 2015

Analysis of land cover change between 1985 and 2015 reveals that urban areas comprised an average of 3.6 % across all study stations by 2015. Forest cover is significantly greater than urban areas, averaging 64.4 %. Most study areas are mountainous upstream catchments with minimal urban land and predominantly forest cover. Overall land cover change was minimal, with the total area undergoing any class transition averaging only 3.3 %. The primary net shifts were urban expansion (from 2.8 % to 3.6 %) and a slight decline in forest cover (from 65.2 % to 64.4 %), while all other classes changed by less than 0.3 %. DRIVE model simulations of land cover change between 1985 and 2015 (Fig. 8a and b) show that the multi-year average flow changes at 90.8 % of observation points fall within ± 3 %, indicating minimal impact compared to temperature variations (Sect. 4.5). Similarly, 82.7 % of gauges show TS-DUH peak changes between ± 5 % (Fig. 8c and d), suggesting relatively stable variation. However, some gauges that show increases exceeding 9 % are influenced by various factors such as snowmelt or urban expansion. Overall, the model suggests that temperature changes have a greater impact than land cover changes.

5 Discussion and conclusions

This study analyzed long-term (30 years) USGS hourly streamflow data from small to medium-sized catchments (≤ 3000 km²) to assess historical trends in flash flood intensity and frequency across the CONUS. It investigated the influences of precipitation, temperature, and land cover change on streamflow by employing multiple methodologies, including the DRIVE, TS-DUH, and random forest regression models, to quantitatively attribute their main drivers and impacts on flash floods. Specifically, we quantitatively assessed the modulatory effects of climate warming on flash flood occurrences and the possible impacts of land cover changes, providing a mechanistic understanding of the primary drivers behind changes in flash flood frequency and intensity. The main conclusions are as follows:

• 1.

  Observations reveal an increase in heavy precipitation events, but this change does not uniformly translate to increase in flood frequency and intensity; 17.3 % of catchments show significant increase in flash flood frequency (p < 0.05), and 6.5 % (3.7 %) show significant increase in POTFL (POTM). Specifically, while approximately 67 % of catchments exhibit a significant upward trend in AFHP6h, only about 23 % of these same catchments show a corresponding significant increase in POTF.

• 2.

  In our model experiments isolating temperature effects, rising temperature tends to decrease both the magnitude and frequency for general floods. DRIVE model simulations demonstrate that since the 1980s, there is a 3.6 % reduction in catchments with significantly increased POTF and an 8.0 % reduction in those with significantly increased POTM.

• 3.

  Land cover changes across the CONUS (1985–2015) are modest overall. DRIVE model simulations indicate that the multi-year average flow changes are within ± 3 % at 90.8 % of the catchments.

The limited translation of widespread precipitation intensification into proportional flood increases likely reflects the combined buffering effects of declining antecedent soil moisture, warming-enhanced evapotranspiration, and catchment storage capacity. In order to examine whether larger catchments show stronger associations with longer-duration precipitation, we computed Spearman rank correlations between each precipitation index and flood metrics, stratified into four drainage-area classes (< 500, 500–1000, 1000–2000, and > 2000 km²). Median correlations remain stable across size classes for all three flood metrics, with ranges across groups typically below 0.05 and no systematic strengthening of longer duration indices in larger catchments (Fig. S1 in the Supplement). Sensitivity analysis across alternative thresholds (p < 0.01, 0.05, 0.10) further confirms that the asymmetry between increasing and decreasing trends is consistent regardless of the threshold adopted (Table S2 in the Supplement).

Several caveats related to statistical robustness merit discussion. Applying trend tests across 294 catchments and multiple variables introduces the potential for false positives; approximately 15 spurious rejections per variable would be expected at p = 0.05 under the null hypothesis. The Benjamini-Hochberg FDR correction was therefore applied (Table S3 in the Supplement). The heavy precipitation frequency trends prove highly robust, with 190 of 195 significant AFHP6h upward trends surviving correction at q = 0.05, and 22 of 51 POTF upward trends remaining significant. Regarding process-level limitations, a long-term decline in ASM indicates a shift toward drier antecedent conditions prior to flood-producing events. However, ASM cannot fully capture the role of antecedent soil moisture in runoff generation. Event-scale attribution linking antecedent moisture windows to individual flood responses therefore warrants further investigation. Furthermore, the study did not consider the evolution of river channels.

The influence of land cover change on flood intensity varies markedly with the degree of urbanization. Since most catchments are located in the less urbanized upstream mountainous areas, the impact of land cover change on flood intensity tends to be low. However, in the rapidly urbanizing catchments, the impact can be substantial. For example, in the Peachtree Creek catchment in Atlanta (Debbage and Shepherd, 2018), the impervious surface fraction increased from 45.2 % in 1985 to 54.8 % in 2015, leading to a 20 % rise (from 57.8 to 69.2 m³ s⁻¹) in TS-DUH peak flow. This highlights that land cover changes can significantly influence flood intensity in certain areas.

Finally, several data and modeling uncertainties should be acknowledged. Data quality, record length, and completeness influence the reliability of trend analyses. Although the May-to-September window was used solely to screen data completeness at the station level, this criterion may favour catchments whose records are most complete during the warm season. Consequently, it may underrepresent basins where rain-on-snow or autumn extratropical storms contribute appreciably to flash flood activity. In the RF attribution, inter-correlated predictors share MDI scores, which dilutes individual variable importance. Aggregating importances into three broad driver categories before catchment classification renders within-category redistribution inconsequential. Cross-category redistribution cannot be excluded but is expected to be minor given the distinct physical nature of the three driver groups. Furthermore, the 0.125° resolution of NLDAS-2 provides limited grid coverage for the smallest catchments, potentially smoothing localized convective peaks and introducing uncertainty in extreme precipitation frequency estimates. This constraint is less consequential for the controlled DRIVE experiments, where identical precipitation forcing is applied across scenarios and only the relative response to temperature or land cover perturbations is compared. Finally, we acknowledge that the inherent limitations of the hydrological models and the uncertainties arising from different data sources warrant further parameter calibration.