METHODS

3.1 Hydrological modeling

In this study, a combination of hydrological model outputs and observations were used to investigate linkages between hydrological characteristics at local scales and climatic teleconnection patterns over a large domain. A combination of modeling and observational data were used to comprehensively investigate the linkages because even intensively instrumented watersheds such as Reynolds Mountain do not record all of the necessary water and energy balance fluxes. For instance, measured precipitation, modeled rainfall, modeled snowfall, and measured streamflow or annual runoff were used. To gain confidence in the model used, its performance in simulating snow was assessed.
The cold regions hydrological modeling (CRHM) platform (Pomeroy et al., 2007) has been widely used and tested in various agricultural and mountainous regions in Canada and other cold areas such as the Canadian Rockies, Qinghai-Tibetan Plateau, Patagonia, the Pyrenees and the Alps (Dornes, Pomeroy, Pietroniro, Carey, & Quinton, 2008; Ellis & Pomeroy, 2007; Ellis, Pomeroy, Brown, & MacDonald, 2010; Fang et al., 2013; Lòpez-Moreno, Pomeroy, Revuelto, & Vicente-Serrano, 2013; Mahmood, Pomeroy, Wheater, & Baulch, 2017; Pomeroy, Fang, & Rasouli, 2015; Rasouli, Pomeroy, Janowicz, Carey, & Williams, 2014; Rasouli et al., 2015; Rasouli, Pomeroy, & Whitfield, 2019a, 2019b; Rasouli, 2017). A strength of CRHM is its ability to capture snow redistribution and sublimation, and thus, it provides an understanding of the spatial and temporal heterogeneity of the hydrological processes in snow dominated watersheds. In this study, components of the hydrological water balance were simulated using the CRHM platform. The model was forced with observation data and run on homogeneous hydrological response units (HRUs) to characterize the surface and near-surface cold regions hydrological processes. The model performance in capturing the spatial variability of snow depth was assessed in Reynolds Mountain.
A set of physically based modules describing the major processes were compiled (using CRHM platform) into a watershed model informed by results from previous scientific findings and modeling experiments in research basins (Carey, Quinton, & Goeller, 2007; Dornes et al., 2008; Flerchinger, Reba, & Marks, 2012; Link, Flerchinger, Unsworth, & Marks, 2004; McCartney, Carey, & Pomeroy, 2006; MacDonald, Byrne, Kienzle, & Larson, 2011; Pomeroy, Hedstrom, & Parviainen, 1999; Pomeroy, Toth, Granger, Hedstrom, & Essery, 2003; Pomeroy et al., 2006; Quinton & Carey, 2008; Reba, Pomeroy, Marks, & Link, 2012; Reba, Marks, Link, Pomeroy, & Winstral, 2014; Williams et al., 2015; Winstral, Marks, & Gurney, 2013). The phase of precipitation, whether it falls as snow or rain, was partitioned based on the psychrometric energy balance of the falling hydrometeors and turbulent transfer equations of blowing snow sublimation (Harder & Pomeroy, 2013). Two energy balance snowmelt and evapotranspiration modules were used in the CRHM model. The snowfall and rainfall intercepted on the canopy were estimated and used to calculate the sublimation and evaporation losses. The turbulent transfer to snow and wind-driven snow transport and redistribution were modeled using CRHM to better close the mass balance in the basin. No calibration was applied for snow simulations. Runoff generated during rain on snow (ROS) events was estimated from an energy balance model and calculating advective heat flux from precipitation when snow is on the ground and the equivalent melt from advected heat to snowpack (Marks, Link, Winstral, & Garen, 2001). In the developed model, vegetation height, density, and stalk diameter control the aerodynamic roughness and play a key role in blowing snow transport. Thus, snow is transported from shorter vegetation HRUs such as sagebrush to taller vegetation HRUs such as riparian forest and aspen. The modules and references to the methods used are listed in Table 1. For more details on each module and each hydrological process refer to Rasouli et al. (2014) and Rasouli (2017).

3.2 The spatial arrangement of the hydrological model

The water and energy balance fluxes were modeled for HRUs with hourly time steps. HRUs that act as control volumes of the calculation were spatially segregated based on surface physiographic information relevant for hydrological model parameterization, including vegetation cover, topography, soil depth, soil layers, the variability of basin attributes, and the level of model complexity. A spatially distributed modeling structure was developed with 22 HRUs (Figure 3). HRU configuration was adapted from previous studies (e.g., Newman, Clark, Winstral, Marks, Seyfried, 2014; Rasouli et al., 2015) and extended based on aspect and slope. Drift HRUs in aspen and sage vegetation are in topographic depressions downwind of slope breaks, where deep snow accumulates, which sustain snowmelt water for the trees in spring and early summer.

3.3 Physical parameterization of the hydrological model

Parameter estimation in Reynolds Mountain was based on previous studies in this basin, other headwater basins in RCEW, and similar snow dominated basins. Model parameters adapted from measurements in the research basins to represent the vegetation characteristics, snowmelt, and blowing snow redistribution. A uniform blowing snow fetch distance was used for all HRUs due to the short upwind distance. Blowing snow was inhibited for the sheltered HRUs. Snow surface roughness length was estimated by Reba et al. (2012). Vegetation heights/density and leaf area index were determined by Link et al. (2004), Seyfried et al. (2009), and Flerchinger et al. (2012). The soil characteristics, including soil water storage capacity and the soil surface saturation, were adapted from Seyfried et al. (2009) and Link et al. (2004), respectively. Initial soil temperature was measured by soil thermocouples prior to the major snowmelt. The thermal conductivity value was taken from Oke (1978).

3.4 Classification of hydroclimatic phases

The water balance cycle of a basin cannot be closed unless the effects of hydrological and meteorological conditions in preceding years are taken into account, in addition to the seasonal precipitation and runoff. There are multiple methods in the literature to break down time series into distinct components that represent different temporal frequencies of the variations (e.g., wavelet transform (Huang et al., 1998) and empirical mode decomposition (Daubechies 1992)). Alexandrov et al. (2012) found the singular spectrum analysis (SSA) method to have the highest performance. We used the SSA method to decompose daily precipitation time series. SSA is typically used for spectral decomposition of a time series, estimating the spectrum of eigenvalues in a singular value decomposition of a covariance matrix. The components of SSA are trend, periodic components, and noise, each having a meaningful interpretation. The basic structure of SSA includes: (1) decomposing the one-dimensional time series into multidimensional series by adding lagged variables and creating a so-called trajectory matrix; and (2) conducting a principal component analysis on the trajectory matrix (Hsieh, 2009).
Precipitation data were used for SSA and to identify the hydroclimatic phases, because it dominates the variability in the water budget of a basin. Precipitation also has seasonality and serial dependencies, which can be decomposed by SSA. The seasonality can be multiple seasons, years, or decades, depending on the nature of the phenomena. The first mode or eigenvector of SSA on precipitation usually shows an intermediate frequency variation, which is ideal for representing hydroclimatic phases with multiple year durations studied herein. The long term average of the first eigenvector of SSA on precipitation is used as the threshold to indicate whether a given hydroclimatic phase shows a wet or dry condition. If the moving average of the first eigenvector of SSA on precipitation is above the long term average for a period of three to eight years, it is defined as a positive hydroclimatic phase, and if it is below the long term average, it is defined as a negative hydroclimatic phase.
Climatic teleconnection patterns represent regional atmospheric and oceanic circulations over multiple years and decades. The teleconnection patterns used in this study are Antarctic Oscillation (AAO), Sea Surface Temperature (SST), Arctic Oscillation (AO), North Atlantic Oscillation (NAO), and Pacific-North American (PNA) pattern. These teleconnection patterns were chosen because large-scale atmospheric circulation changes occur via these teleconnections (Overland & Wang, 2010) and there is evidence of their effects on snow regimes in western North America (e.g., Bao, Kelly, & Wu, 2011). For instance, NAO is one of the most dominant teleconnection patterns in the atmosphere and is similar to AO in describing the difference in pressure gradient between Azores high and Icelandic low in a different time period (Hurrell, Kushnir, Ottersen, & Visbeck, 2003). During the positive phase of NAO, the jet stream blows strongly and consistently from west to east, which locks up the cold and dry Arctic air in the polar region, leading to warm and wet winters in North America. But during the negative phase of NAO, the reduced pressure gradient weakens the westerly wind allowing the cold Arctic air to move to the mid-latitude region (Francis & Vavrus, 2012; Cohen et al., 2014), leading to cold winters in North America. SST in the Niño 3.4 area of the Equatorial Pacific Ocean also can affect ROS events. McCabe, Clark, and Hay (2007) showed that there are more frequent ROS events in the northwestern United States during La Niña conditions than El Niño. Low flow conditions in western North America are also associated with El Niño events and a positive phase of the PNA pattern (Bonsal & Shabbar, 2008). AAO is a large-scale mass transport in the atmosphere between the mid‐latitudes and high latitudes (Gong & Wang, 1999).
The climatic teleconnection patterns and anomalies of observed annual precipitation, mean winter air temperatures, the modeled rainfall to precipitation and runoff to precipitation ratios were compared in each hydroclimatic phase against their long term averages. A correlation analysis was also conducted among the meteorological and hydrological variables to determine the importance of the atmospheric teleconnection patterns in small-scale watersheds. We expect that the detected hydroclimatic phases to be aligned with negative or positive phases of the teleconnection patterns, and snow, ROS runoff, and streamflow to behave differently in each phase. As runoff generated during ROS events contributes significantly to the annual runoff, its spatial and temporal variability was examined in detail under different climate phases.