Interactive comment on “ Observations of seasonal and diurnal glacier velocities at Mount Rainier , Washington using terrestrial radar interferometry ”

Abstract. We present surface velocity maps derived from repeat terrestrial radar interferometry (TRI) measurements and use these time series to examine seasonal and diurnal dynamics of alpine glaciers at Mount Rainier, Washington. We show that the Nisqually and Emmons glaciers have small slope-parallel velocities near the summit ( 10 km range can be used to investigate spatial and temporal variability of alpine glacier dynamics over large areas, including hazardous and inaccessible areas.


Introduction
Direct observations of alpine glacier velocity can help improve our understanding of ice dynamics.Alpine glacier surface velocities are typically dominated by basal sliding, which is tightly coupled to subglacial hydrology (Anderson et al., 2014;Bartholomaus et al., 2007).However, the spatial extent and spatial/temporal resolution of direct velocity measurements are often limited to short campaigns with sparse point measurements in accessible regions (e.g., Hodge, 1974;Driedger and Kennard, 1986).Remote sensing can help overcome many of these limitations.Radar interferometry, a form of active remote sensing, detects millimeter-to-centimeterscale displacements between successive images of the same scene and can see through clouds and fog.In the past few decades, satellite interferometric synthetic aperture radar, or InSAR (e.g., Massonnet and Feigl, 1998;Burgmann et al., 2000), has emerged as an invaluable tool for quantifying glacier dynamics (e.g., Joughin et al., 2010).However, limited data availability and revisit times limit the application of InSAR for the study of many short-term processes.
Terrestrial radar interferometry (TRI), also referred to as ground-based radar interferometry, has recently emerged as a powerful technique for observing glacier displacement that is not prone to the same limitations (Caduff et al., 2014).Sets of radar data acquired at intervals as short as ∼ 1 min from up to several kilometers away allow for observations of velocity changes over short timescales and large spatial extents.Stacking these large numbers of interferogram pairs over longer timescales can significantly reduce noise.Here, we employ this relatively new technique to provide spatially and temporally continuous surface velocity observations for several glaciers at Mount Rainier volcano in Washington State (Fig. 1).Though Rainier's glaciers are among the best-studied alpine glaciers in the USA (Heliker et al., 1984;Nylen, 2004), there are many open questions about diurnal and seasonal dynamics that TRI can help address.Specifically, many aspects of subglacial hydrology and its effects on basal sliding are poorly constrained, especially for inaccessible locations like the Nisqually icefall and ice cliff.Our observations provide new insight into these processes through analysis of the relative magnitude and spatial patterns of surface velocity over diurnal and seasonal timescales.To our knowledge, no other studies have investigated seasonal changes to glacier dynamics using TRI.
Mount Rainier offers an excellent setting for TRI, with several accessible viewpoints offering a near-continuous view with ideal line-of-sight vectors for multiple glaciers, and well-distributed static bedrock exposures for calibration.The ability to image the velocity field of entire glaciers from one viewpoint with minimal shadowing sets this study area apart.Most previous studies only image part of the glaciers under investigation, usually due to less favorable viewing geometries (e.g., Noferini et al., 2009;Voytenko et al., 2015;Riesen et al., 2011).However, the steep topography and local climatic factors at Mount Rainier result in strong atmospheric variability and turbulence -a major source of noise for radar interferometry techniques (Goldstein, 1995).Atmospheric noise is a particular issue for the long ranges (> 10 km) associated with accessible viewpoints at Mount Rainier.To overcome this limitation, we successfully combine, expand on, and evaluate noise reduction techniques such as stacking interferograms (e.g., Voytenko et al., 2015) and deriving atmospheric noise corrections over static control surfaces (bedrock exposures) (e.g., Noferini et al., 2009).We demonstrate that these techniques offer significant uncertainty reductions using a novel bootstrapping approach.
In the following sections, we provide background on Mount Rainier's glaciers and detail our sampling methodology and data processing techniques.We then present TRI results documenting seasonal and diurnal velocity variations for the Nisqually, Wilson, Emmons, and upper Winthrop glaciers and quantify measurement uncertainty.Next we examine the partitioning of observed surface velocities between deformation and basal sliding at different times of year using a simple 2-D flow model and compare our observations to other recent and historical velocity measurements.These comparisons provide ground truth for TRI measurements and new insight into the evolution of the Nisqually Glacier since the late 1960s.

Study area
With a summit elevation of 4392 m, Mount Rainier (Fig. 1) is the largest stratovolcano in the Cascades and is considered the most dangerous volcano in the USA (Swanson et al., 1992).It also holds the largest concentration of glacial ice in mainland USA (Driedger and Kennard, 1986) -87 km 2 was covered with perennial snow and ice in 2008 (Sisson et al., 2011).The steep upper sections of the major glaciers are relatively thin, with typical thicknesses of ∼ 30-80 m (Driedger and Kennard, 1986).Thickness increases at lower elevations, with a maximum of ∼ 200 m for the Carbon Glacier, although these estimates likely provide an upper bound, as these glaciers have experienced significant thinning in recent decades, losing 14 % of their volume between 1970 and 2008 (Sisson et al., 2011).Mass balance stake measurements from 2003 to 2010 show that the average winter balance for Nisqually was 2.4 m water equivalent (m w.e.), average summer balance was −3.5 m w.e., and cumulative net balance was −8.6 m w.e. from 2003 to 2011 (Riedel, 2010;Riedel and Larrabee, 2015).
The glaciers of Mount Rainier have been of interest to geoscientists for over 150 years and have a long record of scientific observation (Heliker, 1984).In this study, we focus on large, accessible, well-documented glaciers in the park: the Nisqually Glacier on the southern flank and Emmons Glacier on the northeastern flank.Additional glaciers in the field of view are also captured, including the Wilson Glacier, which flows into the Nisqually Glacier, the upper Winthrop Glacier, Fryingpan Glacier, upper Kautz Glacier, and Inter Glacier.All glaciers are labeled in Fig. 1.
The Nisqually Glacier is visible from several viewpoints near the Paradise Visitor Center, which is accessible year round.The terminus location has been measured annually since 1918, and three transverse surface elevation profiles have been measured nearly every year since 1931 (Heliker, 1984).Veatch (1969) documented a 24-year history of Nisqually's advances and retreats and other dynamic changes through a meticulous photographic survey from 1941 to 1965.Hodge (1974) conducted a detailed 2-year field study of the seasonal velocity cycle for the lower Nisqually.He found that velocities varied seasonally by about 50 %, with maximum velocities in the spring (June) and minimum in the fall (November).This finding, and the lack of correlation between runoff and sliding speeds, advanced the idea that efficient conduits close as meltwater input decreases in the fall, leading to distributed subglacial storage through the fall, winter, and spring.Increased surface melting in spring and summer leads to increased subglacial discharge and the opening of a more efficient network of conduits capable of releasing some of this stored water (Hodge, 1974).More recently, Walkup et al. (2013) tracked the movements of supraglacial rocks with high precision from 2011 to 2012, yielding velocity vectors for a wide network of points over the lower parts of Nisqually Glacier.
The Emmons Glacier, visible from the Sunrise Visitors Center, has received less attention than Nisqually, despite the fact that it is the largest glacier by area on the mountain (Driedger and Kennard, 1986), mainly because it is not as easily accessible as Nisqually.A large rock fall (∼ 1.1×10 7 m 3 ) from Little Tahoma in December 1963 covered much of the lower Emmons Glacier with a thick debris layer (Crandell and Fahnestock, 1965).The insulating debris cover likely contributed to the advance and thickening of the Emmons Glacier from 1970 to 2008, while all other glaciers on Mount Rainier experienced significant thinning (Sisson et al., 2011).The average 2003-2010 winter balance for Emmons was 2.3 m w.e., the average summer balance was −3.2 m.w.e, and the cumulative net balance was 7.7 m w.e. from 2003 to 2011 (Riedel, 2010;Riedel and Larrabee, 2015).
The National Park Service's long-term monitoring protocols include both the Nisqually and Emmons glaciers and involve regular photographs, annual mass balance measurements, meltwater discharge rates, plus area and volume change estimates every decade (Riedel, 2010;Riedel and Larrabee, 2015).

Instrument description
For this study, we used a GAMMA portable radar interferometer (GPRI) (Werner et al., 2008(Werner et al., , 2012) )  capture millimeter-scale surface displacements.The instrument includes three 2 m antennas mounted on a vertical truss, with one transmit antenna 35 cm above the upper of two receiving antennas, spaced 25 cm apart (Fig. 2).The transmit antenna produces a 35 • vertical beam with 0.4 • width that azimuthally sweeps across the scene to build a 2-D radar image as the truss rotates.The radar operates at a center frequency of 17.2 GHz, with selectable chirp length of 2-8 ms and bandwidth of 25-200 MHz.The radar wavelength is 17.6 mm with range resolution of ∼ 0.75 cm and one-way interferometric change sensitivity of 8.7 mm cycle −1 of phase providing < 1 mm line-of-sight precision.Line-of-sight interferograms are generated by comparing phase differences in successive acquisitions from the same viewpoint.The interval between acquisitions can be as short as ∼ 1 min, allowing for high coherence even in rapidly changing scenes.

Survey Description
We performed four data collection campaigns in 2012 (Table 1).The first campaign occurred on 6-7 July 2012.This timing corresponds to just after the expected peak seasonal glacier velocities at Mount Rainier (Hodge, 1974).Following the success of this study, three subsequent deployments were performed during the late fall and early winter, which should capture near-minimum seasonal velocity (Hodge, 1974).These campaigns were timed to occur before, immediately after, and a few weeks after the first heavy snowfall of the season (2 and 27 November and 10 December 2012, respectively).Distances from the GPRI to the summit were 6.7, 7.6, and 10.8 km from GLPEEK, ROI, and SUNRIZ, respectively.Radar images were continuously collected with a 3 min interval for all surveys.Total acquisition time at each site was dictated by logistics (weather conditions, personnel), with ∼ 24 h acquisitions at SUNRIZ and ROI to capture diurnal variability.
The instrument was deployed on packed snow during the 6 July 2012 GLPEEK and 27 November and 10 December 2012 ROI acquisitions.Over the course of the GLPEEK survey, we noted limited snow compaction and melt beneath the GPRI tripod with total displacement of ∼ 2-4 cm over ∼ 6 h.However, this instrument motion proved to be negligible for the interferogram interval used (6 min).We did not note significant snow compaction under the tripod during the fall/winter surveys.
Weather conditions during the July 2012 surveys were clear with light/variable wind.The 2 November 2012 survey involved high-altitude clouds, passing showers, and brief interruptions in data collection.Weather conditions were clear with sun for the 27 November 2012 campaign and fog with limited visibility on 10 December 2012.

Data processing
All radar data were processed with the GAMMA SAR and Interferometry software suite.Interferograms were generated from single-look complex SLC products with a time separation of 6 min, though sometimes longer if acquisition was interrupted.For example images see Fig. A4.Interferograms were multi-looked by 15 samples in the range direction to reduce noise.A correlation threshold filter of 0.7 and an adaptive bandpass filter with default GAMMA parameters were applied to the interferograms to improve phase unwrapping.Phase unwrapping was initiated in areas with high correlation scores and negligible deformation, such as exposed bedrock or stagnant ice.

Atmospheric noise corrections
Slight changes in the dielectric properties of the atmosphere between the GPRI and target surfaces can lead to uncertainty in the interferometric displacement measurements (Zebker et al., 1997;Werner et al., 2008).Changes in atmospheric humidity, temperature, and pressure can all affect radar propagation velocity (Goldstein, 1995).These variations are manifested as phase offsets in the received radar signal, which must be isolated from phase offsets related to true surface displacements.
This atmospheric noise proved to be significant for the long range (i.e., ∼ 22 km two-way horizontal path at SUN-RIZ), mountainous terrain (i.e., ∼ 2.4 km vertical path from SUNRIZ to summit), and turbulent atmosphere involved with this study, with the magnitude of this noise often exceeding that of surface displacement signals.The scale of the atmospheric noise features we observed in the data was typically much wider than the width of the glaciers, so in order to minimize this atmospheric noise in the individual interferograms, we interpolated apparent displacement values over static control surfaces (e.g., exposed bedrock).To do this, we fit a surface using Delaunay triangulation to a subset (5 %) of pixels over exposed bedrock.The subset of pixels was resampled randomly for each unwrapped interferogram and the interpolated result was smoothed to reduce artifacts and then subtracted from the interferogram.The corrections were applied to all individual interferograms, and the resulting products were stacked to further reduce noise.To stack, we took all the images for a given time period and computed the mean and median at each pixel.This has the effect of augmenting signal and canceling out noise.The median is less affected by outliers and is our preferred result.The median line-of-sight (LOS) velocities from this stack provide a single measurement with a high signal-to-noise ratio for the entire sampling period.
In addition to computing the median LOS velocities for the entirety of each sampling period, we also computed a running mean of the LOS velocities to characterize any shortterm velocity variations in the extended occupation data sets: 7-8 July SUNRIZ (24 h) and 1-2 November ROI (21 h).The running mean was computed every 0.3 h with a 2-h centered (acausal) window, with standard error used to estimate uncertainty.
Interferograms with significant phase unwrapping errors, low correlation, or anomalous noise were excluded from stacking.We only excluded a few images for each site with the exception of SUNRIZ, which produced many images with anomalous noise and unwrapping errors, possibly due to instrument noise and/or the extended range through significant atmospheric disturbance.For this reason, more than half of the data from SUNRIZ were excluded from the analysis (Table 2).For GLPEEK and ROI, interferograms with occasional localized unwrapping errors were preserved during stacking, as they have little influence on the final stack median.However, localized areas with persistent unwrapping errors in the SUNRIZ data were masked using a threshold standard deviation filter of 0.6 m day −1 .
We estimated median LOS velocity uncertainties using a bootstrapping approach (Efron, 1979).This involved resampling the set of images used in the stack with replacement 1000 times for each campaign.Then, for each pixel, the 25th and the 975th ordered values were set as the lower and upper bounds of the 95 % confidence interval.

Conversion from radar coordinates to map coordinates
We developed a sensor model and tools to terrain-correct the stacked GPRI data (in original azimuth, range coordinates) using an existing 2 m pixel −1 airborne LiDAR digital elevation model (DEM) acquired in September 2007/2008(Robinson et al., 2010).While some elevation change has undoubtedly occurred for glacier surfaces between September 2008 and July 2012, the magnitude of these changes (< 20 m) is negligible for orthorectification purposes given the GPRI acquisition geometry.A single control point identified over exposed bedrock in the LiDAR DEM and the multi-look image radar data was used to constrain absolute azimuth orientation information for each campaign.A ∼ 10 m pixel −1 (mean of azimuth and range sample size) grid in UTM 10N (EPSG: 32610) was created for each campaign, with extent computed from the GPRI GPS coordinates, min/max range values, and min/max absolute azimuth values.Each 3-D pixel in this grid was then populated by extracting the radar sample with corresponding range and azimuth.

Correction to slope-parallel velocities
While the line-of-sight vectors for these surveys are roughly aligned with surface displacement vectors (median incidence angles for glacier surfaces are ∼ 22 • for GLPEEK, ∼ 25 • for SUNRIZ, and ∼ 26 • for ROI), glaciological analyses typically require horizontal and vertical velocity components relative to the glacier surface.As each GPRI survey offers only a single look direction, this is not possible.However, we can assume that displacement is dominated by surfaceparallel flow and use the 2007/2008 LiDAR DEM to extract surface slopes needed to estimate 3-D displacement vectors (e.g., Joughin et al., 1998).This approach is intended for relatively smooth, continuous surface slopes over length scales > 2-3× ice thickness.It is therefore possible that the slope-parallel correction can overestimate velocity for steep, high-relief surfaces with significant high-frequency topographic variability (e.g., www.the-cryosphere.net/9/2219/2015/The Cryosphere, 9, 2219-2235, 2015 icefalls).The slope-parallel assumption also begins to break down where the vertical flow velocity component becomes significant.This is expected in the upper accumulation and lower ablation zones, where the submergence and emergence velocities become more significant, respectively, but is less important near the equilibrium-line altitude (ELA) or locations where sliding dominates surface motion.The latter is expected for much of the Nisqually Glacier at least (Hodge, 1974).We implement a slope-parallel correction by first downsampling the 2007/2008 LiDAR DEM to 20 m pixel −1 and smoothing with a 15 × 15 pixel (∼ 300 m), 5 σ Gaussian filter.The slope-parallel velocity (V sp ) is defined as where Ŝ • L is the dot product between the unit vector pointing directly downslope from each grid cell ( Ŝ) and the unit vector pointing from each grid cell to the sensor ( L). Regions where the angle between these two vectors exceeded 80 • were masked to avoid dividing by numbers close to 0 which could amplify noise.

Two-dimensional glacier deformation modeling
Surface flow velocity can be partitioned into internal deformation and basal sliding components.We present a simple, 2-D plane-strain ice deformation model for a preliminary assessment of the importance of basal sliding for the glaciers in our study area.The deformation model uses the shallow ice approximation (SIA) -an approximate solution of the Stokes equations (Greve and Blatter, 2009;Cuffey and Paterson, 2010).The expected surface velocity u s due to internal deformation from the SIA model is where ρ i represents ice density, g represents gravitational acceleration, α represents local surface slope, H represents local ice thickness, A represents an ice softness parameter, and n represents a flow rate exponent.The coordinate system is vertically aligned.The SIA is not well suited for narrow mountain glaciers, so we modify it to simulate the effect of non-local conditions, such as lateral sidewall drag and longitudinal stretching.The ice thickness H and surface slope α are smoothed using a weighting function based on Kamb and Echelmeyer (1986).Kamb and Echelmeyer (1986) calculated a longitudinal coupling length l using a 1-D force balance approach, for each point in their domain.They calculated l to be in the range of one to three ice thicknesses for valley glaciers.We simplified this by using a single value for l over the domain of model.The longitudinal couple length l is used in a weighting function to smooth α and H .The weighting function has where x and y represent the horizontal coordinates of the weight position, and x and y represent the horizontal coordinates of the reference position.Weights are calculated at each point in the model domain, over a square reference window (side length of A w ).H and α are smoothed at the reference position by normalizing weights over the reference window.We choose a coupling length l of ∼ 1.5 ice thickness and an averaging window size of ∼ 3 ice thicknesses, consistent with the usage in Kamb and Echelmeyer (1986).We use a spatially uniform and temporally constant ice softness parameter suitable for ice at the pressure melting point of 2.4 × 10 −24 Pa −3 s −1 (Cuffey and Paterson, 2010, p. 75).
Ice softness can be affected by several factors (e.g., englacial water content and impurities), so we also consider an ice softness parameter up to twice this best estimate in accounting for model uncertainties, as described below.Our best estimates of model input parameters are summarized in Table 3. Surface slope (Fig. A1b) was estimated from the 2007/2008 LiDAR DEM (Robinson et al., 2010).Surface velocities u s are the TRI-derived median slope-parallel velocities.Ice thicknesses H (Fig. A1a) were estimated by differencing the 2007/2008 LiDAR DEM surface elevations and the digitized and interpolated bed topography from Driedger and Kennard (1986).The Driedger and Kennard (1986) bed topography contours were derived from icepenetrating radar point measurements and surface contours from aerial photographs.The published basal contours for Nisqually/Wilson, Emmons, and Winthrop glaciers were digitized and interpolated to produce a gridded bed surface using the ArcGIS Topo to Raster utility.The gridded bed elevations have a root mean squared error of 11 m when compared with the 57 original radar point measurements.A point-toplane iterative closest point algorithm (implemented in the NASA Ames Stereo Pipeline pc_align utility; Shean et al., 2015) was used to coregister the 1986 bed topography to the 2007/2008 LiDAR topography over exposed bedrock on val- ley walls.Mean error over these surfaces was 7.6 m following coregistration, although some of this error can be attributed to actual surface evolution near glacier margins (e.g., hillslope processes) from 1986 to 2008.In addition to these interpolation and coregistration errors, there were likely small changes in ice thickness during the 4-5 years between the 2007/2008 DEM data collection and the 2012 TRI observations, as mass balance measurements suggest that both the Nisqually and Emmons glaciers experienced net mass loss during this time period (Riedel and Larrabee, 2015).Propagation of these uncertainties results in estimated ice thickness uncertainties of ∼ 5-25 %.In order to account for this large uncertainty, we ran the model with ±25 % ice thickness as well as 2× ice softness in order to estimate the possible range of expected deformation velocities.
More sophisticated ice flow models (e.g., Gagliardini et al., 2013;Le Meur et al., 2004;Zwinger et al., 2007) could potentially offer a more realistic picture of the spatial and temporal variability of glacier sliding.However, given the poorly constrained model inputs and observational emphasis for this study, we proceed with the SIA model to obtain approximate estimates for the deformation and sliding components of observed velocities.

Results
The median stacks of surface-parallel velocity for all viewpoints and their respective uncertainty estimates are shown in Figs.3-6.Overall, our results show that repeat TRI measurements can be used to document spatial and temporal variability of alpine glacier dynamics over large areas from > 10 km away.The atmospheric noise removal approach was success-ful in extracting a glacier displacement signal for all campaigns, with excellent results for Nisqually Glacier due to the shorter range from ROI and GLPEEK viewpoints and limited glacier width between control surfaces.Stacking alone was very effective; the velocities of the mean and median stacks with and without the atmospheric noise correction were very similar.The main benefit of the extra step of using stable rock points to subtract an estimate of the atmospheric noise was to significantly reduce the uncertainties and to reduce the noise where velocities are slow.The uncertainties before and after atmospheric correction are compared on Table 2.The median width of the 95 % confidence interval for each corrected, stacked pixel is plotted in Figs.3b  and 5.Note near-zero values over exposed bedrock surfaces used to derive atmospheric noise correction.We were able to reduce uncertainties (half the median confidence interval width) to about ±0.02 to ±0.08 m day −1 over glacier surfaces for some campaigns, with uncertainty dependent on the total number of stacked images, weather conditions, and target range (Table 2).For example, the 6 July 2012 ROI survey had a final confidence interval width of 0.11 m day −1 (∼ ±0.06 m day −1 ) while the 10 December 2012 ROI survey had a final confidence interval width of 0.15 m day −1 (∼ ±0.08 m day −1 ) despite a 50 % increase in stack count.This is likely due to increased local atmospheric variability, as low-altitude clouds obscured the surface during 10 December 2012 survey.The 2 November 2012 ROI survey had the highest stack count (359) with the lowest uncertainty values of ±0.02 m day −1 (Table 2).

July 2012 surface velocities
The 6-7 July 2012 observations show slope-parallel velocities that range from ∼ 0.0 to 1.5 m day −1 for both the Nisqually and Emmons glaciers (Figs.3a, 4, 6).Both display high velocities over their upper and central regions that taper into essentially stagnant (< 0.05 m day −1 ) debris-covered regions near the terminus.In general, slope-parallel velocities near the summit are small (< 0.2 m day −1 ).
On the Nisqually Glacier, a series of local velocity maxima (> 1.0 m day −1 ) are associated with increased surface slopes between local surface highs.Local velocity maxima are also observed for the fast-flowing Nisqually icefall (western branch of upper Nisqually, see Fig. 3) and above the Nisqually ice cliff (eastern branch).A relatively smooth velocity gradient from slow-to fast-moving ice is present upstream of the icefall, while the velocities above the ice cliff display a steep velocity gradient (Fig. 3).The main (south) branch of the Emmons Glacier displays generally increasing velocity from the summit to lower elevations.A large high-velocity region (> 0.7-1.1 m day −1 ) is present over central Emmons, downstream of the confluence of upper branches.These elevated velocities decrease at lower elevations, where ice thickness increases and surface slopes decrease (Fig. A5).A central "core" of exposed ice displays slightly elevated velocities relative to surrounding debris-covered ice within ∼ 1-1.5 km of the terminus.
Velocities exceed 1 m day −1 over the "central" branch of the upper Emmons Glacier, where flow is restricted between two parallel bedrock ridges, with local maxima similar to Nisqually.Velocities at higher elevations within the "central" branch appear slower (< 0.1-0.5 m day −1 ), separated from the fast downstream velocities by a small area that was excluded due to phase unwrapping errors.Photographs show that this area appears heavily fractured with many large blocks indicative of rapid, discontinuous flow (Fig. A3).
Smaller, relatively thin glaciers, such as the Fryingpan, upper Kautz, and Inter glaciers (labeled in Fig. 1), also display nonzero surface velocities of < 0.1-0.2m day −1 , but with limited spatial variability.

Seasonal variability
The repeat observations from the ROI viewpoint provide time series that capture seasonal velocity variability for the Nisqually, Wilson, and upper Kautz glaciers.We observe significant velocity changes during the summer to winter transition and more subtle changes within the winter period.These changes are shown in map view in Fig. 4 and in profile view with corresponding slope and ice thickness in Fig. 6.
These data show a velocity decrease of 0.2-0.7 m day −1 (−25 to −50 %) from July to November 2012 for most of the Nisqually Glacier.This includes central and lower Nisqually and the ice above the ice cliff.The greatest velocity decreases are observed near the crest and lee of surface rises (down-  (Robinson et al., 2010).Refer to Fig. 5 and Table 2 for uncertainty estimates.
stream of data gaps from radar shadows; Fig. 4), where some of the highest velocities were observed in July.In contrast, the area immediately downstream of the ice cliff and the area surrounding the icefall both display an apparent velocity increase for the same time period (Figs. 4,6).While the increase is less than the 95 % confidence interval for most areas, we can confidently state that the icefall and area below the ice cliff do not display the significant decrease in velocity observed elsewhere.
The majority of the Wilson Glacier displays a similar ∼ 0.3-0.7 m day −1 (−40 to −60 %) velocity decrease from July to November.Interestingly, the steep transition where the Wilson merges with the Nisqually displays an apparent velocity increase of ∼ 0.1 m day −1 during this time period (Fig. 4).These data also reveal subtle velocity increases in the debris-covered ice near the Nisqually terminus and the upper Kautz Glacier (Fig. 4), though these increases are statistically insignificant.
The repeat winter observations of Nisqually show relatively constant velocities with some notable variability.Anal- (+130 %) increase over upper Wilson.In the latter case, the 10 December velocities are actually higher than those observed in July.The slowdown over lower Nisqually appears robust, but other trends have amplitudes that are mostly below the 95 % confidence interval for the 27 November and 10 December observational campaigns (Fig. 4).

Diurnal variability
We collected ∼ 21 and ∼ 24 h time series for the Emmons and Nisqually/Wilson glaciers (Table 1) in July and November, respectively, and looked at changes throughout the day.
Although uncertainties are large, we present the time series in Fig. 7.
In general, velocities for these regions remain relatively constant during their respective sampling periods.The Emmons time series shows an apparent decrease in velocity over the central, fast-flowing regions (B, C, D in Fig. 7a) from ∼ 18:00 to 21:00 LT and an apparent increase between ∼ 07:00 and 09:00 LT (Fig. 7a).The Nisqually time series shows an apparent decrease from ∼ 06:00 to 11:00 LT for the icefall and ice cliff and an apparent decrease for several areas of the glaciers followed by an increase (Fig. 7b).How-  Walkup et al. (2013) and TRI velocities at Walkup et al. (2013) sample locations (Fig. 8).Walkup et al. (2013) between 19 July and 11 October 2012 (black) compared to TRI slope-parallel velocities derived from this study at the same locations for two time periods that bracket the time period measured by Walkup et al. (2013).See Table 4 for comparison statistics and Box C in Fig. 1 for context.ever, uncertainties are large and none of these are statistically significant.

Comparison with independent velocity measurements
We now compare our TRI results with independent velocity measurements for an overlapping time period.Walkup et al. (2013) performed repeat total station surveys to document the location of sparse supraglacial cobbles and boulders on the lower Nisqually Glacier from 2011 to 2012.While measurement errors (e.g., cobble rolling/sliding) for these observations are difficult to document, the large sample size and relatively long measurement intervals allow for accurate surface velocity estimates.
Figure 8 shows average velocity vectors measured by Walkup et al. (2013) for the period between 19 July and 11 October 2012, with corresponding surface-parallel velocity vectors from the 7 July and 2 November TRI surveys.This comparison is summarized on Table 4.In general, the velocity magnitudes are similar, with the overall mean of the Walkup et al. (2013) measurements slightly higher on average but often falling between the 7 July and 2 November GPRI magnitudes, as would be expected of a mean velocity spanning approximately the same period.The velocity directions are also relatively consistent, with a median angular difference of 12 • .The greatest deviations are observed near the ice margins and over small-scale local topography (e.g., icecored moraine near western margin), where surface-parallel flow assumptions break down.In general, the two techniques provide similar results and offer complementary data validation.However, since the Walkup et al. (2013) measurements were limited to accessible areas, they cannot be used to validate TRI observations for heavily crevassed areas, icefalls, and other hazardous dynamic areas generally higher on the mountain.

Two-dimensional flow modeling
Figure 9 shows modeled deformation, sliding velocity residual (observations-deformation model), and sliding percent (sliding velocity residual as percentage of total velocity) with best estimate model parameters for Nisqually Glacier in July and November.Figure 10 shows corresponding output for Emmons.The SIA deformation models suggest that most areas of both glaciers are moving almost entirely by sliding.The modeled glacier deformation alone is unable to account for the observed surface velocity during any of the observation periods.Only a median of 1 % of the velocity field over the Nisqually Glacier area can be explained by internal deformation in July and only 2 % in November.If we consider ±25 % ice thickness and up to 2× the ice softness, the possible range of the median deformation contribution is still small: 0.5-7 % in July and 0.5-8 % in November.If we consider only ±25 % ice thickness and do not change the ice softness, the range narrows to 0.5-4 % in both cases.Using stake measurements, Hodge (1974) estimated deforma- tion contributed ∼ 5-20 % of the velocity for the upper third of the ablation area of the Nisqually Glacier.He did not study any areas above the equilibrium line, so to compare directly to Hodge's (1974) numbers, we take the median deformation percentage over approximately the upper third of the ablation area and find a best estimate of 1 % (range 0.3-5 %) for July and 2 % (range 0.5-7 %) for November.These numbers suggest that sliding is even more dominant than Hodge (1974) estimated in this area, though it is difficult to say whether the differences are real (i.e., sliding was higher in 2012 than it was 4 decades ago) or just due to differences in methods and assumptions.
The model results for Emmons suggest that deformation is more important for the Emmons Glacier than for Nisqually.A median of 9 % of the July velocity field of Emmons can be explained by deformation, with a possible range of 3-40 % when considering ±25 % ice thickness and up to 2× the ice softness.When we consider only ±25 % ice thickness, the range narrows to 3-20 %.
There are a few regions where the observed surface velocity can be explained entirely or nearly entirely by internal deformation.These include the area within ∼ 1-2 km of the Nisqually and Emmons Glacier terminus, where ice is relatively thick and observed velocities are small.

Discussion
The continuous coverage of the TRI provides information about the spatial distribution of surface velocities.Several local velocity maxima are apparent along the centerline of the Nisqually Glacier and the central branch of the Emmons Glacier.These velocity maxima are associated with surface crevasses and increased surface slopes, with peak velocities typically observed just upstream of peak slope values (Fig. 6).They are likely related to accelerated flow downstream from local bedrock highs.
However, the local velocity maxima at ∼ 2.1 km in tions in local subglacial hydrology (e.g., reservoir drainage) during this time period.

Icefall and ice cliff dynamics
Terrestrial radar interferometry offers new observations over dynamic, inaccessible areas that have received limited attention in previous studies (e.g., icefalls, ice cliffs).For example, the velocities above the Nisqually ice cliff display an abrupt transition from slow-to fast-moving ice (Fig. 4).This rapid change from slow to fast favors crevasse opening and "detached slab" behavior rather than continuous flow, which is reflected in the heavily crevassed surface at this location.
Our results show that the Nisqually icefall and the icefall at the convergence of the Wilson and Nisqually glaciers show a slight increase in velocity from July to the winter months.This suggests that the icefalls may not be susceptible to the same processes that caused the seasonal velocity decrease over much of the rest of the glacier.This may indicate that there is a lack of local continuity through icefalls, which appears to prevent or dampen propagation of downstream seasonal velocity decreases.It could also indicate that the icefall is relatively well drained year round and is not significantly affected by seasonal changes in subglacial hydrology.A potential explanation for the observed minor increase in velocity could be early winter snow accumulation on blocks within the icefall.
Interestingly, in contrast to the icefall, the hanging glacier above the Nisqually ice cliff displayed a significant velocity decrease from July to November despite similar steep surface slopes and crevasse density.This could potentially be related to the lack of backstress from downstream ice and an increased sensitivity to minor fluctuations in subglacial hydrology.Hanging glaciers are also thought to be the source of some of the repeating glacial earthquakes that are triggered by snow loading (Allstadt and Malone, 2014), which highlights their sensitivity to minor perturbations.

Lack of significant diurnal variability
We expected to see significant variability over the 24-h July time series for Emmons, as atmospheric temperatures varied from 16 to 27 • C at Paradise Visitors Center (∼ 1600 m a.s.l), and skies remained cloud free during data collection.We hypothesized that the resulting increase in meltwater input from late morning through late afternoon might produce an observable increase in sliding velocity.While the results potentially show a slight velocity decrease at higher elevations overnight, and a slight velocity increase in the morning (Fig. 7a, A-D), these changes are not statistically significant nor coincident with times expected to have highest melt input.The lack of a significant diurnal speedup suggests that the subglacial conduits are relatively mature by July and are capable of accommodating the diurnal variations in meltwater flux without affecting basal sliding rates.
We did not expect to see significant diurnal changes in the 21-h November time series for Nisqually (Fig. 7a), as atmospheric temperatures ranged between 2 and 6 • C at Paradise Visitors Center (∼ 1600 m a.s.l.) and skies were partly cloudy to overcast during data collection, so surface meltwater input should have been minimal.Our results show only a minor velocity decrease higher on the glacier in the morning hours but it is not statistically significant and does not occur at times when we would expect increased meltwater.
Though some of the subtle changes in the extended time series may reflect actual diurnal velocity variability, we cannot interpret these with confidence.This suggests that the magnitude of diurnal variability, if it exists, during these time periods is minor when compared to the observed seasonal changes.It also implies that other stacks derived from a subset of the day can be considered representative of the daily mean and can be compared for seasonal analysis.

Seasonal velocity changes
The observed seasonal velocity changes from July to November can likely be attributed to changes in glacier sliding, which in turn are driven by evolving englacial and subglacial hydrology (Fountain and Walder, 1998).During the spring-summer months, runoff from precipitation (i.e., rain) and surface snow/ice melt enters surface crevasses, moulins, and/or conduits near the glacier margins.This water travels through a series of englacial fractures, reservoirs, and conduits and eventually ends up in a subglacial network of channels and reservoirs between the ice and bed.Storage time and discharge rates within the subglacial system are variable, with water finally exiting the system through one or more proglacial streams at the terminus.This dynamic system is continuously evolving due to variable input, storage capacity, and output.In early July, ongoing snowmelt should produce high meltwater discharge that travels through a relatively efficient network of mature conduits.As discharge decreases later in the summer, these subglacial conduits/reservoirs close due to ice creep without high flow to keep them open through melting due to heat from viscous dissipation.By November, there should be little or no surface meltwater input and we would expect to see a minimum in basal sliding velocity (Hodge, 1974).This is consistent with the observed velocity decrease in Fig. 4.However, the deformation modeling results (Fig. 9) show that a significant sliding component is still present for most of the Nisqually Glacier in November and December, when minimum surface velocities are expected.
The spatial patterns of the velocity change observed between July and November can be used to infer the extent of basal sliding.This may provide some insight into subglacial water storage, since the deformation component of surface velocity should remain nearly the same year round.Figure 4 indicates that almost the entire Nisqually Glacier slows down significantly between July and November, sug- gesting that storage is occurring under most of the glacier below the icefall and ice cliff.Significant velocity decreases are observed near local surface rises (Fig. 4), where some of the highest velocities were observed in July.This suggests that there are likely subglacial cavities downstream of these areas with high basal water pressures that can support enhanced sliding during the summer.Hodge (1974) interpreted a delay in both the maximum summer velocity and minimum winter velocity between the terminus and ELA as a propagating "seasonal wave" traveling ∼ 55 m day −1 .While our sampling is limited, the continued 2-27 November slowdown over the lower Nisqually near the terminus (Fig. 4f) could represent a delayed response to the significant slowdown over central Nisqually.This might be expected, as surface velocities near the terminus are dominated by internal deformation and should respond more slowly than areas dominated by basal sliding.

Comparison with historical velocity measurements
As described earlier, Hodge (1972Hodge ( , 1974) ) measured surface velocity for a network of centerline stakes on the lower Nisqually from 1968 to 1970.He documented a significant seasonal cycle with minimum velocities in November and maximum velocities in June.
To put our velocity data in historical context, we digitized Hodge's (1972) July and November 1969 surface velocity data at 19 stake locations along a profile of the lower half of the Nisqually Glacier.We then sampled the 2012 TRI slopeparallel velocities at these locations (Fig. 11).Remarkably, in spite of significant terminus retreat of up to ∼ 360 m and surface elevation changes of approximately −20 m (Sisson et al., 2011), the November 1969 and November 2012 surface velocities are almost identical at stakes 12-20, suggesting that bed properties and local geometry have greater influence over sliding velocity than ice thickness or relative distance from the terminus.In contrast, the July 2012 velocities at stakes 12-20 are 8-33 % faster than the July 1969 velocities.
The ice is mostly sliding at these locations, so the change could be related to a difference in the timing of the peak summer velocities or potentially enhanced sliding in 2012.The nearly identical surface velocities in November 1969 and 2012 suggest that the discrepancy between Hodge's sliding percentage estimates and our estimates (Sect.4.5) is likely related to different methodology and assumptions rather than actual changes in sliding since 1969.
The most notable difference between the profiles is observed closer to the terminus at stakes 7-12.At these locations, the July and November 2012 velocities are both < 0.05 m day −1 , whereas July and November 1969 velocities are ∼ 0.2 and ∼ 0.1 m day −1 , respectively, with significant seasonal variability.This suggests that the ice near the present-day terminus is essentially stagnant and no longer strongly influenced by changes in subglacial hydrology.

Conclusions
In this study, we used repeat TRI measurements to document spatially continuous velocities for numerous glaciers at Mount Rainier, WA, focusing primarily on the Emmons and Nisqually glaciers.We produced surface velocity maps that reveal speeds of > 1.0-1.5 m day −1 over the upper and central regions of these glaciers, < 0.2 m day −1 near the summit, and < 0.05 m day −1 over the stagnant ice near their termini.Novel data processing techniques reduced uncertainties to ±0.02-0.08 m day −1 , and the corrected, surface-parallel TRI velocities for Nisqually display similar magnitude and direction with a set of sparse interannual velocity measurements (Walkup et al., 2013).
Repeat surveys show that Nisqually Glacier surface velocities display significant seasonal variability.Most of the glacier experienced a ∼ 25-50 % velocity decrease (up to −0.7 m day −1 ) between July and November.These seasonal variations are most likely related to changes in basal sliding and subglacial water storage.Interestingly, the steep icefall displays no velocity change or even a slight velocity inwww.the-cryosphere.net/9/2219/2015/The Cryosphere, 9, 2219-2235, 2015 crease over the same time period.We documented no statistically significant diurnal velocity variations in ∼ 24-h data sets for Nisqually and Emmons, suggesting that subglacial networks efficiently handled diurnal meltwater input.Comparisons with 1969 velocity measurements over the Lower Nisqually (Hodge, 1972(Hodge, , 1974) ) reveal similar November velocities in both 2012 and 1969 and faster July velocities in 2012.
Using a simple 2-D ice flow model, we estimate that basal sliding is responsible for most of the observed surface velocity signal except in a few areas, mainly near the termini.The model suggests that about 99 % of the July velocity field for the Nisqually Glacier is due to sliding.Even when we account for the large uncertainties in ice thickness and ice softness, the possible range of sliding percentage is still narrow: 93-99.5 %.Deformation is more important for the Emmons Glacier, where we estimate 91 % of the observed motion is due to sliding, with a much wider possible range of 60-97 % when accounting for uncertainties.
In summary, TRI presents a powerful new tool for the study of alpine glacier dynamics.With just a few hours of fieldwork for each survey, we were able to document the dynamics of several glaciers at Mount Rainier in unprecedented extent and detail from up to 10 km away.TRI is particularly well suited for examining diurnal and seasonal glacier dynamics, especially for areas that are difficult to access directly (e.g., icefalls), like many parts of the glaciers at Mount Rainier.Repeat surveys provide precise surface displacement measurements with unprecedented spatial and temporal resolution, offering new insight into complex processes involving subglacial hydrology and basal sliding.Future studies involving coordinated, multi-day TRI occupations during critical seasonal transition periods could undoubtedly provide new insight into these and other important aspects of alpine glaciology.

Figure 1 .
Figure 1.Glaciers at Mount Rainier and locations of viewpoints used for ground-based radar interferometry.Instrument view angle ranges are indicated by arrows extending away from each viewpoint location.Boxes A-C show zoom areas for later figures.Inset map shows regional location of Mount Rainier.Glacier outlines in this and subsequent figures are from Robinson et al. (2010).

Figure 3 .
Figure 3. (a) Median slope-parallel velocity derived from TRI for GLPEEK and SUNRIZ viewpoints taken on 6-7 July 2011.(b) Width of 95 % confidence interval (high minus low limits for slope-parallel flow field) of slope-parallel velocities for 6-7 July 011 computed by bootstrapping after performing atmospheric noise corrections and stacking.Area shown is indicated by Box A on Fig. 1.

Figure 4 .
Figure 4. (a-d) Median slope-parallel velocities for Nisqually and Wilson glaciers for four different time periods taken from ROI viewpoint.Dashed lines on top left panel show locations of profiles taken to create Fig. 6, markers indicate distance in kilometers.(e-g) Percent change in median slope-parallel velocity for the Nisqually and Wilson glaciers between time periods.Blue indicates a velocity decrease and red indicates a velocity increase relative to the earlier time period; gray polygons indicate areas where velocity change is significant with 95 % confidence.Area shown is indicated by Box B on Fig. 1.

Figure 5 .
Figure 5. Width of 95 % confidence interval (high minus low limits for slope-parallel velocity) over Nisqually Glacier computed by bootstrapping.Shown for four sampling periods from the ROI viewpoint.Note that the color bar is scaled differently than Fig. 3b.

Figure 6 .
Figure 6.(a) and (c): slope-parallel velocity profiles along the two branches of Nisqually Glacier (profile lines shown in map view in Fig. 4a) for all sample time periods and viewpoints.(b) and (d): surface slope and ice thickness along each profile line.Surface slope is smoothed identically to that used for slope-parallel corrections (see text); ice thicknesses are estimated from digitized basal contours from Driedger and Kennard (1986) and surface elevations from the 2007/2008 LiDAR(Robinson et al., 2010).Refer to Fig.5and Table 2 for uncertainty estimates.

Figure 7 .
Figure 7. LOS velocity time series for areas outlined on maps to the right.Shaded region around each line represents ± 1 standard error for a 2-h running mean.(a) 24-h time series at SUNRIZ on 7-8 July 2012, gray box indicates the period with poor data quality (see text for details).(b) 22-h time series at ROI on 1-2 November 2012.

Figure 8 .
Figure 8.Comparison of average azimuth and velocities measured byWalkup et al. (2013) between 19 July and 11 October 2012 (black) compared to TRI slope-parallel velocities derived from this study at the same locations for two time periods that bracket the time period measured byWalkup et al. (2013).See Table4for comparison statistics and Box C in Fig.1for context.

Figure 9 .
Figure 9. Model results for summer (6 July 2012) and a late fall (2 November 2012) time period for Nisqually and Wilson glaciers.(a, d) Modeled surface velocity for internal deformation; (b, e) sliding residual (observed slope-parallel velocity minus the modeled deformation velocity); (c, f) estimate of the sliding percentage (sliding residual divided by total slope-parallel velocity).
Fig. 6   corresponds to a region of decreased surface slopes and increased ice thickness.This location also displayed significant seasonal velocity change, which could be related to varia-

Figure 11 .
Figure 11.July and November 1969 surface velocities measured byHodge (1974; digitized from Hodge, 1972)  at 19 stake locations along lower Nisqually profile (circles), compared with sampled 2012 slope-parallel velocities for corresponding locations/seasons (triangles).Stake locations are labeled and indicated with dotted lines and are shown in map view at right (same map extent as Fig.8).

Nisqually glacier Transmitting Antenna Receiving Antennae Instrument Controller 12V Battery & Inverter Radio- Frequency Unit Rotary Position Motor Figure 2. GPRI
equipment setup during the 27 November 2012 campaign at ROI viewpoint.

Table 2 .
Summary of uncertainty estimates of median stacks.Derived from bootstrapping, 95 % confidence, line of sight velocities; 2 correction refers to removing displacements due to atmospheric noise (interpolated over static control surfaces).

Table 3 .
Constants used in modeling analysis.

Table 4 .
Comparison between