Diel streamﬂow cycles suggest more sensitive snowmelt-driven streamﬂow to climate change than land surface modeling does

. Climate warming will cause mountain snowpacks to melt earlier, reducing summer streamﬂow and threaten-ing water supplies and ecosystems. Quantifying how sensitive streamﬂow timing is to climate change and where it is most sensitive remain key questions. Physically based hydrological models are often used for this purpose; however, they have embedded assumptions that translate into uncertain hydrological projections that need to be quantiﬁed and constrained to provide reliable inferences. The purpose of this study is to evaluate differences in projected end-of-century changes to streamﬂow timing between a new empirical model based on diel (daily) streamﬂow cycles and regional land surface simulations across the mountainous western USA. We develop an observational technique for de-tecting streamﬂow responses to snowmelt using diel cycles of incoming solar radiation and streamﬂow to detect when snowmelt occurs. We measure the date of the 20th percentile of snowmelt days (DOS 20 ) across 31 western USA watersheds affected by snow, as a proxy for the beginning of snowmelt-initiated streamﬂow. Historic DOS 20 varies from mid-January to late May among our sites, with warmer basins having earlier snowmelt-mediated streamﬂow. Mean annual DOS 20 strongly correlates with the dates of 25 % and 50 % annual streamﬂow volume (DOQ 25 and DOQ 50 , both R 2 = 0 . 85), suggesting that a 1 d earlier DOS 20 corresponds with a 1 d earlier DOQ 25 and 0.7 d earlier DOQ 50 . Empirical projections of future DOS 20 based on a stepwise multiple linear

Abstract. Climate warming will cause mountain snowpacks to melt earlier, reducing summer streamflow and threatening water supplies and ecosystems. Quantifying how sensitive streamflow timing is to climate change and where it is most sensitive remain key questions. Physically based hydrological models are often used for this purpose; however, they have embedded assumptions that translate into uncertain hydrological projections that need to be quantified and constrained to provide reliable inferences. The purpose of this study is to evaluate differences in projected end-ofcentury changes to streamflow timing between a new empirical model based on diel (daily) streamflow cycles and regional land surface simulations across the mountainous western USA. We develop an observational technique for detecting streamflow responses to snowmelt using diel cycles of incoming solar radiation and streamflow to detect when snowmelt occurs. We measure the date of the 20th percentile of snowmelt days (DOS 20 ) across 31 western USA watersheds affected by snow, as a proxy for the beginning of snowmelt-initiated streamflow. Historic DOS 20 varies from mid-January to late May among our sites, with warmer basins having earlier snowmelt-mediated streamflow. Mean annual DOS 20 strongly correlates with the dates of 25 % and 50 % annual streamflow volume (DOQ 25 and DOQ 50 , both R 2 = 0.85), suggesting that a 1 d earlier DOS 20 corresponds with a 1 d earlier DOQ 25 and 0.7 d earlier DOQ 50 . Empirical projections of future DOS 20 based on a stepwise multiple linear regression across sites and years under the RCP8.5 scenario for the late 21st century show that DOS 20 will occur on average 11 ± 4 d earlier per 1 • C of warming. However, DOS 20 in colder watersheds (mean November-February air temperature, T NDJF < −8 • C) is on average 70 % more sensitive to climate change than in warmer watersheds (T NDJF > 0 • C). Moreover, empirical projections of DOQ 25 and DOQ 50 based on DOS 20 are about four and two times more sensitive to climate change, respectively, than those simulated by a state-ofthe-art land surface model (NoahMP-WRF) under the same scenario. Given the importance of changes in streamflow timing for water resources, and the significant discrepancies found in projected streamflow sensitivity, snowmelt detection methods such as DOS 20 based on diel streamflow cycles may help to constrain model parameters, improve hydrological predictions, and inform process understanding. al., 2004, 2005) and challenges reservoir operations (Barnett et al., 2005;Immerzeel et al., 2020;Viviroli et al., 2011). Furthermore, ecosystems may evaporate more water as reductions in albedo increase energy inputs (Meira Neto et al., 2020;Gordon et al., 2022), decreasing runoff from upland forested watersheds (Foster et al., 2016;Jepsen et al., 2018;Milly and Dunne, 2020). More than 50 % of mountainous watersheds play essential roles in supporting downstream systems (Viviroli et al., 2007), and snowpack changes are likely to increase lowland agriculture water stress (Immerzeel et al., 2020). However, it remains difficult to predict how much streamflow timing and amount will shift in future climates (Gordon et al., 2022) due to altered snow accumulation patterns (Mote et al., 2018) and melt rates (Musselman et al., 2017) and shifts from snowfall to rainfall (Klos et al., 2014).
Physically based hydrological models are typically used to predict how snow accumulation and melt will interact with the critical zone (CZ) to affect short-term flooding and seasonal water supply (Kopp et al., 2018;Wood and Lettenmaier, 2006). In mountainous regions like the western USA, models need to accurately simulate snow processes across watersheds with varying snowpack conditions (Serreze et al., 1999) and then transport and store that water in the CZ with varying subsurface properties . More precipitation falling as rain instead of snow will result in streamflow dynamics that more closely mirror the amount and timing of rainfall. Precipitation phase (rainfall versus snowfall) is mediated by basin elevation and hypsometry (Jennings et al., 2018;Wayand et al., 2015), which also influences precipitation amounts (Houze, 2012), with higher elevations and steeper watersheds typically having higher precipitation and snowfall. Solar radiation is the primary energy source for snowmelt in snow-dominated montane watersheds (Cline, 1997;Marks and Dozier, 1992). Conversely, cloudiness lowers solar radiation and melt rates (Sumargo and Cayan, 2018). Shallower snowpacks have less cold content and begin to melt earlier when solar radiation is lower (Harpold et al., 2012;Harpold and Brooks, 2018;Musselman et al., 2017), which shifts streamflow earlier (Clow, 2010). Storage and drainage of water in the CZ control the sensitivity of streamflow to earlier rain or meltwater inputs. For example, snowmelt-mediated spring streamflow timing is more sensitive to climate change in watersheds with rapid subsurface drainage than in landscapes with deep groundwater reservoirs that drain slowly (Safeeq et al., 2013). In contrast, slowdraining watersheds have greater sensitivity to snowmeltmediated summer streamflow volume from climate change (Tague and Grant, 2009). The complexity of these storage relationships is exemplified by isotopic evidence showing that the fraction of streamflow that is young water (less than 3 months old) is smaller in steeper watersheds (Jasechko et al., 2016), suggesting that physically modeling interactions between CZ water storage and changing hydrometeorology will be challenging in mountainous areas. In a recent data-driven review, Gordon et al. (2022) proposed a predictive framework composed of three testable and interrelated mechanisms to infer changes to snowmelt-driven streamflow response under warming. Such mechanisms are associated with snow season energy and mass exchanges and the intensity of snow season liquid water input and the synchrony of energy and water availability, and their analysis highlights the complexities in predicting future streamflow in regions where multiple mechanisms interact.
Hydrologists typically apply the following two types of modeling tools to predict streamflow: empirical models and more mechanistically oriented models (conceptual or physically based land surface models). Empirical models assume that long-term and often site-to-site statistical relationships among predictor variables (e.g., precipitation and air temperature) and water fluxes (e.g., evapotranspiration and streamflow) can be used to understand and model their likely changes over time or space. Empirical models used to predict changes over time (sometimes referred to as space-fortime substitutions) have been used to predict responses to climate change in fields such as hydrology (Goulden and Bales, 2014;Jepsen et al., 2018;Sivapalan et al., 2011), biodiversity (Blois et al., 2013), and tree growth (Klesse et al., 2020). Such models use retrospective information from different places (space), typically spanning wide range of conditions (e.g., climate gradients), to predict future changes over time. For example, observed characteristics from warm regions maybe used to infer future changes in cold regions due to global warming. A limitation of this approach is that it neglects non-correlated (or independent) changes in spatially variable factors (Jepsen et al., 2018). For example, heterogeneous patterns of warming, variations in precipitation and vegetation, or changes that occur at different temporal scales (e.g., development of soil properties over hundreds to thousands of years versus shifts from rain to snow over hours) are implicitly neglected in such empirical frameworks.
Conversely, physically based models embed state-of-theart physical understanding of hydrological processes. These models typically require some degree of calibration or validation to observations (e.g., daily streamflow) to improve and assess their predictive skill. The current generation of regional weather models using the Weather Research and Forecasting model (WRF; Skamarock et al., 2008) coupled to the Noah-Multiple Parameterization land surface model (Noah-MP;Niu et al., 2011), which we refer as NoahMP-WRF, has shown promising results for modeling atmospheric and snow processes in the contiguous USA (He et al., 2019;Liu et al., 2017;Musselman et al., 2017;Scaff et al., 2020). For example, snow simulations have been used to quantify mountain snowmelt and streamflow response to climate change (Musselman et al., 2017(Musselman et al., , 2018. These simulations use a pseudo global warming approach, which perturbs the historical climate with a climate change signal from an ensemble of global climate models (GCMs); using this perturbation avoids systemic biases in the GCMs and avoids issues re-S. A. Krogh et al.: Diel streamflow cycles suggest more sensitive snowmelt-driven streamflow to climate change 3395 lated to their interannual variability . Comparisons between land surface models and empirically based predictions of future streamflow are rare but valuable (Jepsen et al., 2018) and could help to diagnose modeling deficiencies and improve predictions.
New observations of streamflow generation during snowmelt could be key to improving current hydrological models. Determining whether streamflow response was produced by rainfall or snowmelt is an important but difficult task (Weiler et al., 2018). Few simple, low-cost observational tools are available to separate rainfall-driven from snowmelt-driven contributions to streamflow or to separate this year's snowmelt from the previous years' melt and storage. One method that can be straightforwardly applied to existing long-term observations is based on coupled diel cycles in solar radiation, snowmelt, and streamflow Lundquist and Cayan, 2002). Diel (24 h) cycles in streamflow and shallow groundwater levels can result from daily cycles in snow-/ice melt and evapotranspiration, which are both ultimately driven by solar radiation inputs . This mechanistic response has been used to study watershed properties like kinematic wave celerity , the impact of snowpack variability on streamflow timing (Lundquist and Dettinger, 2005), groundwater fluctuations (Loheide and Lundquist, 2009), and transitions from snowmelt to evapotranspirationdominated streamflow fluctuations Mutzner et al., 2015;Woelber et al., 2018). More recently,  combined local observations and remote sensing to show that streamflow diel response was tightly controlled by the timing of snowpack disappearance. However, it remains unknown whether information embedded in the diel streamflow response following snowmelt events can be used to inform streamflow predictions under climate change and whether such projections are consistent with current state-of-the-art hydrological modeling. The purpose of this research is to evaluate whether land surface hydrology model simulations and a new diel streamflowbased empirical model yield similar projected end-of-century changes in streamflow volume timing across mountainous western USA headwater watersheds. To this aim, we extend the diel cycle index approach of  using diel streamflow observations to detect days when streamflow is coupled to snowmelt inputs (i.e., a snowmelt-dominated streamflow event) and investigate their contributions to historical variability in streamflow volume timing. We then compare empirical diel streamflow-based projections by the end of the century under a RCP8.5 pseudo global warming scenario against predictions from a state-of-the-art land surface model (under the same climate scenario) across 31 mountainous watersheds in the western USA to answer the following questions: 1. Do historical diel streamflow cycles indicate earlier snowmelt in warmer watersheds and years, and can we use diel observations of snowmelt to predict the timing of streamflow volume?
2. In which watersheds is the timing of snowmelt the most sensitive to climate change, as projected by an empirical diel streamflow-based model?
3. Do historical streamflow volume timings and future empirical diel streamflow-based projections diverge from commonly used, state-of-the-art land surface models?
A list with the abbreviations used in this study is presented in Table 1.

Study domain and data
We studied 31 mountainous watersheds in the western USA (Table 2), spanning snow fractions from 0.27 to 0.78 ( Fig. A3a), aridity index values from 0.22 to 2.86 (Addor et al., 2017), and soil depths from 0.27 to 2.52 m (Addor et al., 2017;Pelletier et al., 2016; Table 2). These watersheds are part of the CAMELS (Catchments Attributes and MEteorology for Large-sample Studies) dataset (Addor et al., 2017;Newman et al., 2015), which provides daily streamflow and meteorological forcing, among other observed and simulated hydrometeorological variables at the watershed scale. These watersheds were chosen because their streamflows are unregulated, they have relatively small drainage areas (< 250 km 2 ), and they are at relatively high elevations (> 1000 m a.s.l. -above sea level). This last criterion was introduced to focus on watersheds with snowmeltdriven streamflow regimes. The names, locations, elevations, slopes, drainage areas, and other key characteristics of the 31 watersheds are presented in Table 2. The data used in this analysis include hourly streamflow, incoming shortwave radiation, mean daily relative humidity, air temperature, and precipitation. Hourly streamflow was obtained from the U.S. Geological Survey (USGS). Hourly incoming shortwave radiation is from phase 2 of the National Land Data Assimilation System (NLDAS-2; Xia et al., 2012) at the nearest grid point to the watershed outlet. Mean daily relative humidity, air temperature, and precipitation at the watershed scale are from CAMELS, based on the Daymet dataset (https://daymet.ornl.gov/, last access: 20 June 2022), which in turn is interpolated from existing ground observations. Available hourly streamflow records vary significantly across watersheds, extending back to 1986 for some sites. Figure A1a shows the number of years that have more than 70 %, 80 %, and 90 % of days with hourly records for the period between 1 December and 1 August. Based on this preliminary analysis, we selected water years with more than 80 % of days with hourly streamflow records. This threshold for data availability results in most watersheds having more than 5 years to analyze (except for site ID 10 and 30 with 4 years).

Snowmelt and streamflow diel coupling
To identify days when solar radiation-driven snowmelt is coupled to the streamflow response, hereafter called snowmelt days for simplicity, we calculated the correlation between hourly values of solar radiation and lagged stream- flow ( Fig. 1). A snowmelt day is defined as a day in which the Spearman correlation between hourly solar radiation and lagged streamflow is statistically significant (p value ≤ 0.01) and exceeds a given cutoff. Due to the lagged diel streamflow response after snowmelt, we lagged diel streamflow from solar radiation between 6 and 18 h, computed the correlation of all combinations, and kept those statistically significant correlations that were above a predefined correlation cutoff. Although having both a correlation cutoff and a statistical significance criterion may be redundant, we used both to guarantee significant correlations above different correlation cutoffs. We tried several correlation cutoffs (r > 0.5, 0.6, 0.7, 0.8, and 0.9; see Fig. 1 for r > 0.6) to assess their effects on the detection algorithm (Fig. A2). The preliminary lag window of 6 to 18 h was used to avoid confounding snowmelt signals with evapotranspiration (ET)-induced streamflow diel responses Mutzner et al., 2015;Woelber et al., 2018). ET-induced streamflow diel response can positively correlate with solar radiation with lags below 6 h, due to the previous day's ET, and above 18 h, due to the next day's ET diurnal signal . However, this preliminary lag window may incorrectly select days with a rainfall-induced streamflow diel response or rain-on-snow events. To minimize this, we further restricted the lags that could be selected based on optimum lags from snowmelt days with clear skies. Clear-sky days were defined as days with solar radiation greater than 80 % of the clear-sky solar radiation value (gray areas in left panels on Fig. 1). This lag window was defined on a monthly and watershed basis and was calculated as the lags between the 10th and 90th percentile of clear-sky days with Spearman correlations above 0.8. This second filter also helped to avoid the incorrect selection of ET-induced streamflow diel response, as it minimized the chance of selecting 18 h lags that can be associated with ET. Despite efforts to select only snowmelt-driven streamflow diel responses, this methodology does not guarantee that rainfall-driven streamflow diel changes with lags within our lag window will always be excluded. Excluding such cases would require hourly precipitation observations, which are unavailable at some of our study watersheds. However, we believe that any such cases will minimally affect the results of our analysis.
To better assess the potential impact that rainfall may have on our proposed diel analysis, particularly on the effect of rain-on-snow events, we analyzed which days classified as snowmelt days also had rainfall. We assessed daily rainfall using the daily precipitation time series from CAMELS based on the Daymet product for each watershed. A false detection rate metric was computed for each watershed, in which every day classified as a snowmelt day with daily precipitation above 5 mm and a mean daily air temperature above 2 • C was assumed to be misclassified (Fig. 2). A false detection rate of 100 % means that all snowmelt days were misclassified and 0 % means that no days had significant rainfall. On average, the false detection rate was estimated at 7 %, with a standard deviation of 5 %, and only watershed Figure 2. Percentage of days that were classified as having snowmelt following the diel streamflow cycle analysis that also had daily precipitation above 5 mm and a mean daily air temperature above 2 • C. Symbols are associated with the mean annual percentage of snowmelt days under clear-sky conditions. Sunny sites (circles) have > 90 %, clear-sky snowmelt days, partly cloudy sites (squares) have between 70 % and 90 %, and cloudy sites (diamonds) have < 70 % clear-sky snowmelt days. Clear-sky snowmelt days are defined as those with more than 80 % of the potential clear-sky solar radiation.
ID 24 and 31 (located in WA and OR, respectively) exceeded 15 %, with 21 % and 29 %, respectively. This suggests that the effect of potential rainfall-induced diel streamflow cycles (including rain-on-snow events) in most watersheds is low (except for watershed ID 24 and 31), supporting further analysis. We also assessed the mean cross-site false detection rate for precipitation thresholds of 1 and 10 mm and found reasonable values of 12 % and 3 %, respectively. However, we believe that 1 mm is not a reasonable threshold, as a 1 mm rainfall event would be unlikely to produce a distinguishable diel streamflow signal and could represent error/noise in the Daymet product.

The empirical diel streamflow-based model
We defined the day when the 20th percentile of the snowmelt days (as defined in Sect. 2.2) occurs (DOS 20 ) as a new metric to characterize the seasonality of early snowmelt for each water year and watershed. However, other metrics such as the 5th, 10th, and 30th percentiles (presented in the Appendix) were also investigated to assess the impact of this choice on the analysis. We chose this metric because we expected it to be associated with the timing of streamflow volume, and the choice of slightly earlier or later snowmelt day metrics (e.g., DOS 10 or DOS 30 ) would not substantially change our results. We fitted a stepwise multiple linear regression model (MLR; p value < 0.01; Eq. 1) to reconstruct historical DOS 20 across all watersheds and years ( Fig. 7) using four climate variables as predictors, i.e., total precipitation, air temperature, relative humidity, and solar radiation, as follows: where x 1 is cumulative air temperature (i.e., degree day; • C), x 2 is cumulative precipitation (mm), x 3 is mean relative humidity (%), x 4 is mean solar radiation (W m −2 ), and the β i are regression coefficients. Mean annual climate variables were calculated for the period between 1 November and DOS 20 (i.e., between late fall and the metric representing the date of early snowmelt events). As a result, DOS 20 is present in both sides of Eq. (1); therefore, the stepwise MLR requires an iterative solution when used in a predictive mode (i.e., for the climate change analysis when DOS 20 is unknown). The MLR model is the basis of our empirical diel streamflowbased model which is used to assess changes in DOS 20 due to climate change (i.e., changes in x 1 , x 2 , x 3 , and x 4 in Eq. 1). We verified the stepwise MLR assumptions, namely linear relationships between each predictor and DOS 20 , normally distributed residuals, homoscedasticity, and the absence of strong multicollinearity (as suggested by a variance inflation factor < 3). We also tested other metrics related to the timing of early snowmelt events. These included the first snowmelt day, the first three consecutive snowmelt events, and the 5th, 10th, and 30th percentiles of snowmelt days (DOS 5 , DOS 10 , and DOS 30 , respectively). All metrics were also computed using each of the different Spearman correlation cutoffs (Tables A1-A5), but the main analysis presented here focuses on DOS 20 based on snowmelt days calculated with hourly Spearman correlations > 0.8. We predict changes to DOS 20 based on the stepwise MLR model and end-of-the-century mean climate change forcing from NoahMP-WRF . NoahMP-WRF was run under a high-emission scenario (RCP8.5), using the pseudo global warming approach for the end of the century. Overall, it projects a warmer (4-5.2 • C), wetter (0 %-20 % increase in precipitation) climate (Figs. A4 and A5). These mean annual changes in climate were applied to the predictors in the stepwise MLR model to predict changes in DOS 20 . As previously mentioned, predictors used in the stepwise MLR were calculated for the period between 1 November and DOS 20 ; therefore, as we do not know the value of DOS 20 in the future, an iterative solution is required to solve for DOS 20 in Eq. (1). We find a numerical solution, using a 2 d convergence threshold between iterations, so that |DOS20 i+1 -DOS20 i | ≤ 2 d, where i is the number of the iteration.

Streamflow volume timing from a land surface model
Historical NoahMP-WRF simulations include the period 2001-2013 over the contiguous USA at 4 km spatial resolution and the period 2071-2100 under pseudo global warming S. A. Krogh et al.: Diel streamflow cycles suggest more sensitive snowmelt-driven streamflow to climate change 3399 . NoahMP-WRF simulations include an improved Noah configuration, which aims to better represent the snow physics. These improvements include the following : the rain-snow transition is based on a microphysics partitioning approach as opposed to a subjective temperature-based approach, patchy snowpack are allowed in the calculation of the surface energy balance, the heat transport from rainfall to the ground is included, and the snow depletion curve is vegetation dependent. These improvements allow for a better representation of the surface energy balance and the simulation of snow accumulation and melt processes. We used daily watershed-scale outputs of surface and subsurface runoff from historical and future NoahMP-WRF simulations to estimate the date of 25 % and 50 % of annual streamflow volume (DOQ 25 and DOQ 50 , respectively). Given the range of the watershed drainage areas (4-236 km 2 ; Table 2), watersheds covering several grid cells use the total surface and subsurface runoff for their corresponding grid cells. Small watersheds are represented by only the single nearest NoahMP-WRF grid cell. The way NoahMP-WRF is implemented within WRF lacks a streamflow routing scheme such as the one in WRF-Hydro (Gochis et al., 2020); therefore, we used the sum of surface and subsurface runoff to estimate DOQ 25 and DOQ 50 . We also repeated the analysis using surface runoff only, leading to similar results (Fig. A7). Given the relatively coarse NoahMP-WRF spatial resolution (4 km) compared to the watershed drainage areas (4-236 km 2 ), we assume that mean streamflow timing metrics are not significantly affected by the lack of streamflow routing.

Empirical relationships between DOS 20 , climate, and streamflow
Mean annual DOS 20 (the date of the 20th percentile of snowmelt days) has a strong regional variability that is reasonably captured by a negative linear correlation (R 2 = 0.48) with the mean winter air temperature (November to February; T NDJF ) in watersheds with T NDJF < −3 • C, whereas warmer watersheds do not follow the same pattern (Figs. 3a and 4a). Warmer sites (T NDJF > −3 • C) have a more variable mean DOS 20 , ranging from mid-January to early May, whereas the coldest sites (T NDJF < −8 • C) have a later and less variable DOS 20 around mid-to late May. On average, the regression suggests that a 1 • C warming of results in 7.2 d earlier DOS 20 . A relationship between later DOS 20 and colder T NDJF is also found in the year-to-year variations at most watersheds (21 out of 31; Fig. 3b). A strong negative linear relationship was found between the date of the 25 % of the annual streamflow volume (DOQ 25 ) and T NDJF (Fig. 3c). Warmer watersheds (T NDJF > 0 • C) generate streamflow earlier (DOQ 25 between mid-December and early March) com-pared to the coldest watersheds (T NDJF < −8 • C) where DOQ 25 is between early and late May (Fig. 3c). On average, the cross-site regression shows that each 1 • C warmer T NDJF produces a 13 d earlier DOQ 25 . For most watersheds (25 out of 31), interannual regressions show a similar pattern, with warmer years having earlier DOQ 25 ; however, these interannual regressions have shallower slopes than the cross-site relationships ( Fig. 3b and d). Previous work by Stewart et al. (2005) also related seasonal meteorological patterns with the spring onset and streamflow timing and found similar relationships (e.g., warmer watersheds have earlier spring onset and streamflow timing). However, the definition of the spring onset was based on the cumulative hydrograph (the day when the cumulative departure from the mean streamflow was the minimum), as opposed to our more mechanistic diel streamflow analysis. Other definitions for spring onset based on streamflow, snow pillows, and air temperature are presented by Lundquist et al. (2004).
Strong correlations between DOS 20 and both DOQ 25 and DOQ 50 (the date of 50 % of the annual streamflow volume; R 2 = 0.85; Fig. 5a and c) suggest connections between the timing of snowmelt and streamflow generation across watersheds and years. On average, sites that melt earlier are associated with earlier DOQ 25 (Fig. 5a) and a lower ratio of snowfall to total precipitation (snow fraction < 0.5). The relationship between DOS 20 and DOQ 25 closely follows the 1 : 1 line (Fig. 5a), although three sites in Washington and Oregon (site ID 24, 25, and 31; see Table 2 and Fig. 6a) deviate substantially from this pattern, perhaps because they receive relatively little of their precipitation as snow. Similar watershed-level relationships using interannual variability in DOQ 25 were found for most watersheds, with statistically significant slopes varying between 0.4 and 2.5 d d −1 (Fig. 5b). DOS 20 also predicts DOQ 50 well, with 10 d earlier snowmelt producing 7 d earlier DOQ 50 on average (Fig. 5c) and similar watershed-level interannual relationships (Fig. 5d). The same three relatively rainy watersheds have DOQ 50 prior to the DOS 20 (Figs. 5c and 6b), suggesting that early snowmelt timing is not an important predictor of DOQ 50 in such places.

Diel streamflow-based sensitivity of snowmelt timing (DOS 20 ) to climate change
We fitted a stepwise MLR with four climate variables (air temperature, precipitation, relative humidity, and solar radiation) to predict the diel streamflow-based DOS 20 metric across watersheds and years. A total of 333 watershed-year combinations of DOS 20 and climate variables were used to train the stepwise MLR model. The watershed-year relationship between observed and MLR predictions has a relatively high R 2 of 0.83, a root mean square error (RMSE) of 17.5 d, and normally distributed residuals (p < 0.01) off the predictions of inter-watershed mean annual DOS 20 (Fig. 7b) is also strong (R 2 = 0.83 and RMSE = 13.2 d) and follows the 1 : 1 line. Similarly, when we look at interannual values, represented by the lines overlapping the circles in Fig. 7b, we find a good agreement with most slopes close to 1 : 1 (see the inset in Fig. 7b). This analysis demonstrates that the MLR model can reasonably represent both the mean annual DOS 20 values at each watershed and their interannual variability. Table A4 shows the standardized beta coefficients that indicate the importance of each climate variable in the stepwise MLR. For the 0.8 correlation cutoff, we found that incoming shortwave radiation has the greatest importance (β = 0.75), followed by relative humidity (β = 0.37) and air temperature (β = −0.31).
Empirical diel streamflow-based projections under climate change show earlier mean annual DOS 20 in all watersheds (i.e., earlier snowmelt initiation), with significant variability from site to site (Fig. 8a). Most watersheds show significant end-of-century changes in DOS 20 , ranging from up to 3 months earlier in cold sites where, historically, snowmelt under clear-sky conditions dominates (circles in Fig. 8a) to as little as 20 d earlier in warm sites under historically cloudier conditions. The cross-site average change in DOS 20 is 55.3 d, with a standard deviation of 21.8 d. In many watersheds, the mean projection of DOS 20 under climate change is within the historically observed variability in DOS 20 (Fig. 8a). The empirical model predicts that, on average, colder watersheds (T NDJF ≤ −8 • C) are about 70 % more sensitive to climate change (13.7±4.6 d • C −1 ) than warmer watersheds (T NDJF > 0 • C) (8.1 ± 6.2 d • C −1 ), as represented by the change in the DOS 20 per degree of warming (Fig. 8b). Site ID 24 (South Fork Tolt River, WA) shows almost no change in its DOS 20 , which can be attributed to its weaker climate change signal compared to the other watersheds (about +4 • C, 5 % precipitation increase and virtually no change in humidity and solar radiation; Fig. A4). The diel streamflow-based analysis suggests an average sensitivity of DOS 20 to a climate change of 11.1 ± 4.2 d • C −1 across all watersheds.   ; Fig. 9c), whereas empirical diel streamflowbased projections range between early January and late March (mean in mid-February; Fig. 9a), averaging about −60 d ( DOQ 25 ; Fig. 9c). These results indicate that empirical diel streamflow-based projections of DOQ 25 are about 4 times more sensitive to climate change than those from NoahMP-WRF. Historical DOQ 50 is reasonably well represented by NoahMP-WRF under the current climate (blue symbols; Fig. 9b) with a mean difference against observations of 7 d; however, future changes of about −20 d are projected, which are roughly half of the −40 d predicted by the empirical streamflow-based projections ( DOQ 50 ; Fig. 9d). Empirical diel streamflow-based projections of DOQ 50 range between mid-February and early April, whereas NoahMP-WRF projections range between mid-March and mid-May, suggesting later estimates of streamflow volume by the land surface model. Watersheds with the largest disagreement between the empirical model and NoahMP-WRF projections for streamflow volume timing are those where DOS 20 is the most sensitive to warming (represented by the orange and yellow symbols in Fig. 9c and d). These watersheds are characterized by historical cold winter temperatures (T NDJF < −6 • C), with snowmelt occurring mostly under sunny conditions (circle symbols) in the Rocky Mountains.

Discussion
The new DOS 20 metric based on the diel streamflow analysis quantifies the timing of early snowmelt events and suggests that shifts towards earlier snowmelt will generate larger shifts toward earlier streamflow in colder, sunnier watersheds than in warmer, cloudier watersheds where snowmelt is more interspersed with rain. Despite the intuitive connections between snowmelt and streamflow, empirically linking changes in earlier snowmelt rates (Harpold and Brooks, 2018;Musselman et al., 2017) with changes in streamflow amount (Barnhart et al., 2016) and timing (Stewart et al., 2004) has been challenging (Weiler et al., 2018). This study represents of the first empirical analysis of streamflow-induced snowmelt change across a regional climate gradient not relying only on streamflow volume. Understanding these connections is challenging due to the representative scales at which snow (point scale) and streamflow (watershed-scale) are typically measured and analyzed. For example, evidence of snowmelt at Snow Telemetry (SNOTEL) sites in the USA has shown more intermittent snowmelt events at sites with higher humidity, and future modeling suggests lowerhumidity sites will experience slower, earlier snowmelt (Harpold and Brooks, 2018; Musselman et al., 2017). However, the cascading effects of earlier and slower snowmelt on streamflow amount and timing remain relatively unexplored (e.g. Berghuijs et al., 2014) and are potentially affected by surface and subsurface hydrological connectivity, vegetation water use, and other processes that are not easily measured or parameterized. Our diel streamflow analysis has limitations in places dominated by rainfall, as evidenced by higher false detections in areas with low snow fractions (Fig. 2) and by the small (or nonexistent) interannual correlation between DOS 20 and the metrics DOQ 25 and DOQ 50 (Fig. 5a and c) in those places. Conversely, the colder and sunnier watersheds, primarily in the intermountain region, have strong interannual correlations between DOS 20 and DOQ 25 (Figs. 5a and 6a), reflecting the importance of snowmelt (instead of rain) in controlling streamflow volume timing. Because the diel streamflow analysis does not require the many assumptions that are embedded in physically based models, it is an independent tool that can be used to verify historical streamflow simulations from subdaily resolved hydrological models. For example, land surface models could be benchmarked against observed snowmelt days based on the diel streamflow analysis or metrics like DOS 20 to better represent processes associated with snowmelt-driven streamflow generation. The diel streamflow analysis is also easier to implement than detailed process-based models because it only requires observed hourly streamflow data and solar radiation. If measured solar radiation is not available, it can be reliably represented by land surface models like NLDAS-2 (Luo et al., 2003) that assimilate field observations and remotely sensed radiation (including the effects of clouds) into an atmospheric modeling framework. In our analysis, we tested the sensitivity of some modeling decisions, such  as the correlation cutoff between hourly solar radiation and streamflow used to detect snowmelt days and metrics for snowmelt timing and found similar sensitivities of DOS 20 to climate change across different correlation cutoffs and snowmelt timing percentiles (Table A5). Metrics like the first snowmelt day or the first 3 consecutive snowmelt days showed less consistent results (Table A5), likely due to individual early or midwinter melt events that do not necessarily represent the seasonal watershed behavior. The diel streamflow analysis has the following four main limitations that need to be examined in future work: (1) it requires a steep enough stage-discharge relationship so that daily stream-flow cycles can be detected across the flow regime, (2) it focuses on snowmelt driven by solar radiation (and energy fluxes synchronized with it), (3) it is sensitive to assumptions about the lag time between solar radiation and streamflow, and (4) it is sensitive to assumptions about evapotranspiration losses. A steep stage-discharge relationship, in which small changes in discharge are associated with large changes in stage, is ideal to observe small diel streamflow changes with sufficient precision. The second limitation originates from the assumption that the majority of snowmelt is correlated with solar radiation, which is supported by the dominant role of solar radiation in process-based studies of mar- itime and continental snowpacks (Cline, 1997;Jepsen et al., 2012;Marks and Dozier, 1992). Because our method allows the lag time between solar radiation and streamflow to vary within a predefined window, we expect it to capture the effects of other important energy fluxes, such as sensible heat, that often lag the diel patterns of solar radiation by several hours (Ohmura, 2001). Rain-on-snow events are particularly challenging to detect with our analysis, as days with a lower percentage of incoming shortwave radiation (< 80 % of clear sky) are filtered out to avoid issues with potential rainfall-dominated diel signals. It may also misclassify rainfall-driven diel streamflow cycles, although we checked for rainfall-induced cycles and found that these accounted for only a small fraction (7 % on average; Fig. 2) of our inferred snowmelt days. The relationships between streamflow timing (i.e., DOS 20 , DOQ 25 , and DOQ 50 ) and meteorological drivers in rainier sites showed cross-site and interannual relationships that are consistent with those in colder, more snow-dominated places (except for watershed ID 24, 25, and 31; e.g., Fig. 3a and c). The third limitation is that the spatiotemporal variability in snowpack, surface and subsurface storage, and evapotranspiration will change the magnitude and lag time of the diel streamflow response Lundquist and Cayan, 2002;Lundquist and Dettinger, 2005), which we address by allowing variable watershedand month-specific time lags. However, lag times greater than 24 h, which are associated with large watersheds or large subsurface storage, will make this method impossible to apply. The method may also miss early snowmelt-driven diel cycles in watersheds with dry soils, as the diel signal will be buffered by the subsurface storage capacity before generating a measurable streamflow response. Our empirical diel streamflow-based model implicitly assumes that other variables not included in the analysis vary together with the pre-dictive variables (climate) and neglects watersheds' physical (e.g., soil storage) and biological (e.g., vegetation) properties that do not necessarily covary with climate. The fourth limitation is that evapotranspiration losses must be small relative to snowmelt inputs, which is necessary because the effect of evapotranspiration is out of phase with the effect of snowmelt . Evapotranspiration effects are minimized by focusing on early snowmelt periods, when evapotranspiration losses are small (Bowling et al., 2018;Cooper et al., 2020;Winchell et al., 2016).
Hydrological modeling in land surface models attempts to physically represent snowpack storage, snowmelt, subsurface storage, and its release to the streamflow, which is challenged by uncertain forcing data and simplified and uncertain model parameters. For example, snowmelt modeling in complex terrain is challenged by steep climate gradients and by the lack of adequate forcing data (precipitation, temperature, wind, etc.). Characterizing precipitation phase and timing in steep watersheds remains challenging in rainto-snow transition zones (Harpold et al., 2017;Jennings et al., 2018;Wayand et al., 2015), which will presumably increase in extent in the future (Klos et al., 2014). Complex terrain affects radiation fluxes, which are hard to estimate at kilometer spatial scales (Müller and Scherer, 2005) used in most land surface models. Most of our study sites are forest covered, which exerts a strong control on the snowpack mass and energy balance (Lundquist et al., 2013;Pomeroy et al., 1998;Safa et al., 2021) with spatially heterogeneous effects on snow accumulation and melt that remain challenging to model (Broxton et al., 2015;Krogh et al., 2020). The presence of preferential flow paths through the snowpack impacts the timing of melt release (Leroux and Pomeroy, 2017) and is not typically included in hydrological models. Once snowmelt is released from the snowpack, simulating Figure 9. Changes to DOQ 25 and DOQ 50 due to climate change under a RCP8.5 pseudo global warming climate scenario by the end of the century. Panels (a) and (b) compare historical against projected values between NoahMP-WRF and the empirical diel streamflow-based model. Panels (c) and (d) compare the projected change in streamflow timing (future minus historical) between NoahMP-WRF and the empirical diel streamflow-based model, colored by the sensitivity of DOS 20 to climate change as projected by the empirical diel streamflowbased model (Fig. 8b). Symbols surrounded by black circles indicate sites that were excluded from the regression analysis in Fig. 5 (rainier site ID 24, 25, and 31). Symbols represent the historical mean annual percentage of clear-sky snowmelt days, where sunny sites (circles) have > 90 % clear-sky snowmelt days, partly cloudy sites (squares) have between 70 % and 90 %, and cloudy sites (diamonds) have < 70 %; clear-sky snowmelt days are defined as those with more than 80 % of the potential clear-sky solar radiation. We make no inference about the cloudiness condition of snowmelt days under the RCP8.5 climate scenario; however, red symbols (upper panels) follow the same symbology for easier interpretation.
(and validating) what fraction flows as subsurface and surface runoff remains difficult. Decades of tracer studies (e.g., Godsey et al., 2010;Kirchner, 2003) have shown that streamflow during and after hydrologic events (i.e., snowmelt or rainfall events) is typically old water that has been stored in the watershed for months to years. Land surface models like NoahMP-WRF lack realistic groundwater stores to represent old water and lack hillslope and near-stream processes (Fan et al., 2019). For example, previous work at Sagehen Creek (site ID 23) suggests that streamflow remains ∼ 80 % groundwater even during the snowmelt freshet (Urióstegui et al., 2017), despite a strong snowmelt diel response caused by pressure changes induced by infiltrating snowmelt. Innovative observations that give new physical insights, like the diel streamflow analysis, could bring new information to modeling beyond what is possible with typical daily discharge resolution (Kirchner, 2006).
The diel-based analysis of snowmelt-driven streamflow to changing climate gives unique insights over previous efforts using daily and seasonal streamflow volumes (Berghuijs et al., 2014;Stewart et al., 2005) and retrospective hydrological modeling (Barnhart et al., 2016). Empirical projections of DOS 20 under the pseudo global warming scenario (Fig. 8b) show that colder, drier, and sunnier sites (typical of the Rocky Mountains) are about twice as sensitive to warming as warmer, more humid, and cloudier sites (typical of the Pacific Northwest). Humid and warmer sites have lower snow fractions (< 0.5; more rainfall effects) and, thus, a smaller snowmelt signal in the diel streamflow observations. In contrast, Harpold and Brooks (2018) showed that winter ablation at SNOTEL sites in humid places, like the Pacific Northwest, are more sensitive to warming than less humid places, like the southwestern USA. However,  showed general agreement between SNOTEL snowmelt response and the snowmelt-induced diel streamflow signal at the warm Sagehen Creek watershed (site ID 23). The sensitivity of the early snowmelt timing metric (DOS 20 ) to climate change is a function of changes in precipitation phase (rainfall versus snowfall), snowpack ablation (changes in the patterns of melt and sublimation), and hydrological partitioning to streamflow versus evaporative loss. Due to the empirical basis of our analysis, these sensitivities are not easy to disentangle, but the diel analysis is a new source of information that could help in that effort. The reliability of the empirical diel streamflow-based projections partially depends on whether climate projections are within or outside the range of observed climate conditions across the large climatic gradient found in the western USA. Under the pseudo global warming scenario, cold, sunny watersheds like those in the Rocky Mountains (site ID 9 and 10) will shift toward more humid, warmer conditions (Fig. A6), like those observed in Southern Idaho (site ID 29) and the northern Sierra Nevada (site ID 23). In contrast, the pseudo global warming scenarios for places like the Pacific Northwest, particularly those involving changes in atmospheric humidity above 5 g m −3 (Fig. A4), have not been observed in the historical record and therefore are more uncertain. Determining reasonable conditions to apply empirical models that use observed differences in sites to predict future changes (often called spacefor-time models), like the presented diel streamflow analysis, has been posed as one of the 23 unsolved problems in hydrology (Blöschl et al., 2019).
The sensitivity of historical snowmelt-mediated streamflow volume timing (DOQ 25 and DOQ 50 ) to climate change differs substantially between the empirical diel streamflowbased approach and a land surface model, raising questions about current state-of-the-art projections of early season streamflow timing from NoahMP-WRF, particularly in cold watersheds ( Fig. 9c and d). The observed data used in the diel streamflow-based approach have larger and more variable streamflow timing responses to climate change (10-17 d • C −1 ) in cold, dry, sunny places that are representative of small, high-elevation Rocky Mountain watersheds (Fig. 8b). The historical diel streamflow analysis suggests that NoahMP-WRF may be systematically underpredicting the sensitivity of streamflow volume timing to earlier snowmelt-induced streamflow in colder and sunnier places (Fig. 9c) that are most likely to have increased temperature and increased cloudiness in the future. The same mean annual future climate scenarios were applied to both approaches; however, important differences in the streamflow timing response were found between NoahMP-WRF and diel streamflow-based projections ( Fig. 9c and d). NoahMP-WRF underpredicts historical DOQ 25 (Fig. 9a) across most sites, whereas DOQ 50 is much better represented. It is worth noting that when DOQ 25 simulated by NoahMP-WRF is calculated using surface runoff alone (Fig. A7a), rather than subsurface plus surface runoff, it performs better against observed DOQ 25 . However, NoahMP-WRF projected sensitivity in streamflow timing to climate change remains significantly lower than predictions based on the diel-streamflow analysis (Fig. A7c). We used these simulations in the analysis because Noah-MP underlies the U.S. National Water Model, and thus, its relevance to policy and research is high. There are many differences in the way that NoahMP-WRF and the empirical diel streamflow-based approach function. NoahMP-WRF can track the hourly covariance in precipitation, temperature, and humidity to estimate precipitation partitioning between rain and snow. It is also able to represent hourly radiative and turbulent energy at the snowpack, and the cold content needed to predict snowmelt. Its physical hydrology is also advanced and able to consider antecedent conditions and allow evapotranspiration losses that also modulate streamflow. Despite the advantages of land surface models like NoahMP-WRF in constraining processes for future projections, the simplicity of diel streamflow-based analysis also provides several advantages. One of the main advantages is that it is derived from observations, and thus, it is well constrained by the observed spatial and temporal variability of snowmelt across watersheds and years (Fig. 7b). Also, it does not assume anything about the complex spatial distribution of snowpacks and precipitation or subsurface properties, which are major constraints to physically based models (Baroni et al., 2010;Christiaens and Feyen, 2001;Wilby et al., 2002). While the empirical diel streamflowbased model is not a replacement for land surface models like NoahMP-WRF, partly because the underlying streamflow datasets are not available everywhere, there is added value in including new benchmarks like the proposed DOS 20 to further constrain modeling decisions and improve model fidelity required for reliable and accurate hydrological predictions.

Conclusions
Water management in the western USA requires accurate predictions of how both short-term climate variability and long-term climate change will alter snowmelt and streamflow. Differences in predictions of snowmelt-induced streamflow between empirical diel streamflow-based projections and a land surface model (NoahMP-WRF) raise important questions about the sensitivity of streamflow timing to climate change, particularly in cold regions, and its impact on water planning. Significant differences exist in the way diel streamflow-based and land surface models predict changes to snowmelt and streamflow timing, with both approaches having strengths and weaknesses; however, the land sur-face model misrepresents historical patterns in streamflow response that are more accurately estimated by the empirical model. We show that DOS 20 is a strong predictor of the early season hydrograph response, particularly in cold, sunny areas where the NoahMP-WRF streamflow timing simulations lack sensitivity to climate change. Rigorously validating future model predictions is impossible, but snowmelt and streamflow timing, inferred from diel streamflow cycles, could be used to refine land surface models and better determine the risk to valuable snow water resources (Barnett et al., 2005;Sturm et al., 2017;Viviroli et al., 2007), particularly in cold regions. Our novel approach can complement the benchmarking or calibration of physically based hydrological models beyond typical benchmarking against daily streamflow or snow accumulation metrics. For example, the snowmelt timing metric DOS 20 based on diel streamflow observations could be used to test how well land surface models, running at subdaily scales and fine spatial resolution, can reproduce the historical snowmelt regime across watersheds and years. As land surface models move towards real application for water management (Kopp et al., 2018), the hydrology community must seek ways to test and improve them using widely available datasets if we are to meet the grand water management challenges posed by climate change in mountainous regions.   In red are the perturbed mean climate variables under the RCP8.5 pseudo global warming scenario by the end of the century. This analysis suggests that most of the climate change signal from NoahMP-WRF pseudo global warming is within the observed climate variability, except for air temperature and atmospheric humidity in some watersheds. Blue symbols (circle, square, and diamond) associated with historical values represent the mean annual percentage of clear-sky snowmelt days, where sunny sites have > 90 % clear-sky snowmelt days, partly cloudy have between 70 % and 90 %, and cloudy have < 70 %; clear-sky snowmelt days are defined as those with more than 80 % of the potential clear-sky solar radiation. We make no inference about the cloudiness condition of snowmelt days under the RCP8.5 pseudo global warming scenario, and thus, we use a five-point star (in red) for the future scenario. Figure A5. Mean annual climate changes projected by WRF under a RCP8.5 pseudo global warming scenario by the end of the century. Panel (a) shows changes in precipitation against air temperature. Panel (b) shows incoming shortwave against absolute humidity. Numbers represent the gauge IDs as presented in Table 2. Figure A6. (a) Principal component analysis for historical precipitation (Pp), air temperature (AT), absolute humidity (AH) and shortwave radiation (SWR) at each watershed and the changes associated with the pseudo global warming as simulated by WRF. Panel (b) shows the same analysis but excluding precipitation from the analysis. Blue symbols (circle, square, and diamond) associated with historical values represent the mean annual percentage of clear-sky snowmelt days, where sunny sites have > 90 % clear-sky snowmelt days, partly cloudy have between 70 % and 90 %, and cloudy have < 70 %; clear-sky snowmelt days are defined as those with more than 80 % of the potential clear-sky solar radiation. We make no inference about the cloudiness condition during snowmelt days under the RCP8.5 pseudo global warming scenario, and thus, we use a five-point star (in red) for the future scenario. Numbers next to blue symbols represent the gauge IDs as presented in Table 2. Figure A7. Same as Fig. 9 but using streamflow timing metrics from NoahMP-WRF for a RCP8.5 pseudo global warming scenario, calculated using surface runoff only instead of using surface plus subsurface runoff (as in Fig. 6). Note the improved fit in historical DOQ 25 ; however, this analysis yields very similar results to those in Fig. 6, with NoahMP-WRF streamflow simulations being much less sensitive to climate change than the empirical diel streamflow-based model suggests. Panels (a) and (b) compare historical against projected values between NoahMP-WRF and the empirical diel streamflow-based model. Panels (c) and (d) compare the projected change (future minus historical) between NoahMP-WRF and the diel streamflow-based model, colored by the sensitivity of DOS 20 to climate change as projected by the empirical diel streamflow-based model (Fig. 5b). Symbols surrounded by black circles indicate sites that were excluded from the regression analysis in Fig. 3 (rainier site ID 24, 25, and 31). Symbols (circle, square, and diamond) represent the historical mean annual percentage of clear-sky snowmelt days, where sunny sites have > 90 % clear-sky snowmelt days, partly cloudy have between 70 % and 90 %, and cloudy have < 70 %; clear-sky snowmelt days are defined as those with more than 80 % of the potential clear-sky solar radiation. We make no inference about the cloudiness condition of snowmelt days under the RCP8.5 pseudo global warming climate scenario; however, red symbols (upper panels) follow the same symbology for easier interpretation.  Table A2. Root mean square error (RMSE) and coefficient of determination (R 2 ; in parentheses) associated with several stepwise multiple linear regressions (similar to the one in Eq. 1) using different early snowmelt timing metrics (e.g., Eq. 1 uses DOS 20 ) and correlation cutoffs (r) between hourly solar radiation and streamflow used to define snowmelt days. DOS xx represents the date when the xxth percentile of snowmelt days occurs. Bolded numbers are associated with the stepwise MLR in Eq. (1), which is also shown in Fig. 7a.
Early snowmelt timing metrics r > 0.5 r > 0.6 r > 0.7 r > 0.8  Table A3. Coefficient of determination (R 2 ) for the site-averaged stepwise multiple linear regression, analogous to that presented in Fig. 7b, for different modeling decisions (correlation cutoff between hourly solar radiation and streamflow, r, and early snowmelt days metrics). DOS xx represents the date when the xxth percentile of snowmelt days occurs. The bolded number is associated with the stepwise MLR in Eq.
Early snowmelt timing metrics r > 0.5 r > 0.6 r > 0.7 r > 0.  Table A4. Standardized beta coefficients for the stepwise MLR associated with the different correlation cutoffs (r) between hourly solar radiation and streamflow and different early snowmelt metrics. These stepwise MLR models follow the same structure as that of Eq. (1); however, in this case, predictors were standardized to estimate their relative importance. Note: AT -air temperature; Pp -precipitation; RH -relative humidity; SWR -incoming shortwave radiation. DOS xx represent the date when the xxth percentile of snowmelt days occurs. The asterisk * indicates rows that do not meet all the MLR assumptions. Bolded numbers are associated with the modeling decisions used in the result and discussion sections.
Early snowmelt timing metrics r > 0.5 r > 0.6 r > 0.7 r > 0.8 Code and data availability. Data from NoahMP-WRF simulations can be accessed through their public website at https://doi.org/10.5065/49SN-8E08 (Rasmussen et al., 2021). Hourly shortwave radiation can be accessed online at https://doi.org/10.5067/6J5LHHOHZHN4 (Xia et al., 2009). Hourly streamflow from the USGS database can be accessed online at https://waterdata.usgs.gov/nwis/sw (US Geological Survey, 2016). The code used to process and analyze the data presented in the study is available upon request to the corresponding author.
Author contributions. SAK and AH designed the study. SAK performed all the analyses, prepared the figures, and drafted the first version of the paper. JWK developed the diel cycle index, which served as the initial idea for the presented snowmelt detection method. GS collected and preprocessed USGS hourly streamflow data and NLDAS-2 solar radiation. LS preprocessed daily surface and subsurface runoff from the WRF CONUS-I simulations. All the authors reviewed and contributed to the final version of the paper.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.