Figure 1. Location of CAMELS catchment outlets across the CONUS. Each catchment was classified into one of the six groups based on the annual average fraction of annual precipitation falling as snow (fsnow ; categories were defined at 0.1 intervals, with all catchments having fsnow of at least 0.5 grouped into one category). The 18 major water resources regions over the CONUS are also shown (grey polygons bounded by whites lines) for reference (see also Figure S1 in Supporting Information).
2.2 Identifying streamflow and precipitation extremes
Our approach to quantifying rainfall and streamflow extremes is based on the annual maxima (AMAX) index, one of the most common indices for assessing temporal changes in hydro-climatic extremes [Do et al. , 2017; Kundzewicz et al. , 2004; Ledingham et al. , 2019; Villarini and Smith , 2010; Westra et al. , 2013]. We first processed streamflow data to obtain the magnitude (QMAX index) and the timing (QDOYMAXindex; defined as the ordinal date of annual floods) for each AMAX streamflow event. To reduce the chance of misattributing flood events, we omitted from our data set any years missing more than 15% of daily values (to have a compromise between data coverage and data quality).
We then processed daily precipitation to derive three sets of variables, each representing a different approach to quantifying precipitation extremes (a threshold of 15% missing data was also used to remove years with insufficient data.). The first variable representing precipitation extremes is AMAX precipitation (P), which is defined using the same approach to that of AMAX streamflow. The second precipitation variable allows us to assess whether constraining the timing of precipitation extremes to seasons where the catchments are wet, and when floods are more likely to occur [Ivancic and Shaw , 2015], could improve the relationship between precipitation extremes and floods. This second variable is calculated as the annual maximum precipitation based only on months in which soil moisture was above-average. The third precipitation variable is intended to take into account catchment saturation and snow dynamics, and how they affect the propagation of precipitation intoeffective precipitation [Berghuijs et al. , 2016;Berghuijs et al. , 2019]. We calculated this variable using a coupled soil-snow routine [Berghuijs et al. , 2016; Hock , 2003; Stein et al. , 2020; Woods , 2009] based on daily precipitation, temperature, and evapotranspiration (all readily available in the CAMELS dataset). Details of this routine is provided in the Supporting Information; for further reading, see Stein et al.[2020].
Finally, we calculate the intensity and timing of each of the three precipitation AMAX variables, leading to a total of six new precipitation indices; PMAX and PDOYMAX(for the first precipitation variable), Psm.MAX and Psm.DOYMAX (for the second), and Peff.MAX and Peff.DOYMAX (for the third). Note that evapotranspiration is only available from October 1980 onward, thus Peff.MAX and Peff.DOYMAXare not available for 1980.
2.3 Assessing the correlation between the spatial pattern of changes in precipitation extremes and the spatial pattern of changes in floods
We calculated temporal changes in the magnitude of AMAX streamflow (QMAX) and changes in precipitation extreme intensity (PMAX, Psm.MAX, and Peff.MAX) using normalized Theil–Sen slope [Gudmundsson et al. , 2019; Stahl et al. , 2012] as follows:
\(\tau_{c}=\text{median}\left(\frac{x_{j}-x_{i}}{j-i}\right)\)(1)
\(T_{c}=\frac{\tau_{c}\times 10\text{years}}{{\overset{\overline{}}{x}}_{c}}\times 100\)(2)
where \(\tau_{c}\) is the Theil-Sen slope estimator for catchmentc , which is defined as the median of the average annual difference in AMAX values (x ) between all possible pairs of years. The indices i and j represent year numbers such that i \(\in\) [1, nc -1], j \(\in\)[2, nc ], i < j , andnc is the number of years in the data record (after the screening process described above) for each catchment.\(T_{c}\) is the normalized trend, expressed as a percentage of change per decade relative to the mean of all AMAX values in a catchment (\({\overset{\overline{}}{x}}_{c}\)). This approach leads to fourTc values for each catchment, one for QMAX and three for precipitation intensity (PMAX, Psm.MAX, and Peff.MAX).
To evaluate whether the spatial pattern of changes in floods can be explained by the spatial pattern of changes in precipitation extremes, we calculated the coefficient of determination R2[Rao , 1973], which is the square of the correlation between the Tc values of QMAX and the Tc values of a precipitation extreme metric (e.g., Tc of PMAX). The value of R2 ranges from 0 to 1, and a high R2indicates a strong correlation.
2.4 Assessing the causal relationship between precipitation extremes and floods
To quantify the strength of a potential relationship between precipitation extremes and floods, we calculated the probability that an annual precipitation extremes would be followed closely (in time) by an annual flood extreme in each catchment (referred to as the co-occurrence probability hereafter). We matched the timing of AMAX precipitation events (represented by PDOYMAX, Psm.DOYMAX, and Peff.DOYMAX indices) and the timing of annual floods (represented by QDOYMAXindex) and derive the fraction of annual flood events that can be directly attributed to a preceding precipitation extremes. To account for travel time required for precipitation to reach a catchment outlet, we adopted a previous approach [Ivancic and Shaw , 2015] and allowed a lag of up to 5 days. Specifically, we presume that there is a causal relationship between AMAX precipitation and AMAX streamflow if 0 ≤ QDOYMAX - PDOYMAX ≤ 5. If precipitation extremes and floods were independent, the random chance of a match, on average, would be less than 2%, which is the random chance of QDOYMAX having a value between PDOYMAX and PDOYMAX + 5 (six days) of all possible days in a non-leap year (365 days).
Given the important relationship between the precipitation type (i.e. snow or rain) over much of the CONUS [Berghuijs et al. , 2016], as well as other parts of the world [Blöschl et al. , 2017; Do et al. , 2020a], we also assessed whether relationships between precipitation extremes and floods vary by precipitation type. To achieve this objective, we used the annual average proportion of precipitation that falls as snow, readily available in Addor et al. [2017a], and is denoted as fsnow . Each catchments was classified into one of the six categories; the first five are defined by fsnow values between 0.0 and 0.5 at intervals of 0.1; the sixth category includes all catchments with anfsnow value between 0.5 and 1.0 (see Figure 1). We then assessed whether there are significant differences in the co-occurrence probability of precipitation extremes and floods acrossfsnow classification categories.
2.5 Assessing temporal correlation between the intensity of precipitation extremes and flood magnitude
We identified catchments with similar causal relationship between precipitation extremes and flood by grouping catchments into seven groups according to the co-occurrence probability at 0.1 intervals (note that all catchments with at least 0.6 co-occurrence probability grouped into one category). We then measured the co-variation between the intensity of precipitation extremes (e.g., PMAX index) and the magnitude of floods (QMAX index) at each catchment using the coefficient of determination R2. The R2 values were then analyzed alongside the co-occurrence probability to quantify the extent of which changes in precipitation extremes are useful to infer changes in flood magnitude.
3 Results and Discussions
3.1 A low correlation between the spatial pattern of changes in floods and the spatial pattern of changes in precipitation extremes
Figure 2 shows temporal changes in the magnitude of annual floods and precipitation extremes across the CAMELS catchments. We note that the effective precipitation (Peff) could be equal to zero throughout the year wherever precipitation could not make the catchment fully saturated. Specifically, there are 110 catchments having zero Peff.MAX over more than 20 years, leading to a zero Thiel-Sen slope estimated as shown in Figure 2h (see also Supplementary Figure S2). We also removed 33 catchments that have Peff.MAX equal to zero across all years from our analyses, leading to a sample size of 638 catchments for Peff.MAX assessment.
Over the reference period, more CAMELS catchments (53%) experienced a decrease in QMAX index (Figure 2a), consistent with recent investigations [Do et al. , 2017; Do et al. , 2020b; Gudmundsson et al. , 2019; Hodgkins et al. , 2019;Hodgkins et al. , 2017]. On the contrary, PMAXindex shows an increasing trend (Figure 2b) over the majority of catchments (67%), although some interior water resources regions exhibited a more prominent decreasing trend (e.g. Missouri Region; see Figure S2 in the Supporting Information for trends in annual floods and precipitation extremes over individual regions). There is a low correlation between the spatial pattern of changes in PMAX and the spatial pattern of changes in QMAX (Figure 2f) with an R2 of 0.11, indicating that only 11% of the spatial variation of trends of QMAX can be explained using trends of PMAX.
The spatial pattern of Psm.MAX trends (Figure 2c) is generally consistent with that of PMAX trends, while the spatial pattern of Peff.MAX trends (Figure 2d) shows a substantial difference relative to that of PMAX trends, and appears to be more consistent with the spatial pattern of QMAX trends. The coefficient of determination between precipitation extreme trends and QMAX trends support this finding, with an R2 of 0.06 and 0.17 for Psm.MAX trends and Peff.MAX trends respectively (Figure 2g and Figure 2h). These results are generally expected, as the snow-soil routine underlying Peff.MAXcan be seen as a simple conceptual model that takes into account several catchment processes.