We used globally available climate datasets to increase the generalizability of the research. Catchment-scale annual and long-term average Snow Persistence (SP), Aridity Index (AI), and Seasonality Index (SI) were calculated for each catchment. These indices’ spatial maps and summary statistics are provided in Figure 1 and Supplementary Table S1.

To calculate SP at each catchment during the study period, the 500 m resolution—MODIS/Terra 8-Day L3 Version 6— snow cover product (MOD10A2) was used (Hall and Riggs, 2016; Hammond, 2020b; a). This data source features improved snow cover detection accuracy and represents an iteration of NASA MODIS data that attempts to correct previously discovered issues (Riggs et al., 2016). At each 500 m resolution grid, the annual SP value was determined as the fraction of days with snow on the ground from January 1-July 3, generally covering peak snow accumulation and complete snow ablation for many parts of North America (Moore et al., 2015). Gridded annual SP data were then clipped using a catchment boundary polygon (Section 2.2) to calculate catchment-scale annual SP. The long-term average catchment-scale SP value for each catchment was then calculated by temporally averaging over the study period.

The aridity index (AI) is a measure of the relative “dryness” of a given region. It is defined by the ratio between potential evapotranspiration and total precipitation (Budyko, 1974). It can be interpreted by focusing on the critical value of 1. As defined in Budyko (1974), when AI is above 1, a region is considered “arid”. The aridity index has been shown in previous continental-scale work to be a key factor in classifying hydrologic behavior (Kuentz et al., 2017). The seasonality index (SI) is a dimensionless measure of the synchronicity between daily precipitation and temperature in a single variable, mathematically defined in Eq. 14 of Woods (2009). At the extremes, −1 (out-of-phase) signifies winter-dominant precipitation and +1 (in-phase) signifies summer-dominant precipitation (Berghuijs et al., 2014). Catchment-scale annual (and long-term average) aridity and seasonality indices were calculated using catchment boundary polygons (Section 2.2) with monthly and daily temperature and total precipitation data from the ERA5-Land database (Muñoz-Sabater et al., 2021) and monthly potential evapotranspiration data released as part of the TerraClimate database (Abatzoglou et al., 2018).

A well-distributed set of stream gauges across North America were selected for this study. We started with daily streamflow observations at 671 and 1,602 stream gauges across the United States and Canada, respectively. These daily streamflow observations were analyzed to generate annual and long-term values of streamflow signatures (see Section 2.4). The United States set of stream gauges and their corresponding catchment boundary polygons were obtained from the CAMELS dataset (Newman et al., 2015; Addor et al., 2017). Addor et al. (2017) developed and validated catchment boundary polygons for each gauge; with each catchment having two sets of area estimates based on two different area estimation methodologies. We queried the catchments listed in the CAMELS dataset to get updated data from the U.S. Geological Survey (2020). The Canadian gauges were derived from the HYDAT dataset released by Environment and Climate Change Canada (2021). We derived flow direction and flow accumulation and ultimately catchments’ area and boundary polygons for the Canadian gauges using the Watershed Tool in ArcGIS. The calculated areas of Canadian catchments were then compared against area information from topographic maps by the Water Survey of Canada (Environment and Climate Change Canada, 2021). Thus, like the United States catchments, there were two sets of area estimates available for the Canadian catchments. Note that differences in area estimation between two alternative approaches stem from different elevation data sources, as well as different algorithms to calculate the boundary polygons. Catchments used from both countries generally feature minimal human disruption (e.g., stream diversions and flow regulations).

From the original set of catchments, only catchments with the following criteria were selected and used in the statistical analyses.

a) Catchments with 15 or more years of acceptable streamflow records during the study period. “Acceptable” years have fewer than 10% of days with missing streamflow data.

Following this screening, a total of 1,187 out of 2,273 original catchments were selected, featuring a wide range of catchment sizes (from under 10 km² to over 273,000 km²), elevations (from about 14 m to over 3,500 m), and precipitations (average annual precipitations of just below 0.3 m/year to about 3.17 m/year) from across North America (Janssen and Ameli, 2021). Most importantly, these catchments represent a variety of North American climatic regions across wide gradients of snow persistence, aridity, and seasonality (Figures 1A–C).

Shape-based streamflow signatures represent the shape and flashiness of streamflow hydrographs, providing insights into streamflow duration and variability at different streamflow levels (McMillan, 2020; Janssen and Ameli, 2021; McMillan, 2021). In this study, seven shape-based streamflow signatures were used to assess the ability of SP to inform streamflow shape and flashiness at low-flow and high-flow conditions. Baseflow index (BFI), normalized Q₅, the slope of the flow duration curve at low-flow (Low-FDC), and low-flow event duration are used to represent low-flow conditions, while normalized Q₉₅, the slope of the flow duration curve at high-flow (High-FDC), and high-flow event duration are used to represent high-flow conditions. Each signature reflects one aspect of low-flow or high-flow behavior, therefore, together these seven signatures reflect the form and functionality of hydrographs at high and low flows.

The baseflow index (BFI) describes the proportion of total flow that come from baseflow, where baseflow is defined as flow that is maintained between precipitation events (Addor et al., 2017). Note that an α of 0.925 was used to calculate the baseflow index using the Lyne and Hollick digital filter algorithm developed by Ladson et al. (2013), with the relevant methodological details contained in their paper. Having meaningful baseflow, especially during dry periods and meteorological droughts, supports ecosystems (Poff et al., 1997; Boulton, 2003) and human activities, especially in arid regions (Saedi et al., 2022).

Q₅ represents the 5% quantile of flow (i.e., extreme low-flow) and Q₉₅ represents the 95% quantile of flow (i.e., extreme high-flow). Here, these values were divided by Q_mean to produce normalized Q₅ and normalized Q₉₅, respectively. Both signatures serve as dimensionless shape-based metrics that describe the relative values of extreme low-flow and extreme high-flow while being agnostic to the average flow magnitude of the catchment (Janssen and Ameli, 2021). Conceptually, they reflect the contributions from the deepest (or shallowest) flow paths to streamflow relative to the contributions from the remaining flow paths of a catchment. Thus, they measure the extreme behavior of catchments, directly translating to events such as droughts or floods in a given stream.

The flow duration curve (FDC) plots recorded streamflow values against their exceedance probability (Searcy, 1959). The slope of the FDC can then be taken within a given range of flow percentiles to describe the variability of streamflow within the specified range (Sawicz et al., 2011). We used two slopes to represent the flow variability during low-flow (Low-FDC) and high-flow (High-FDC). The former represents flow variability between the 5th and 30th flow percentiles and the latter represents flow variability between the 70th and 95th flow percentiles (see Supplementary Table S1 for further details). These flow percentiles reflect the behavior of low-flow and high-flow without capturing overlapping information from the middle portion of the flow duration curve (Sawicz et al., 2011). Large FDC slopes indicate more variable, unstable, and flashier streamflow within the given range of flow, whereas smaller slopes indicate damped and stable streamflow variability. Larger Low-FDC indicates a flashy low-flow hydrograph, implying that streamflow quickly becomes extremely low during dry periods. At the same time, a larger value of High-FDC reflects a flashy high-flow hydrograph, meaning that streamflow rises rapidly after water input events, leading to high-flow instability during wet periods.

High-flow and low-flow event durations express the average temporal length of high and low-flow events, respectively (Addor et al., 2017). We quantified the number of consecutive days with flow larger (smaller) than a flow threshold to calculate high-flow (low-flow) event duration in each year and along the entire time-series. Here, we used a threshold of 0.4x Q_mean to detect low-flow events and a threshold of 4x Q_median to detect high-flow events, where Q_mean refers to average streamflow and Q_median refers to median streamflow. Note that our high-flow detection algorithm is not limited to the detection of only extreme high-flow events. Similarly, our low-flow detection algorithm is not limited to the detection of only extreme low-flow events. Indeed, events with flow magnitudes much smaller than Q₉₅ and much larger than Q₅ can still be detected as high-flow events and low-flow events, respectively.

We calculated annual and long-term values of each signature for each catchment during the study period (2001–2019) using the code and methodology developed by Addor et al. (2017). Details of the mathematical/methodological basis of the code from Addor et al. (2017) are summarized in Table 3 of their study. The spatial maps of the long-term values of streamflow signatures used in our study are shown in Figure 1. Supplementary Table S1 reports some statistics on these signatures.

We conducted three major sets of statistical analyses in this paper. First, through an inter-catchment (spatial) analysis across the 1,187 study catchments, long-term average values of snow persistence (SP), aridity index (AI), seasonality index (SI) and long-term values of streamflow signatures were used to explore how the spatial linkages between SP and shape-based signatures at stream low and high-flow conditions differ among distinct climatic regions (Section 2.5.1). Secondly, through an inter-annual correlational analysis within each snow-rich study catchment and along the study period (2001–2019), annual values of SP and streamflow signatures were used to quantify the catchment-level temporal association between inter-annual variation of SP and inter-annual variation of each streamflow signature (Section 2.5.2). Thirdly, through a standardized inter-catchment and inter-annual correlational analysis across all snow-rich catchments and along the study period (2001–2019), we explored how in a “typical” snow-rich catchment, inter-annual variation of SP is associated with inter-annual variation of streamflow signatures, and how such associations differ between relatively arid/wet years or years with relatively in-phase/out of phase seasonalities (Section 2.5.3). As part of the Supplementary Material, we ranked the associations between SP and streamflow signatures compared to the associations between AI (and SI) and streamflow signatures. In doing so, correlational analyses with annual-level streamflow signatures and climatic indices (SP, AI, and SI) across all study catchments (∼20,000 datapoints for each streamflow signature/climatic index pair) were conducted (see Supplementary Figure S1). We believe that this combination of large sample sizes, climatic diversity, and a suite of statistical analyses can provide meaningful hydrological insights on the patterns between climatic indices and streamflow behavior (Sawicz et al., 2011; Addor et al., 2018; Janssen and Ameli, 2021; Anderson et al., 2022). All analyses were performed in R version 4.2.0 (R Core Team, 2022).

Spatially, SP has a generally positive and statistically significant (p < 0.05) Spearman correlation with BFI across all study catchments (ρ = 0.45), across arid in-phase catchments (ρ = 0.52), and across wet out-of-phase catchments (ρ = 0.58; Figure 2A). When scaled, the impact of SP on BFI is larger across arid in-phase catchments (ρ _scaled = 0.44), compared to wet out-of-phase catchments (ρ _scaled = 0.25; Figure 2B). Indeed, in both climatic regions, catchments with higher SP sustain larger BFI values compared to the catchments with lower SP, while such SP’s effect is stronger in arid in-phase catchments.

Further exploration at the catchment-level through the inter-annual correlational analysis, reveals that a positive association between SP and BFI does not exist temporally in many catchments. For example, in snow-rich high-elevation catchments across the Rocky Mountains and southeastern Canada, temporal associations between SP and BFI are negligible to strongly negative (Figure 4A). The standardized correlational analysis further emphasizes that in a typical snow-rich catchment, there is a negative temporal correlation between SP and BFI, with consistent negative temporal correlation during relatively dry and wet years as well as during years with relatively in-phase and out-of-phase seasonalities (Figure 5A). In a typical snow-rich catchment, the correlation coefficients between standardized SP and standardized BFI are statistically significant (p < 0.01) and slightly differ between relatively dry and relatively wet years (R = −0.31 versus −0.24) as well as between years with relatively in-phase and out-of-phase seasonalities (R = −0.29 versus −0.26).

This evidence suggests that catchments with longer SP usually have larger baseflow indices, regardless of the aridity and seasonality of the region. However, this positive spatial association cannot be extrapolated over time in most catchments, especially in high-elevation snow-rich catchments, as they usually experience a smaller BFI during years with relatively higher SP, regardless of the relative dryness and seasonality of those years.

Generally, across all study catchments, SP has a positive and statistically significant (p < 0.05) Spearman correlation with normalized Q₅ (ρ = 0.31; Figure 2A). In arid in-phase climates (ρ = 0.35) and wet out-of-phase climates (ρ = 0.38; Figure 2A), the correlation is also statistically significant. The overall impact of SP on normalized Q₅ is slightly larger across arid in-phase catchments (ρ _scaled = 0.22) than wet out-of-phase catchments (ρ _scaled = 0.12; Figure 2B). Arid out-of-phase catchments also have a positive spatial correlation between SP and Q₅, however, it is small and insignificant.

The catchment-level inter-annual correlational analysis reveals that the positive spatial association between SP and normalized Q₅ does not exist temporally in most snow-rich catchments (Figure 4B). Although there are a handful of catchments in the Pacific Northwest with positive temporal correlations, and a handful of catchments in central Canada with near zero temporal correlation, almost all other catchments have a strongly negative temporal correlation. The standardized correlational analysis further clarifies that in a typical snow-rich catchment, in a year with longer SP, a lower normalized Q₅ is more likely, especially when considering wetter years (Figure 5B). In a typical snow-rich catchment, the correlation coefficient between standardized SP and standardized normalized Q₅ is −0.29 for relatively wet years compared to −0.18 for relatively dry years. The difference in correlation is less pronounced when partitioning by seasonality index (i.e., R = −0.28 for relatively in-phase years versus R = −0.25 for relatively out-of-phase years).

The combined evidence depicts that higher SP is spatially associated with a larger extreme low-flow to average flow ratio. However, this positive spatial association cannot be extrapolated over time in almost all catchments. For example, in high-elevation snow-rich catchments across the Rocky Mountains, years with relatively higher SP exhibit a smaller normalized Q₅. Indeed, in snow-rich catchments during years with relatively higher SP, average flows may increase more than extreme low-flows.

Spatially, there are statistically significant negative (p < 0.05) Spearman correlations between SP and Low-FDC across all study catchments (ρ = −0.47), across arid in-phase catchments (ρ = −0.41), and across wet out-of-phase catchments (ρ = −0.54; Figure 2A). When scaled, the overall impact of SP on Low-FDC in arid in-phase catchments (ρ _scaled = −7.64) is larger than in wet out-of-phase catchments (ρ _scaled = −4.34; Figure 2B). This suggests that there is a general association between SP and decreased low-flow variability in North American catchments, with a particularly notable relationship between SP and decreased low-flow variability in arid in-phase catchments. While the studied relationship has the same sign for arid out-of-phase catchments, the relationship is not significant at the p < 0.05 level (ρ = −0.21).

The catchment-level inter-annual correlational analysis in snow-rich regions reveals that the temporal correlation between inter-annual variability of SP and Low-FDC does not seem to have a clear geographical pattern, with many areas having a mix of small positive and small negative temporal correlations (Figure 4C). Our standardized correlational analysis suggests that in a typical snow-rich catchment, there is a very weakly negative Pearson correlation between SP and Low-FDC, with R = −0.13 for relatively dry years and R = −0.09 for relatively wet years (and similarly for in-phase versus out-of-phase years; Figure 5C), though the correlations are still statistically significant.

These results suggest that SP is associated with mitigating low-flow variability over space, particularly in arid in-phase and wet out-of-phase regions. Our temporal findings generally agree with the spatial analysis, however, the temporal correlations are weak and appear to be insignificant and random across most individual catchments.

Spatially, there is a statistically significant positive (p < 0.05) Spearman correlation between SP and low-flow event duration across all study catchments (ρ = 0.46), across arid in-phase catchments (ρ = 0.33) and across arid out-of-phase catchments (ρ = 0.41), with weaker spatial correlation across wet out-of-phase catchments (ρ = 0.23; Figure 2A). When scaled, the overall impact of SP on low-flow event duration is the largest across arid in-phase catchments (ρ _scaled = 45.91), followed by arid out-of-phase catchments (ρ _scaled = 32.58), and wet out-of-phase catchments (ρ _scaled = 21.58; Figure 2B). This finding suggests that longer SP is linked with an increase in the average length of low-flow duration events across a wide diversity of climates in North America,.

The catchment-level inter-annual correlational analysis also reveals that the temporal Pearson correlation between inter-annual SP and low-flow event duration is strongly positive in most snow-rich catchments in North America, with some exceptions in the Pacific Northwest (Figure 4D). The standardized correlational analysis shows a statistically significant (p < 0.01) weakly positive Pearson correlation between standardized SP and standardized low-flow event duration (R = 0.16–0.23) that is slightly stronger in wetter/in-phase years (Figure 5D).

The combination of several different lines of evidence suggests that SP is associated with the lengthening of low-flow event durations, particularly in arid climates. Note that in the extraction of low-flow events, our algorithm is not limited to the extraction of only extreme low-flow events since events with flow magnitudes of consecutive days much larger than Q₅ (but still smaller than the threshold of 0.4x Q_mean) can be detected as low-flow events (see Section 2.4 for more details). In addition to a strong positive spatial association, positive temporal associative effects also exist in snow-rich catchments, albeit with some exceptions, implying a longer average length of low-flow durations in snow-rich years.

Spatially, there is a statistically significant positive (p < 0.05) Spearman correlation between SP and normalized Q₉₅ across catchments in arid out-of-phase climates (ρ = 0.51) and a negative spatial correlation in wet out-of-phase climates (ρ = −0.3; Figure 3A). When scaled, the positive overall impact of SP on normalized Q₉₅ in arid out-of-phase climates is 1.56 and the negative impact of SP on normalized Q₉₅ in wet out-of-phase climates has a ρ _scaled of −0.69 (Figure 3B). Despite the negative spatial correlation between SP and normalized Q₉₅ in arid in-phase climates, the correlation is small and not statistically significant.

The catchment-level inter-annual correlational analysis reveals that the temporal correlation between SP and normalized Q₉₅ is strongly positive in most snow-rich catchments, with the strongest positive temporal association between SP and normalized Q₉₅ across the Rockies and through parts of eastern Canada (Figure 4E). The standardized correlational analysis suggests that in a typical snow-rich catchment, there is a statistically significant (p < 0.01) positive association between SP and normalized Q₉₅ with minor variations across drier/wetter years (R = 0.35 for drier years, R = 0.33 in wetter years) and minor variations across in-phase/out-of-phase years (R = 0.35 for more in-phase years, R = 0.30 for more out-of-phase years; Figure 5E).

The results generally imply that in snow-rich catchments, years with relatively higher SP may exhibit a larger normalized Q₉₅ compared to years with relatively lower SP. Indeed, for snow-rich catchments, during years with relatively higher SP, extreme high-flow could increase to a larger extent than average flow, due to a larger contribution of surface flow (or shallow subsurface flow) to streamflow relative to the contribution of deep groundwater flow. However, over space, the associations between SP and normalized Q₉₅ could be positive or negative depending on the climatic aridity.

Spatially, there is a statistically significant negative (p < 0.05) Spearman correlation between SP and High-FDC across catchments located in arid in-phase climates (ρ = −0.29) and in wet out-of-phase climates (ρ = −0.41; Figure 3A). This suggests that SP has a role in reducing the variability of catchment high-flows in both arid in-phase and wet out-of-phase regions. When scaled, the overall impact of SP on High-FDC is slightly stronger in arid in-phase regions (ρ _scaled = −3.44) than in wet out-of-phase regions (ρ _scaled=−2.02; Figure 3B).

The catchment-level inter-annual correlational analysis reveals a strong positive temporal Pearson correlation between SP and High-FDC in most snow-rich catchments, with the strongest positive correlations appearing in the Rockies and Eastern Canada (Figure 4F). Similarly, the standardized correlational analysis suggests that in a typical snow-rich catchment, there is a statistically significant (p < 0.01) positive correlation between standardized SP and standardized High-FDC with little difference between relatively drier and wetter years (R = 0.36 for drier years, R = 0.31 in wetter years) and with no difference between relative seasonality levels (R = 0.31 for both relatively in-phase and out-of-phase years; Figure 5F).

These results suggest that a longer presence of snow on the ground from one catchment to another is associated with reduced high-flow variability, in both arid in-phase and wet out-of-phase regions. Over time, in snow-rich catchments, years with relatively higher SP may exhibit a larger variability of high-flow compared to years with relatively lower SP. Again, this might be due to a potentially higher contribution of fast and highly variable surface flow to streamflow relative to the contributions from other compartments to streamflow, in years with higher SP.

Spatially, there is a statistically significant positive (p < 0.05) Spearman correlation between SP and high-flow event duration across all study catchments (ρ = 0.49), across arid in-phase catchments (ρ = 0.31), across arid out-of-phase catchments (ρ = 0.71), and across catchments located in wet out-of-phase climates (ρ = 0.26; Figure 3A). The overall impact of SP on increasing the high-flow event duration is stronger in arid climates (for both in-phase, ρ _scaled = 17.59, and out-of-phase, ρ _scaled = 35.02, seasonalities) than in wet climates with out-of-phase seasonality (ρ _scaled = 9.04; Figure 3B).

Like the spatial analysis, the catchment-level inter-annual correlational analysis reveals a strong positive temporal Pearson correlation between SP and high-flow duration across much of North America with the strongest (Figure 4G) correlations in Western Canada and the Rockies, albeit with a few exceptions to this pattern in parts of the Pacific Northwest. Further exploration through the standardized correlational analysis depicts that in a typical snow-rich catchment, there are statistically significant (p < 0.01) positive Pearson correlations between SP and high-flow duration with only some differences between relatively drier and relatively wetter years (R = 0.31 in wetter years versus R = 0.24 in drier years) as well as between relatively in-phase and relatively out-of-phase years (R = 0.26 relatively in phase and R = 0.27 relatively out-of-phase; Figure 5G).

Generally, SP is associated with lengthening the duration of high-flow events over space and time. In particular, the spatial correlation and overall impact of SP on the lengthening of high-flow duration is stronger in arid regions than in wet regions. Note that in the extraction of high-flow events, our algorithm is not limited to the extraction of only events with extreme high-flow, since events with flow magnitudes much smaller than Q₉₅ (but still larger than the threshold of 4x Q_median) can be detected as high-flow events (see Section 2.4 for more details).

Our results suggest that SP is spatially and temporally associated with most shape-based signatures related to low-flows across North America, extending links previously explored between annual SP and annual water yield for sites in Colorado and the western United States (Hammond et al., 2018b; Harrison et al., 2021). Our inter-catchment spatial analyses suggest that a higher persistence of snow on the ground could lead to a damper (less flashy) low-flow hydrograph through larger baseflow indices, larger extreme low-flow to average flow ratios, smaller low-flow variability, and more sustained low-flow event durations (Figure 2A). These analyses further show that snow persistence becomes a more impactful factor in increasing the ratios of baseflow (or extreme low-flow) to average flow, in reducing the variability of low-flow, and in lengthening low-flow events, within arid/in-phase regions (e.g., the central United States and Canada; Figures 1, 2) compared to wet/out-of-phase regions. Generally, our results suggest that, in arid/in-phase and wet/out-of-phase climates a catchment with a larger snow persistence may have a damper low-flow hydrograph and a longer duration of low-flow events. This effect is particularly strong when comparing catchments with different degrees of snow persistence within arid/in-phase climatic regions. This is similar to links found between SP and annual water yield for dry catchments in Hammond et al. (2018b).

Over time, the magnitudes of baseflow and extreme low-flow generally increase during years with higher SP. However, years with a higher SP may show smaller ratios of baseflow (or extreme low-flow) to average flow compared to years with a lower SP (Figures 4A, B, 5A, B). A disproportionate increase in the magnitude of high-flows during snowy years may decrease the baseflow index and decrease the ratio of extreme low-flow to average streamflow while still increasing overall low-flow magnitudes. This discussion is supported by our high-flow related findings which suggest that the ratio of extreme high-flow to average streamflow and high-flow variability increases during years with higher SP in typical snow-rich catchments (Figures 4E, F, 5E, F). Our results are in agreement with prior studies that have estimated greater surface and shallow subsurface partitioning of snowmelt, as compared to deeper flow-path partitioning, during rapid snowmelt events in particularly snowy years (Barnhart et al., 2016; Fuss et al., 2016; Gordon et al., 2022; Hale et al., 2022). Indeed, a typical snow-rich catchment—which can generally be characterized with a damped and smooth hydrograph at low-flow—can have a moderate swing toward a relatively flashier low-flow hydrograph during its very snowy years.

Results also suggest that SP is, to some extent, associated with shape-based streamflow signatures related to high-flows over both space and time. A larger snow persistence is spatially associated with a lower High-FDC (i.e., damper high-flow hydrographs) in both arid/in-phase and wet/out-of-phase climates (Figure 3A). In all three regions considered in our spatial analyses, a larger snow persistence is associated with a longer duration of high-flow events with flow magnitudes continuously larger than the 4× Q_median threshold (Figure 3A). These analyses further show that snow persistence becomes a more impactful factor in decreasing high-flow variability (damping high-flow hydrograph), and in lengthening the duration of high-flow events, within arid/in-phase climatic regions (e.g., the central United States and Canada; see Figure 1) compared to wet/out-of-phase regions (Figure 3B). Generally, our results suggest that, in arid/in-phase and wet/out-of-phase climates a catchment with a larger snow persistence may have a damper high-flow hydrograph and a longer duration of events with large flow magnitudes. This effect is particularly strong when comparing catchments with different degrees of snow persistence within arid/in-phase climatic regions. On the other hand, our temporal analysis suggests that in a typical snow-rich catchment, snowier years are associated with a larger ratio of extreme high-flow to average flow and a larger high-flow variability, both reflecting a more flashy high-flow hydrograph. Indeed, in arid/in-phase and wet/out-of-phase climates, a typical snow-rich catchment—which can generally be characterized with a relatively damped high-flow hydrograph as explained above—can have a moderate swing toward a relatively flashier high-flow hydrograph during its snowier years. Please see the Supplementary Material and Supplementary Figure S2 for how SP could be linked to flow paths (and processes) that generate variability in high-flow and low-flow conditions.

Our results also suggest a stronger correlation between SP and most shape-based low-flow and high-flow signatures compared to AI and SI (Supplementary Figure S1). This simple analysis suggests that SP has stronger predictive capabilities compared to other widely used climatic indices. However, a more thorough assessment of the relative importance of these three climatic indices requires the use of a more sophisticated variable importance analysis such as those used in Addor et al. (2018), Li and Ameli (2022), or Janssen et al. (2023). Such a ranking is beyond the scope of this paper.

In sum, our findings validate that SP could be used as a new catchment descriptor for explaining shape-based low-flow and high-flow signatures over both space and time, and our findings suggest that SP could be useful globally in a variety of climate regimes. As a simple and widely available index of snow accumulation and ablation, SP can be used 1) in the development of parsimonious, large-scale, and spatially regionalized catchment hydrologic functioning frameworks and 2) in the temporal assessment of the impact of snow loss on low-flow and high-flow characteristics. Our findings also help to identify catchments with low-flow and high-flow vulnerabilities. For example, in arid/in-phase climates a catchment with a shorter persistence of snow may experience a relatively larger variability of both high-flows and low-flows and a relatively shorter duration of high-flow and low-flow events compared to a catchment with a higher persistence of snow. The former catchment would thus be less stable and more vulnerable in response to climate driven alterations. The same is true in wet/out-of-phase climates but with smaller differences among catchments in low-flow and high-flow vulnerabilities. In wet/out-of-phase climates, low and high-flows could be more dominantly sourced by rainfall, reducing the role that snow can play in maintaining sufficient and stable flow during high-flow and low-flow conditions.

Catchments in North America (and worldwide) are highly vulnerable to snow loss from climate change (Huning and AghaKouchak, 2018; IPCC, 2021; Shrestha et al., 2021; Siirila-Woodburn et al., 2021). Recent research studies suggest that within the next 35–60 years, snow-poor conditions resulting from unabated climate change could persist across large parts of western North America (Siirila-Woodburn et al., 2021). Furthermore, snowmelt timing on an annual level is expected to be earlier in the year (Musselman et al., 2021). All this paints a picture of a future with less persistence of snow on the ground if current warming trends continue. The growing aridification of large parts of North America (Overpeck and Udall, 2020) and potential changes in the annual precipitation timing (or seasonality) of large regions of the continent (Easterling et al., 2017; Gershunov et al., 2019; Zhang et al., 2019), along with the changes mentioned above in snow presence on the ground, emphasize the need for regionally-sensitive analytical tools built on the relationships among persistence of snow, aridity, seasonality, and streamflow characteristics. Such tools could help identify when and where snowpack alterations are most hydrologically influential (Siirila-Woodburn et al., 2021). Our paper takes a step in this regard by identifying the aridity and seasonality-dependent temporal associations between snow persistence and hydrograph shape (and variability) across North America.

Broadbent et al. (2021) recently showed that snowmelt is predicted to occur 50–130 days earlier in alpine regions due to climate change by 2100. Or, as Harpold and Brooks (2018) showed, Colorado snowpacks are melting as much as a month earlier than the historical norm. As a hypothetical example, an Alpine climatic catchment within North America (e.g., across the Rocky Mountains) could have a long-term snow persistence value of 0.75 (See Figure 1A). An illustrative ∼2-month decline in snow presence on the ground would imply a change in snow persistence from 0.75 to 0.40 (i.e., halving the time that snow is on the ground). While previous work might have been only able to comment that such a drop in SP would result in lower average flow (Hammond et al., 2018b; Eurich et al., 2021) and reduced peak flows (Kampf and Lefsky, 2016), we can now say that such a drop in SP can result in a substantial change in the hydrograph shape through alterations in the hydrograph’s flashiness, variability, and duration of low-flow and high-flow events.

For example, a typical snow-rich catchment such as those across the Rocky Mountains has relatively damp low-flow and high-flow hydrographs and relatively long durations of high-flow events with flow magnitudes larger than a prescribed flow threshold (see Section 4.1). Large climate-driven reductions in snow persistence for a typical snow-rich catchment may lead to significantly damped low-flow and high flow hydrographs over time and shortened durations of high-flow events (see Section 4.1 and Figure 5). This implies that climate-driven reductions in snow persistence could alter the historical functional shape and variability of streamflow hydrographs of typical snow-rich catchments. Thus, the historical functional form that surrounding ecosystems and human activities have adapted to may change drastically. Such alterations can be disruptive to ecosystems, water resources management practices and human life (Parasiewicz et al., 2019; Clarke et al., 2022; IPCC, 2022; Rodell and Li, 2023; Rohde, 2023).

In previous studies, such as with Zhou et al. (2015), it was thought that relationships from spatial analyses can be interchangeably used for temporal analyses (i.e., there is a symmetry between spatial and temporal associations). Berghuijs and Woods (2016), showed that such symmetry is not valid for assessing the association between the aridity index and the partitioning of precipitation into streamflow and evaporation (within the context of the Budyko framework). Our results also suggest that this symmetry is not valid for the linkages between snow persistence and many shape-based streamflow signatures. Although inter-catchment spatial analysis suggests a positive association between SP and BFI (and normalized Q₅; Figure 2), inter-annual temporal analysis suggests a negative association between SP and BFI (and normalized Q₅) in most snow-rich catchments (Figure 4) and for a typical snow-rich catchment (Figure 5). While inter-catchment spatial analysis suggests a moderately strong negative association between SP and Low-FDC, inter-annual analysis suggests a very weak temporal association between SP and Low-FDC in snow-rich catchments (Figures 2, 4). Such findings corroborate previous work and serve as another warning to practitioners seeking to switch the relationships across space and time. Because of the coevolution of catchment attributes (e.g., climate, topography, geology, soil, or land cover), focusing strictly on spatial analyses and extrapolating the spatial findings over time may be misleading unless the confounding effects are explicitly accounted for.

Despite the strong catchment-level temporal correlations within many snow-rich catchments (Figure 4), the correlation values in our spatial analyses and spatiotemporal “typical catchment” analyses are not strong between SP and some shape-based signatures (e.g., normalized Q₅ and Q₉₅). This was not surprising, as recent predictive models developed across a large sample of catchments also showed that shape-based signatures are very hard to predict (Addor et al., 2018). Janssen and Ameli (2021) used random forest modeling along with detailed climatic, geological, and topographical attributes to predict the shape-based signatures used in our paper across hundreds of catchments in the USA and Canada. They obtained low predictive capabilities (R ² < 0.40) when predicting shape-based signatures over space. As the goal of our paper was not prediction or identifying relative predictor importance, statistically significant p-values (Figures 2–5) allowed us to detect the linkages between SP and relevant streamflow signatures. The broad nature of our work—incorporating the climatic diversity present across North America with over a thousand catchments and almost 20,000 annual datapoints—helps to lend meaningful credence to the observed statistically significant results. Echoing previous work, large sample size studies are important for finding hydrological patterns (Sawicz et al., 2011; Addor et al., 2018; Janssen and Ameli, 2021; Anderson et al., 2022), as they guard against spurious correlations (Jackson and Somers, 1991; Bonett and Wright, 2000; Yue et al., 2002). Future work can leverage current and upcoming advances in open and accurate data, such as with catchment geology (Smelror, 2020), topography (Amatulli et al., 2018; Amatulli et al., 2020), or climatic data (Overpeck et al., 2011; Tang et al., 2021; Tang et al., 2022), and include snow persistence, along with other explanatory factors of stream low/high-flow behavior to improve the prediction of shape-based signatures.

Echoing previous work on the snow persistence product, it is known that cloud cover can disrupt snow detection (especially in winter months) of daily snowfall data by up to 22%–50% (Molotch and Meromy, 2014). Similarly, forest cover can disrupt snow detection in densely forested regions (Rittger et al., 2013). To reduce these impacts, this study used newly developed remotely sensed data that attempts to better algorithmically correct snowfall detection errors (Riggs et al., 2016). Further research and corrections may be required to create higher-quality global SP products, particularly for regions such as Asia, Europe, or South America (Hammond et al., 2018a). Such improved data could be used to evaluate the association between snow persistence and streamflow signatures at the global scale across a larger spectrum of aridity and seasonality than was considered in our paper.

There are limitations with using BFI in snow-dominated catchments. Notably, there are known issues with the estimation of baseflow values, and BFI calculations can be highly sensitive to parameterization (Voutchkova et al., 2019). This paper sought to use consistent methodologies to calculate BFIs for comparability between catchments, and we mitigated the bias in our results by using other low-flow behavior signatures (e.g., normalized Q₅, Low-FDC, and low-flow event duration). Sensitivity analyses with a variety of parameterizations can be done in the future to dig deeper into the exact relationships between SP and BFI.

Our results are limited to the associations between SP and streamflow signatures at perennial streams. We acknowledge it may be a fruitful area of future research to identify the effects that snowpack variability may have on streamflow signatures in intermittent streams across different climatic regions of North America (Eng et al., 2016). Recent research has shown that snow water equivalent was an essential descriptor of non-perennial streamflow signatures in the western United States (Hammond et al., 2021), but further research is needed to evaluate these relations in greater detail and in the context of declining snowpacks. Indeed, given that SP itself is a relatively new measure of snow behavior for hydrologic applications (Kampf and Lefsky, 2016; Hammond et al., 2018a; b;Eurich et al., 2021), such research into intermittent streams can be very insightful for inferring generalizable knowledge.

Future work can also use snow persistence, along with ecologically important metrics of catchment biological viability (McMillan, 2020), to inform catchment management efforts (Siirila-Woodburn et al., 2021) by creating thresholds, models, or even practical guidelines for levels of snow persistence required to sustain catchment ecosystems and anthropogenic activities. This can provide direct and actionable advice and insights for practitioners responsible for managing such natural resources. Another fruitful area of research could be to explore the relationships between snow persistence and hydrologic behavior in response to higher-order impacts of climate change. For example, greater wildfire risk or tree mortality in many catchments in North America (Westerling et al., 2006; Whitman et al., 2015; Abatzoglou and Williams, 2016) are likely to have effects on forest snow ablation (Gleason et al., 2013), snow persistence, and ultimately on streamflow characteristics (Siirila-Woodburn et al., 2021). Such interactions could prove informative for generating knowledge of climate change’s effects on hydrology.