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.