Analyzing Glacier Surface Motion Using LiDAR Data

Understanding glacier motion is key to understanding how glaciers are growing, shrinking, and responding to changing environmental conditions. In situ observations are often difficult to collect and offer an analysis of glacier surface motion only at a few discrete points. Using light detection and ranging (LiDAR) data collected from surveys over six glaciers in Greenland and Antarctica, particle image velocimetry (PIV) was applied to temporally-spaced point clouds to detect and measure surface motion. The type and distribution of surface features, surface roughness, and spatial and temporal resolution of the data were all found to be important factors, which limited the use of PIV to four of the original six glaciers. The PIV results were found to be in good agreement with other, widely accepted, measurement techniques, including manual tracking and GPS, and offered a comprehensive distribution of velocity data points across glacier surfaces. For three glaciers in Taylor Valley, Antarctica, average velocities ranged from 0.8–2.1 m/year. For one glacier in Greenland, the average velocity was 22.1 m/day (8067 m/year).


Introduction
Glaciers are dynamic and in constant flux, and understanding glacial motion provides important information about their growth and retreat [1].However, the size and extent of glaciers, along with challenging environmental conditions, can severely limit data collection in the field or, in particularly hazardous conditions, field work may be precluded entirely [2][3][4].As remote sensing data becomes more widely available, it is becoming a primary data collection technique in the cryosphere, particularly in areas that are inaccessible to traditional field methods [2].
Traditional methods of collecting glacier velocity data rely on stakes drilled into the ice and/or GPS deployed on the glacier surface [5][6][7][8].The optical and geodetic surveys require little post-processing but only provide data at discrete points, usually covering a fraction of the entire glacier.Remote sensing techniques, including synthetic aperture radar (SAR) and multispectral imagery, are increasingly being used to detect and monitor changes in the cryosphere [4,[9][10][11][12].In difficult to reach areas of the Himalayas, SAR has been used to measure flow velocities at the Kangshung and Khumbu Glaciers [11].SAR was also used to measure flow velocities at the Shirase Glacier, Antarctica, and at Helheim Glacier, in Greenland [9,10,12].Often, these methods rely on satellite based platforms, which offer regular, long duration coverage, but they can be inhibited by severe topography and viewing angle [11,12].In addition, SAR tends to have low backscattering intensity over dry snow, reducing the volume of data collected in some environments [11,12].
Terrestrial and airborne laser scanning (TLS and ALS, respectively) are two active remote sensing techniques that use light detection and ranging (LiDAR), typically in the near-infrared wavelengths, to provide precise 3D elevation models of surfaces.TLS and ALS have both been used to survey glaciers to create highly detailed digital elevation models (DEMs) for the purpose of determining interannual variability in surface elevation in order to estimate mass balance or long-term volume change [2,3,13,14].
In this study, we use LiDAR DEMs, collected over time, to calculate glacier surface velocity using two different methods-particle image velocimetry (PIV) and manual tracking.Applying PIV, an image processing technique [15][16][17][18], to LiDAR data provides an opportunity to measure velocity across entire glaciers rapidly and the ability to measure high resolution nuances in the flow field that may be missed using other techniques.When compared to other methods such as feature extraction, and other image correlation techniques, PIV offers an approach that is sensitive to small scale, locally variable changes.This is especially important on glacier surfaces where common features may deform slightly between the collection of repeat data.Successful application of this method will make another technique available for assessing velocity of glaciers and other slow moving landforms.The PIV results are compared to manual tracking results and, where available, in situ data.The new method is tested on ALS and TLS point clouds of six glaciers in Antarctica and Greenland, with data covering a range of spatial and temporal resolutions.

Study Sites
The horizontal surface motion for six glaciers-five from Taylor Valley, Antarctica (Canada, Commonwealth, Rhone, Suess, Taylor) and one from Greenland (Helheim)-was analyzed in the study.Taylor Valley (Figure 1) is one of the McMurdo Dry Valleys (MDV) located in East Antarctica (77.5 • S, 163 • E).The valley landscape is composed of sandy, gravelly valley bottoms with expanses of exposed bedrock and is populated with perennially ice-covered lakes and ephemeral streams originating from glaciers that flow into the valleys from the surrounding mountains [19].Air temperatures average about −17 • C and summer temperatures typically hover just below freezing [20].
offer regular, long duration coverage, but they can be inhibited by severe topography and viewing angle [11,12].In addition, SAR tends to have low backscattering intensity over dry snow, reducing the volume of data collected in some environments [11,12].
Terrestrial and airborne laser scanning (TLS and ALS, respectively) are two active remote sensing techniques that use light detection and ranging (LiDAR), typically in the near-infrared wavelengths, to provide precise 3D elevation models of surfaces.TLS and ALS have both been used to survey glaciers to create highly detailed digital elevation models (DEMs) for the purpose of determining interannual variability in surface elevation in order to estimate mass balance or longterm volume change [2,3,13,14].
In this study, we use LiDAR DEMs, collected over time, to calculate glacier surface velocity using two different methods-particle image velocimetry (PIV) and manual tracking.Applying PIV, an image processing technique [15][16][17][18], to LiDAR data provides an opportunity to measure velocity across entire glaciers rapidly and the ability to measure high resolution nuances in the flow field that may be missed using other techniques.When compared to other methods such as feature extraction, and other image correlation techniques, PIV offers an approach that is sensitive to small scale, locally variable changes.This is especially important on glacier surfaces where common features may deform slightly between the collection of repeat data.Successful application of this method will make another technique available for assessing velocity of glaciers and other slow moving landforms.The PIV results are compared to manual tracking results and, where available, in situ data.The new method is tested on ALS and TLS point clouds of six glaciers in Antarctica and Greenland, with data covering a range of spatial and temporal resolutions.

Study Sites
The horizontal surface motion for six glaciers-five from Taylor Valley, Antarctica (Canada, Commonwealth, Rhone, Suess, Taylor) and one from Greenland (Helheim)-was analyzed in the study.Taylor Valley (Figure 1) is one of the McMurdo Dry Valleys (MDV) located in East Antarctica (77.5°S, 163°E).The valley landscape is composed of sandy, gravelly valley bottoms with expanses of exposed bedrock and is populated with perennially ice-covered lakes and ephemeral streams originating from glaciers that flow into the valleys from the surrounding mountains [19].Air temperatures average about −17 °C and summer temperatures typically hover just below freezing [20].The MDV receive very little precipitation annually (3-50 mm water-equivalent) [21] and the glaciers move only 0.3-20 m/year, though the fastest regions were not in any of the areas surveyed in this study [7,8].The terminus of Taylor Glacier, which is the area examined in this study, has been found to have surface velocities ranging from 0.3-10 m/year, based on surveys from 1993-2001 and 2002-2004 [7,8].Velocities at the Canada Glacier terminus were found to be slower, moving at a rate of up to 5 m/year and slowing rapidly near the glacier edges, based on the 1993-2001 survey [8].
Helheim Glacier (Figure 2) is an outlet glacier of the Greenland Ice Sheet (66.38 • N, 38.8 • W).Annual variations in temperature range between ±15 • C, with summer temperatures averaging around 5 • C [9].In contrast to glaciers in Taylor Valley, Helheim Glacier moves quite rapidly with velocities near the centerline of the terminus reaching up to 25 m/day and slowing up glacier (10-15 m/day) and near the glacier edges (0-5 m/day) [9,10,17].
Remote Sens. 2017, 9, 283 3 of 13 The MDV receive very little precipitation annually (3-50 mm water-equivalent) [21] and the glaciers move only 0.3-20 m/year, though the fastest regions were not in any of the areas surveyed in this study [7][8].The terminus of Taylor Glacier, which is the area examined in this study, has been found to have surface velocities ranging from 0.3-10 m/year, based on surveys from 1993-2001 and 2002-2004 [7-8].Velocities at the Canada Glacier terminus were found to be slower, moving at a rate of up to 5 m/year and slowing rapidly near the glacier edges, based on the 1993-2001 survey [8].
Helheim Glacier (Figure 2) is an outlet glacier of the Greenland Ice Sheet (66.38°N, 38.8°W).Annual variations in temperature range between ±15 °C, with summer temperatures averaging around 5 °C [9].In contrast to glaciers in Taylor Valley, Helheim Glacier moves quite rapidly with velocities near the centerline of the terminus reaching up to 25 m/day and slowing up glacier (10-15 m/day) and near the glacier edges (0-5 m/day) [9][10]17].

Methods
The Taylor Valley LiDAR point clouds were collected during two aerial surveys.NASA flew the aerial surveys during the austral summer of 2000-2001 and the data was processed by the US Geological Survey to create a detailed DEM of the valley bottom [22][23].NASA's Airborne Topographic Mapper (ATM), which has a green laser with a scan angle of ±15° and a scan frequency

Methods
The Taylor Valley LiDAR point clouds were collected during two aerial surveys.NASA flew the aerial surveys during the austral summer of 2000-2001 and the data was processed by the US Geological Survey to create a detailed DEM of the valley bottom [22,23].NASA's Airborne Topographic Mapper (ATM), which has a green laser with a scan angle of ±15 • and a scan frequency of 20 Hz, was used for data collection.The NASA point cloud density was, on average, one point per 2.7 m 2 [22,23].The region was resurveyed in the summer of 2014-2015 by the National Center for Airborne Laser Mapping (NCALM), University of Houston, under contract from Portland State University, for the purpose of estimating landscape change over the intervening 14 years since the NASA survey [24].The 2014 data were collected using a Teledyne Titan MW with three independent lasers operating at 532, 1064, and 1550 nm at an angle of ±30 • and a scan frequency of 20 Hz [25].An Optech Gemini ALTM, operating at 1064 nm with an angle of ±25 • and a frequency of 35 Hz, was also used for some areas that required higher altitude collections due to extreme topography.Point cloud density ranged from 2 points per m 2 to >10 points per m 2 and averaged 4.7 points per m 2 [26].
The Helheim Glacier surface was surveyed 23 times between 11 July and 12 July 2014 using a terrestrial laser scanner.The scanner was a RIEGL VZ-6000 operating at a wavelength of 1064 nm.For the purpose of this study, three TLS scans were selected for analysis.The first two scans were taken 51 min apart and the third scan was taken 12 h later, providing a test of fine time scale surface motion and overall surface change at Helheim.To examine the effect of varying point cloud density on the PIV analysis, two separate areas on Helheim glacier were studied in detail (Figure 2).Region 1 was located further from the scanner and had an average point cloud density of 1 point per m 2 .Region 2, located near the scanner, had an average point cloud density of 5 points per m 2 .GPS data collected coincidentally with the TLS data at multiple points on the glacier surface is available to verify the velocity measurements derived from PIV.
LiDAR point clouds cannot be used directly in PIV analysis, which is an image-based change detection technique [15][16][17].The point clouds were first converted to greyscale images and colored by elevation using the Terrascan software package within MicroStation.Interpolation in regions with little or no data was limited to three pixels and interpolated regions accounted for less than 10% of each study area; though data gaps larger than three pixels exist in some places, they were not interpolated.The interpolation was done by taking an average of the elevation values from the surrounding pixels.The image generation process grids the point cloud data and the resolution for this process was chosen based on the manual tracking results.The Taylor Valley glacier images were produced at a resolution of 2 m and 0.5 m was used for Helheim Glacier.In both cases, the resolution was finer than the motion being measured, based on the results of the manual tracking measurements.
Following the image conversion process, image pairs were analyzed using the PIVLab software, an open source MATLAB GUI [27][28][29].Particle image velocimetry cross-correlates subsections of two images by searching for sets of features common to both images to solve for the component velocity of each feature [15,16].When applied to glaciers, these trackable features are most commonly the peaks and valleys present on the glacier surfaces.PIVLab can complete multiple interpolation passes on a set of images at increasingly fine spatial resolution.The resolution of each pass was carefully chosen because a resolution smaller than the minimum displacement creates spurious results.A two pass interpolation was used; the resolutions of the first and second passes were scaled to be three times and two times, respectively, larger than the expected displacement, based on manual tracking results.For example, if the expected displacement was on the order of 1 m/year, the first pass would search the data in 3 m × 3 m grid cells and the second pass would use 2 m × 2 m grid cells.The grid resolution used to translate point clouds into images and the resolution of the PIV interpolation passes were two important factors in producing high resolution results and these values should always be chosen with respect to the temporal and spatial resolution of the data.
Given the time interval between images and the vector field of displacement, the glacier surface velocities were calculated.Raw coordinates and u-and v-velocity vectors were exported from PIVLab in a text file format and plotted in MATLAB.Rather than plot thousands of individual vectors, the raw results were gridded to show the average motion in a given region.The Taylor Valley glaciers were gridded in 200 m × 200 m cells and the Helheim Glacier results were gridded in 100 m × 100 m cells.The grid sizes were chosen such that the data in each cell was normally distributed and the standard deviation for each was calculated.
Manual tracking was also completed to validate the PIV results.Thirty to 40 points were chosen to measure displacement on each glacier by hand.These measurements were done directly on the point clouds in Microstation and only clearly defined sets of features, with unique identifying traits, were chosen for the manual tracking analysis.Displacement was measured either from peak to peak or trough to trough, depending on which had a higher density of LiDAR points.Error for the manual tracking results was calculated using the standard deviation.

Results
Commonwealth and Rhone Glaciers, in Taylor Valley, both presented unique challenges.Manual tracking could not be completed at either glacier.At Commonwealth Glacier, which has lower roughness on average (0.3 m) than any of the other glaciers examined, identifying features that can be used in manual tracking were typically smaller than the resolution of the NASA data (Figure 3a).Rhone Glacier presented the opposite challenge to manual tracking, with an average roughness of 0.7 m (Figure 3b).Though Taylor Glacier has a higher surface roughness, with an average of 1.7 m, valleys and peaks on Rhone are more chaotic and the lower resolution of the NASA data makes it difficult to identify the unique features that permit manual tracking with confidence.Commonwealth and Rhone Glaciers were excluded from further consideration with PIV because manual tracking measurements for reference were not available at either site.
Remote Sens. 2017, 9, 283 5 of 13 Manual tracking was also completed to validate the PIV results.Thirty to 40 points were chosen to measure displacement on each glacier by hand.These measurements were done directly on the point clouds in Microstation and only clearly defined sets of features, with unique identifying traits, were chosen for the manual tracking analysis.Displacement was measured either from peak to peak or trough to trough, depending on which had a higher density of LiDAR points.Error for the manual tracking results was calculated using the standard deviation.

Results
Commonwealth and Rhone Glaciers, in Taylor Valley, both presented unique challenges.Manual tracking could not be completed at either glacier.At Commonwealth Glacier, which has lower roughness on average (0.3 m) than any of the other glaciers examined, identifying features that can be used in manual tracking were typically smaller than the resolution of the NASA data (Figure 3a).Rhone Glacier presented the opposite challenge to manual tracking, with an average roughness of 0.7 m (Figure 3b).Though Taylor Glacier has a higher surface roughness, with an average of 1.7 m, valleys and peaks on Rhone are more chaotic and the lower resolution of the NASA data makes it difficult to identify the unique features that permit manual tracking with confidence.Commonwealth and Rhone Glaciers were excluded from further consideration with PIV because manual tracking measurements for reference were not available at either site.The flow field at Canada Glacier was highly coherent closer to the glacier terminus, where surface roughness increased (Figure 4).However, the flow field away from the rougher terminus region, towards the smoother ice up glacier showed a higher degree of scatter.The average velocity, as determined by PIV, on Canada Glacier was 1.41 ± 0.3 m/year.Manual tracking results yielded an average velocity of 1.62 ± 0.5 m/year, well within the uncertainty of the PIV results.When the outer 200 m of the glacier were considered separately from the center region, manual tracking produced faster velocities than PIV in both regions.The manual tracking center velocity estimate, 1.75 ± 0.6 m/year, fell within the error of the PIV results, 1.60 ± 0.4 m/year.The results for the edge region, 0.97 ± 0.1 m/year from PIV and 1.50 ± 0.3 m/year from manual tracking, did not agree within error though The flow field at Canada Glacier was highly coherent closer to the glacier terminus, where surface roughness increased (Figure 4).However, the flow field away from the rougher terminus region, towards the smoother ice up glacier showed a higher degree of scatter.The average velocity, as determined by PIV, on Canada Glacier was 1.41 ± 0.3 m/year.Manual tracking results yielded an average velocity of 1.62 ± 0.5 m/year, well within the uncertainty of the PIV results.When the outer 200 m of the glacier were considered separately from the center region, manual tracking produced faster velocities than PIV in both regions.The manual tracking center velocity estimate, 1.75 ± 0.6 m/year, fell within the error of the PIV results, 1.60 ± 0.4 m/year.The results for the edge region, 0.97 ± 0.1 m/year from PIV and 1.50 ± 0.3 m/year from manual tracking, did not agree within error though there were far fewer points to analyze in this region when compared to the center.The results for the Canada Glacier and the other Taylor Valley glaciers are reported in Table 1.
Remote Sens. 2017, 9, 283 6 of 13 there were far fewer points to analyze in this region when compared to the center.The results for the Canada Glacier and the other Taylor Valley glaciers are reported in Table 1.The results for Suess Glacier were more coherent in areas of higher roughness and showed more scatter in the smoother region in the center of the glacier (Figure 5).The average velocity at Suess Glacier was determined to be 0.86 ± 0.3 m/year, which is lower than the manual tracking average of 1.25 ± 0.7 m/year, though the two values agree within the stated error.Similar to the results for Canada Glacier, the manual tracking velocities for the glacier center and edge regions were faster than the PIV velocities in the same areas.The PIV center velocity, 1.00 ± 0.4 m/year, was significantly slower than the manual tracking velocity, 2.45 ± 0.3 m/year.However, there was agreement, within error, between the edge velocities (0.75 ± 0.2 m/year from PIV, and 1.03 ± 0.5 m/year from manual tracking).The results for Suess Glacier were more coherent in areas of higher roughness and showed more scatter in the smoother region in the center of the glacier (Figure 5).The average velocity at Suess Glacier was determined to be 0.86 ± 0.3 m/year, which is lower than the manual tracking average of 1.25 ± 0.7 m/year, though the two values agree within the stated error.Similar to the results for Canada Glacier, the manual tracking velocities for the glacier center and edge regions were faster than the PIV velocities in the same areas.The PIV center velocity, 1.00 ± 0.4 m/year, was significantly slower than the manual tracking velocity, 2.45 ± 0.3 m/year.However, there was agreement, within error, between the edge velocities (0.75 ± 0.2 m/year from PIV, and 1.03 ± 0.5 m/year from manual tracking).
Most of Taylor Glacier is too smooth to produce velocity results either through PIV or manual tracking of surface features.However, the region near the terminus is sufficiently rough to use both methods with success (Figure 6).The PIV results produced an average velocity of 1.98 ± 0.6 m/year.Manual tracking results produced an average velocity of 3.6 ± 1.3 m/year.The center and edge results from manual tracking (3.68 ± 1.0 m/year and 3.42 ± 1.3 m/year, respectively) were both faster than those for PIV (2.03 ± 0.9 m/year and 1.71 ± 0.6 m/year, respectively) but, in both cases, the results agreed within the error.Most of Taylor Glacier is too smooth to produce velocity results either through PIV or manual tracking of surface features.However, the region near the terminus is sufficiently rough to use both methods with success (Figure 6).The PIV results produced an average velocity of 1.98 ± 0.6 m/year.Manual tracking results produced an average velocity of 3.6 ± 1.3 m/year.The center and edge results from manual tracking (3.68 ± 1.0 m/year and 3.42 ± 1.3 m/year, respectively) were both faster than those for PIV (2.03 ± 0.9 m/year and 1.71 ± 0.6 m/year, respectively) but, in both cases, the results agreed within the error.Most of Taylor Glacier is too smooth to produce velocity results either through PIV or manual tracking of surface features.However, the region near the terminus is sufficiently rough to use both methods with success (Figure 6).The PIV results produced an average velocity of 1.98 ± 0.6 m/year.Manual tracking results produced an average velocity of 3.6 ± 1.3 m/year.The center and edge results from manual tracking (3.68 ± 1.0 m/year and 3.42 ± 1.3 m/year, respectively) were both faster than those for PIV (2.03 ± 0.9 m/year and 1.71 ± 0.6 m/year, respectively) but, in both cases, the results agreed within the error.The PIV results for Helheim Glacier show predominantly west to east motion on the glacier surface with a slight SE motion near the eastern boundary, which is less pronounced in the 51 min results (Figure 7) than the 12 h results (Figure 8).On average, PIV velocities ranged from 17.4 ± 0.5 The PIV results for Helheim Glacier show predominantly west to east motion on the glacier surface with a slight SE motion near the eastern boundary, which is less pronounced in the 51 min results (Figure 7) than the 12 h results (Figure 8).On average, PIV velocities ranged from 17.4 ± 0.5 -28.3 ± 2.5 m/day across the two regions and the two time frames examined (Table 2).These results were in good agreement with those from GPS, 20.7 ± 0.03-22.2± 0.16 m/day (Table 3), and manual tracking, 18.6 ± 3.1-33.5± 4.5 m/day (Table 2).The PIV results for Helheim Glacier show predominantly west to east motion on the glacier surface with a slight SE motion near the eastern boundary, which is less pronounced in the 51 min results (Figure 7) than the 12 h results (Figure 8).On average, PIV velocities ranged from 17.4 ± 0.5 -28.3 ± 2.5 m/day across the two regions and the two time frames examined (Table 2).These results were in good agreement with those from GPS, 20.7 ± 0.03-22.2± 0.16 m/day (Table 3), and manual tracking, 18.6 ± 3.1-33.5± 4.5 m/day (Table 2).

Discussion
The results from PIV on image pairs from four glaciers show good agreement with the results of other techniques, including manual tracking and, in the case of Helheim Glacier, GPS measurements.Overall, Canada Glacier showed the closest agreement between the manual tracking and PIV results, with an average difference of just 0.21 m/year between the two measurements.This is most likely due to the well distributed across-glacier peaks and valleys present on Canada Glacier.Importantly, these features are more present even in areas of low roughness than on either Suess or Taylor Glaciers.The manual tracking and PIV results agreed, within error, in every case except for the edge region of Canada Glacier and the center region of Suess Glacier.These two discrepancies may be the result of bias inherent in manual tracking since the method requires features to be clear and well defined in order to be accurately tracked.Very smooth or very rough regions become difficult to track manually with confidence but the PIV algorithm, which tracks brightness as a function of elevation, is able to track more subtle patterns of motion.
Despite the good agreement between PIV and manual tracking overall, PIV derived velocities were consistently lower than those determined through manual tracking.Particle image velocimetry is more sensitive to smaller glacier motions than manual tracking.Glaciers do not move at one speed across their entire surface and PIV can more easily track smaller, slower features that appear when an area of motion is examined as a whole but are not significant enough to detect when looking for individual features [7,8,17].As an additional method of verification, GPS data collected at ablation stakes on glacier surfaces can be used where they are available.
Along with the GPS data collected at Helheim Glacier (Table 3), GPS data are also available for Canada and Taylor Glaciers [30].The time period represented in this data is from 1995 to 2001 so it is not contemporaneous with the LiDAR data collection but it is still a valuable point of comparison.Canada Glacier was found to be moving at an average of 0.77 ± 0.45 m/year over the whole glacier surface and at a slower rate of 0.45 ± 0.12 m/year near the terminus, which agrees well with the average and edge velocities determined using PIV.Taylor Glacier showed an average velocity of 1.0 ± 0.35 m/year with a velocity near the terminus of 1.11 ± 0.21 m/year, which also shows good agreement with the PIV results.Published averages for Taylor Glacier as a whole tend to be higher than these estimates, ranging from 0.3-25 m/year [6][7][8] but maps of Taylor Glacier published by Fountain et al. (2006) and Kavanaugh et al. (2009) show the highly variable nature of velocity across the glacier, with higher velocities located outside of the area included in this study.
In situ data was not available for Suess Glacier.The results from PIV and manual tracking agreed for the average and edge glacier velocities but not for the center region.Suess Glacier's surface characteristics, seen more clearly at a resolution of 100 m × 100 m (Figure 9), offer some clues to this discrepancy.The SW center region has a higher roughness, and correspondingly clearer surface features, than the NE center region.There were relatively few features available for manual tracking in the NE region, which the PIV results show to be moving significantly slower than the SW center.The manual tracking results for Suess Glacier are, consequently, biased towards higher velocities.At a resolution of 100 m × 100 m, the distribution of data for Suess Glacier is still predominantly Gaussian.Though the 200 m × 200 m resolution results were used to compute the statistics shown in Table 1, the finer resolution results, coupled with the glacier roughness, more fully describe the reason for the discrepancy between PIV and manual tracking velocities in the glacier center.
Surface roughness and the shape of surface features both play a key role in the utility of PIV and manual tracking to analyze glacier LiDAR data.Taylor Glacier has a higher surface roughness than any of the other glaciers considered in the study, with an average surface roughness of 1.7 m.However, PIV and manual tracking both worked well on Taylor Glacier but not on Rhone Glacier, with an average surface roughness of 0.7 m.While surface roughness provides a useful first pass test for interpreting PIV results and understanding possible sources of error, the nature of the surface features themselves are also important.Features on Taylor Glacier are more easily defined in a visual analysis and they have distinct shapes in both the NASA and NCALM LiDAR point clouds.Despite having lower geometric roughness, features on the surface of Rhone Glacier are more closely packed and less distinct.Surface roughness and the shape of surface features both play a key role in the utility of PIV and manual tracking to analyze glacier LiDAR data.Taylor Glacier has a higher surface roughness than any of the other glaciers considered in the study, with an average surface roughness of 1.7 m.However, PIV and manual tracking both worked well on Taylor Glacier but not on Rhone Glacier, with an average surface roughness of 0.7 m.While surface roughness provides a useful first pass test for interpreting PIV results and understanding possible sources of error, the nature of the surface features themselves are also important.Features on Taylor Glacier are more easily defined in a visual analysis and they have distinct shapes in both the NASA and NCALM LiDAR point clouds.Despite having lower geometric roughness, features on the surface of Rhone Glacier are more closely packed and less distinct.
A comparison of Helheim Glacier to the Taylor Valley glaciers shows that spatial and temporal resolution is also vital to the use of PIV and manual tracking.Spatially, on average, all three of the Helheim Glacier TLS scans used here have higher resolution than the Antarctica ALS scans.While the NCALM data collected in 2014 has a similar average resolution, the earlier NASA data collected in 2001 has a substantially lower resolution on average (<1 point/m 2 compared to ~5 points/m 2 ).This significantly limits the size of detectable features.Additionally, while large features on the Taylor Valley glaciers remain easily recognizable, smaller features could be obscured by changes on the glacier surface.The high frequency and resolution of the Helheim Glacier TLS scans make even very small features trackable manually and digitally.

Conclusions
Particle image velocimetry has been shown to be a useful technique for analyzing glacier surface velocity using remotely collected LiDAR data.Four glaciers-Helheim, Canada, Suess, and Taylorwere analyzed in this study using PIV, manual tracking, and, where available, GPS.Strong agreement was found between average PIV and manual tracking velocities at all of the glaciers and these findings were verified using GPS data at Canada, Taylor, and Helheim Glaciers.Most of the glaciers also showed strong agreement between PIV and manual tracking velocities when examining the glacier center and edges separately, a strong indicator that PIV can help to determine full, acrossglacier velocity fields, rather than providing a few discrete data points.A comparison of Helheim Glacier to the Taylor Valley glaciers shows that spatial and temporal resolution is also vital to the use of PIV and manual tracking.Spatially, on average, all three of the Helheim Glacier TLS scans used here have higher resolution than the Antarctica ALS scans.While the NCALM data collected in 2014 has a similar average resolution, the earlier NASA data collected in 2001 has a substantially lower resolution on average (<1 point/m 2 compared to ~5 points/m 2 ).This significantly limits the size of detectable features.Additionally, while large features on the Taylor Valley glaciers remain easily recognizable, smaller features could be obscured by changes on the glacier surface.The high frequency and resolution of the Helheim Glacier TLS scans make even very small features trackable manually and digitally.

Conclusions
Particle image velocimetry has been shown to be a useful technique for analyzing glacier surface velocity using remotely collected LiDAR data.Four glaciers-Helheim, Canada, Suess, and Taylor-were analyzed in this study using PIV, manual tracking, and, where available, GPS.Strong agreement was found between average PIV and manual tracking velocities at all of the glaciers and these findings were verified using GPS data at Canada, Taylor, and Helheim Glaciers.Most of the glaciers also showed strong agreement between PIV and manual tracking velocities when examining the glacier center and edges separately, a strong indicator that PIV can help to determine full, across-glacier velocity fields, rather than providing a few discrete data points.
Surface roughness, surface feature shapes and distribution, LiDAR point cloud density, and temporal resolution were all found to be important when applying PIV to glacier surface motion.While PIV can be applied to glaciers that are difficult to examine using manual tracking (Commonwealth, Rhone), there is no way to check that the results are reasonable.In addition, the PIV results for both Commonwealth and Rhone Glaciers showed a great deal more scatter in the data than at any of the other four glaciers used in the study, suggesting that the same surface feature geometry that inhibits manual tracking also inhibits PIV.However, on glaciers with distinct features, PIV offers a method to rapidly and systematically collect thousands of velocity vectors, characterizing glacier motion.Ideally, PIV should be complementary to other types of glacier velocity data.Regardless though, it has the potential to significantly expand our understanding of the complex dynamics that characterize glacier flows by offering a wider, more comprehensive data distribution than other currently available techniques.
Whether these techniques are used to measure glacier flow or other moving landforms, special attention needs to be paid to the time scale of data collection and the spatial scales used to rasterize and analyze LiDAR data.At the data collection stage, LiDAR data acquisitions must be temporally spaced to allow movement beyond the uncertainty of the data while still maintaining common features.For example, a rapid warming event or large volume of snowfall could significantly change the appropriate timescale for repeat measurements, since both events have the potential to occlude glacier surface features.Once data is collected, the key to rasterizing for use with PIV is to choose a scale that is larger than the uncertainty but smaller than the expected motion.Finally, the interpolation windows selected for PIV analysis should be larger than the anticipated motion so that faster moving regions will not be falsely truncated.These guidelines should be applicable to LiDAR-based PIV studies across a number of disciplines.

Figure 2 .
Figure 2. Helheim Glacier near the terminus.The blue triangle indicates the TLS scan position, the orange diamonds indicate the position of the two GPS units on the glacier surface during the TLS data collection, and black boxes indicate the bounds of study regions 1 and 2. The image was retrieved from the online Earth Explorer, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https://earthexplorer.usgs.gov/.

Figure 2 .
Figure 2. Helheim Glacier near the terminus.The blue triangle indicates the TLS scan position, the orange diamonds indicate the position of the two GPS units on the glacier surface during the TLS data collection, and black boxes indicate the bounds of study regions 1 and 2. The image was retrieved from the online Earth Explorer, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https://earthexplorer.usgs.gov/.

Figure 3 .
Figure 3. (a) Commonwealth Glacier roughness from the 2014 NCALM data is shown in the left figure.The rougher (light blue) area indicates the extent of the glacier.White areas contain no data and are mostly located on the valley floor, rather than the glacier surface.(b) Rhone Glacier roughness from the 2014 NCALM data is shown on the right.The extent of the glacier is still fairly clearly described by the higher values of roughness (light blue).

Figure 3 .
Figure 3. (a) Commonwealth Glacier roughness from the 2014 NCALM data is shown in the left figure.The rougher (light blue) area indicates the extent of the glacier.White areas contain no data and are mostly located on the valley floor, rather than the glacier surface.(b) Rhone Glacier roughness from the 2014 NCALM data is shown on the right.The extent of the glacier is still fairly clearly described by the higher values of roughness (light blue).

Figure 4 .
Figure 4. Canada Glacier PIV Results.Surface roughness (average 0.4 m) of Canada Glacier is overlaid by PIV results, averaged over 200 m × 200 m cells.

Figure 4 .
Figure 4. Canada Glacier PIV Results.Surface roughness (average 0.4 m) of Canada Glacier is overlaid by PIV results, averaged over 200 m × 200 m cells.

Figure 5 .
Figure 5. Suess Glacier PIV results for 200 m × 200 m gridded data.The color scale shows the glacier surface roughness (average 0.4 m).Glacier surface velocities are higher in the center of the glacier and decrease towards the edges.

Figure 6 .
Figure 6.PIV derived surface flow field at the Taylor Glacier terminus overlaid on a map of surface roughness (average 1.7 m).The raw results were averaged over 200 m × 200 m regions.

Figure 5 .
Figure 5. Suess Glacier PIV results for 200 m × 200 m gridded data.The color scale shows the glacier surface roughness (average 0.4 m).Glacier surface velocities are higher in the center of the glacier and decrease towards the edges.

Figure 5 .
Figure 5. Suess Glacier PIV results for 200 m × 200 m gridded data.The color scale shows the glacier surface roughness (average 0.4 m).Glacier surface velocities are higher in the center of the glacier and decrease towards the edges.

Figure 6 .
Figure 6.PIV derived surface flow field at the Taylor Glacier terminus overlaid on a map of surface roughness (average 1.7 m).The raw results were averaged over 200 m × 200 m regions.

Figure 6 .
Figure 6.PIV derived surface flow field at the Taylor Glacier terminus overlaid on a map of surface roughness (average 1.7 m).The raw results were averaged over 200 m × 200 m regions.

Figure 7 .
Figure 7. Region 1 results on the left and Region 2 results on the right, both for a time interval of 51 minutes.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m, respectively.

Figure 8 .
Figure 8. Region 1 results on the left and Region 2 results on the right, both for a time interval of 12 hours.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m respectively.

Figure 7 .
Figure 7. Region 1 results on the left and Region 2 results on the right, both for a time interval of 51 min.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m, respectively.

Figure 7 .
Figure 7. Region 1 results on the left and Region 2 results on the right, both for a time interval of 51 minutes.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m, respectively.

Figure 8 .
Figure 8. Region 1 results on the left and Region 2 results on the right, both for a time interval of 12 hours.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m respectively.

Figure 8 .
Figure 8. Region 1 results on the left and Region 2 results on the right, both for a time interval of 12 h.The Helheim Glacier PIV results have been plotted on top of a contour plot of point density, rather than roughness.It is important to note that point density is scaled differently for Regions 1 and 2, in order to highlight the attributes of each.Raw PIV results were averaged at a resolution of 100 m × 100 m.Though not shown, average surface roughnesses for Regions 1 and 2 were 0.8 m and 1.0 m respectively.

Figure 9 .
Figure 9. Suess Glacier PIV results for 100 m × 100 m gridded data.The color scale shows the glacier surface roughness (average 0.4 m).Glacier surface velocities are higher in the center of the glacier and decrease towards the edges.

Figure 9 .
Figure 9. Suess Glacier PIV results for 100 m × 100 m gridded data.The color scale shows the glacier surface roughness (average 0.4 m).Glacier surface velocities are higher in the center of the glacier and decrease towards the edges.

Table 1 .
Results for Canada, Suess, and Taylor Glaciers in Taylor Valley.All refers to the entire glacier surface that was used in the PIV analysis.Edge refers to the 200 m of glacier nearest to the glacier edge and center refers to everything other than the edge.Glacier Region PIV (m/year) Error Manual Tracking (m/year) Error

Table 1 .
Results for Canada, Suess, and Taylor Glaciers in Taylor Valley.All refers to the entire glacier surface that was used in the PIV analysis.Edge refers to the 200 m of glacier nearest to the glacier edge and center refers to everything other than the edge.

Table 2 .
PIV and manual tracking results for Helheim Glacier.

Table 2 .
PIV and manual tracking results for Helheim Glacier.

Table 2 .
PIV and manual tracking results for Helheim Glacier.