Extracting recent short-term glacier velocity evolution over Southern Alaska from a large collection of Landsat data

Southern Alaska from a large collection of Landsat data Bas Altena1, Ted Scambos2, Mark Fahnestock3, and Andreas Kääb1 1Department of Geosciences, University of Oslo, Blindern, 0316 Oslo, Norway 2National Snow and Ice Data Center (NSIDC), University of Colorado, Boulder, CO 80303, USA 3Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK 99775, USA Correspondence to: Bas Altena (bas.altena@geo.uio.no)

Glacier velocity monitoring through satellite remote sensing has proven to be a useful tool to observe velocity change on a basin scale.Several studies have focused on dynamics of individual glaciers in Alaska, at an annual or seasonal resolution (Fatland and Lingle, 2002;Burgess et al., 2012;Turrin et al., 2013;Abe and Furuya, 2015;Abe et al., 2016).Such studies can give a better understanding of the specific characteristics of a glacier, and which circumstances are of importance for this behaviour and response.Region-wide annual or "snapshot" velocities also have been estimated over the Saint Elias Mountain range in previous studies (Burgess et al., 2013;Waechter et al., 2015).Their results give a first-order estimate of the kinematics at hand.
With frequent satellite data coverage, it is possible to detect the time of glacier speed-ups to within a week (Altena and Kääb, 2017b), although the study did not include an automated approach.In the most recent work, regional analyses have been conducted with over sub-seasonal (Moon et al., 2014;Armstrong et al., 2017) and multi-decadal (Heid and Kääb, 2012a;Dehecq et al., 2015) periods.With such data one is able to observe the behaviour of groups of glaciers that experience similar climatic settings.Consequently, surges and other glacier-dynamical events can be put into a wider spatio-temporal perspective.
Since the launch of Landsat 8 in 2013 a wealth of high-quality medium-resolution imagery is being acquired over the cryosphere on a global scale.Onboard data storage and rapid ground-system processing have made it possible to almost continuously acquire imagery.The archived data has enormous potential to advance our knowledge in glacier flow.Extraction of glacier velocity is one of the stated mission objectives (Roy et al., 2014).However, the data rate far exceeds the possibilities for manual interpretation.Fortunately, automatically generated velocity products are now available (Scambos et al., 2016;Rosenau et al., 2016), though at this point sophisticated quality control and post processing methods are still being developed.
Up to now, most studies of glacial velocity have had an emphasis on either spatial or temporal detail.When temporal detail is present, studies focus on a single or a handful of glaciers (Scherler et al., 2008;Quincey et al., 2011;Paul et al., 2017).On the other hand, when regional assessments are the focus, the temporal resolution ranges from a single time stamp up to annual resolution (Copland et al., 2009;Dehecq et al., 2015;Rosenau et al., 2015).Furthermore, most studies rely on filtering in the post-processing of vector data by using the correlation (Scambos et al., 1992;Kääb and Vollmer, 2000) or through median filtering within a zonal neighborhood (Skvarca, 1994;Paul et al., 2015).Some sophisticated post-processing procedures are available (Maksymiuk et al., 2016), but rely on the coupling with models based upon the Navier-Stokes equations.Also geometric properties can be taken into account during the matching to improve robustness and reduce post-processing efforts, such as reverse-correlation (Scambos et al., 1992;Jeong et al., 2017) or triangle closure (Altena and Kääb, 2017a).
Thus glacier velocity data is increasingly available, but in general post-processing is not at a sufficient level to directly exploit the full information content within these products.In this study we aim at the construction of a post-processing chain that is capable of extracting temporal information from stacks of noisy velocity data.Our emphasis is on discovering patterns over a mountain-range scale.Analysis of the individual details of the glacier-dynamical patterns identified by the processing will be considered in later work.For a single glacier, a manual selection of low-noise, good-coverage velocity data sets is possible.
However, such a strategy will not be efficient when multiple glaciers or mountain ranges are of interest.In this contribution we want to develop the methodological possibilities further and try to construct glacier velocities at a monthly resolution.Therefore, our implementation focuses on automatic post-processing, without the help of expert knowledge or human interaction.
The presented method retains spatial detail present in the data and does not simplify the flow structure to flowlines.Consequently, with this methodology we want to generate products that can improve our knowledge about the influence and timing of tributary and neighbouring ice flow variations.
In this study, we discuss the data used and provide background on the area under study.We then introduce the spatiotemporal structure of the data, followed by an explaination of our process for vector "voting" and vector field smoothing.The next section highlights our results and our validation and assessment of the performance of our method.
2 Data and study region

GoLIVE velocity fields
The velocity fields used in this study are based upon repeat optical remote sensing imagery and are distributed through the National Snow and Ice Data Center (NSIDC) (Scambos et al., 2016).These velocity fields are derived from finding displacements between pairs of Landsat 8 imagery, using the panchromatic band with 15 meter resolution.A high-pass filter of one kilometer spatial scale is applied before processing.Normalized cross-correlation is applied between the image pairs on a sampling grid with 300 meters spacing (Fahnestock et al., 2016;Scambos et al., 1992) and a template size of 20 pixels (or 300 m).The resulting products are grids with lateral displacements, the absolute correlation value, signal-to-noise ratio and ratio between the two best matches.At the time of writing, displacement products can cover a time interval from 16 days up to 96 days.For a detailed description of the processing chain see Fahnestock et al. (2016).
The Landsat 8 satellite has a same-orbit revisit time of 16 days and a swath width of 185 km.Only scenes which are at least 50% cloud-free are used (as determined by the provided estimate in the metadata for the scenes).Consequently, not every theoretical pair combination is matched, and also pairs across tracks are neglected (cf.(Altena and Kääb, 2017a)) to avoid more complicated viewing geometry adjustments.Georeference errors are compensated by the estimation of a polynomial bias surface through areas outside glaciers (i.e., assumed stable).The glacier mask used for that purpose is from the Randolph Glacier Inventory (RGI) (Pfeffer et al., 2014).The resulting grids come in Universal Transverse Mercator (UTM) projection and if orthorectification errors are minimal, displacements for precise georeferencing require only horizontal movement of a few meters (generally <10 m).In total we use twelve Landsat path/row tiles to cover our study area (Figure 1). -145

Study region
The region of interest covers the Saint Elias, Wrangell and Kluane Mountain ranges, as well as some parts of the Chugach range.These ranges host roughly 42 000 km 2 glacier area, whereby roughly 22% of the glacier area is connected to marine terminating fronts draining into the Gulf of Alaska (Molnia, 2008).The glaciers in this area are diverse, as a wide range of thermal conditions (cold and warm ice) and morphological glacier types (valley, icefields, marine terminating) occur in these mountain ranges (Clarke and Holdsworth, 2002).This might partly be due to the large precipitation gradient over the mountain range.
The highest amount of precipitation falls in summer or autumn.The study area covers mountain ranges that have two different clusters of climate.Along the coast one finds a maritime climate with a small annual temperature range.These mountains function as a barrier, and the mountain ranges behind, in the interior, have therefore a more continental climate (Bieniek et al., 2012).

Methodology
GoLIVE and other velocity products are composed products from at least two satellite acquisitions.When imagery from multiple time instances are used, combinations of displacements, with different (overlapping) time intervals can be constructed.In order to be of use for time-serie analysis, such velocity fields with different time spans need to be combined into a dataset with regular timesteps.To reduce the noise, the temporal configuration of these products can be used to synthesize a multi-temporal   velocity field.

Temporal network configuration
At the high latitude of the Southern Alaska scenes from adjacent tracks have an overlap of 60%.Looking at only one track (or satellite path), the 16-day revisit makes several matching combinations of integer multitudes of 16 days possible.For example, over a 64-days period (∆t), five images are acquired in one satellite track and their potential pairing combinations can be illustrated as a network (Figure 2).In this network, every acquisition (at time t) is a node, and these nodes are connected through an edge that represents a matched pair leading to a collection of displacements (d) with an associated similarity measures (ρ).
When velocities over different timespans are estimated, this network has in theory a great amount of redundancy.However in practice this is complicated, as combinations of images are not processed when there is too much obstruction by clouds.Furthermore, individual displacements can be gross errors, as an image match was not established due to surface change or lack of contrast and thus loss of similarity.Consequently, when data from such a network is combined to synthesize one consistent velocity time series the estimation procedure needs to be able to resist multiple outliers or be able to identify whether displacement estimates could be extracted at a reliable level at all.The network shown in figure 2 can be seen as a graph; nodes correspond to timestamps and edges to matched image pairs.Such a graph can be transformed into an adjacency matrix (A G , see Figure 2).In this matrix the columns and rows represent different timestamps.The edges can be directed, meaning it can assign which acquisition is the master (reference) or the slave (search) image during the matching procedure.For the GoLIVE data, the oldest acquisition is always the reference image, hence within the matrix only the upper triangular part has filled entities.The spacing of the timesteps is 16 days and the amount of days is set into the corresponding entries when a time step is covered by an edge.Individual days are specified instead of a binary entity, to be able to merge adjacency matrices from different tracks which have different acquiring dates.If partial overlap of an edge occurs, then the time steps are proportionally distributed.For example, for a small network of three nodes, velocity (v) can be estimated through least-square adjustment of the displacements (d) through the following systems of equations (Altena and Kääb, 2017a), Here the subscript denotes the timespan given by the starting and ending timestamps of the interval.This relational structure of displacements is similar to a leveling network.When the adjacency matrix is converted to an incidence matrix, then this matrix is the design matrix (A) (Strang and Borre, 1997).This makes the generation of such network adjustments easily implemented.
This formulation of the temporal network makes it possible to estimate the unknown parameters, i.e. the temporal components of the velocity time series, through different formulations.This is illustrated in Figure 2, where a selection of a velocity (between t 2 & t 3 ) is estimated.Displacements that fall between the two images can be used for the estimation (here blue), which we here call a "closed" network.But as can be seen in the figure as red connections, other displacements from outside the time span are over-arching and stretching further than the initial time interval.Such measurements can be of interest as they can fill in gaps or add redundancy, but the glacier dynamics obtained will be smeared compared to the real ones.Consequently, we call such a network configuration an "open" network (here red).

Voting
The velocity dataset we use (like any) contains a large amount of incorrect or noisy displacements.Typically, the distribution of displacements has a normal distribution but with long tails.Moreover, a least-square adjustment is very sensitive to outliers contained in the data to be fitted.Therefore, direct least-square computation of velocity through the above network is not easily possible and some selection procedure is needed to exclude gross errors.Outlier detection within a network such as in equation 1 can be done through statistical testing (Baarda, 1968;Teunissen, 2000), assuming measurements (d) are normally distributed.However, such procedures are less effective when several gross errors are present within the set of observations.
Extracting information from highly contaminated data is therefore an active field of research.For example, robust estimators change the normal distribution to a heavy tailed distribution.Nevertheless, such estimations typically still start with normal least-squares adjustment based on the full initial set of observations, and only in the next step the weights are iteratively adjusted according to the amount of misfit.Hence, such methods are still restricted to robust a-priori knowledge or a data-set with relatively small amounts of contamination by gross errors.
Another common approach to cope with the adjustment of error-rich observations is through sampling strategies such as least-median of squares (Rousseeuw and Leroy, 2005), or random sampling and consensus (RANSAC) (Fischler and Bolles, 1981).A minimum amount of observations are picked randomly to solve the model.The estimated parameters are then used to assess how the initial model fits in respect to all observations.Then the procedure is repeated with a new set of observations.
The sampling procedure is stopped when a solution is within predefined bounds, or executed a defined amount of times after which the best set is taken.Such methods are very popular as they can handle high contamination of data (up to 50%) and still result in a correct estimate.Put differently, the break-down point is .5 (Rousseeuw and Leroy, 2005).However, we use a different approach as these methods implement polynomial models.Our data set benefits from including conditional equations as well.
The equations that form the model can be seen as individual samples that populate the parameter space.In such a way the individual relations within the equation propagate into points, lines, or surfaces depending on their dimensions and relation given by the equation.Hence, measurements can be transformed into a shape that is situated within the parameter space (which has a finite extent and resolution).The collection of shapes will be scattered through out this parameter space, but such shape converge at a common point which is most likely the correct parameter values.This transformation is the Hough transform and is commonly used in image processing for the detection of lines and circles.
However, the Hough transform is discrete while the measurements are noisy, thus the shapes should also be fuzzy.In this study we follow this latter direction and discretize the displacement-matching search-space after which we exploit a voting strategy.But the shapes in parameter space are blurred, so a fuzzy Hough transform (Han et al., 1994) is implemented.Our matching search space is simply the linear system of equations of the network described above.To illustrate the system, an example of a network with four displacements (d) is shown in Figure 3.These displacements can be placed within a discrete grid that spans a large part of the parameter space (in this example two dimensional).For the short time span displacements, this results in vertical or horizontal lines (blue, purple and yellow).The lines originate from the ambiguity, which is also in the entries of the design matrix (A) as these are 0 entries.For the overarching displacement (in red) the displacement is a combination of both velocities.Only the total displacement is known, thus the line results in a diagonal orientation (5 + 2 can be a solution as well as 1 + 6, or a non-integer combination).
In this toy example, the lines concentrate around 3 for both velocities, with a slight preference to the upperleft, which indicates a slight speed-up between these two timespans.However, the lines are not perfectly intersecting, as measurement noise is present.To simulate this noise, the lines are perturbed with a Gaussian function, as for the toy example shown in figure 4a.
In this case, every observation will fill the parameter space with a discretized weighting function.This fuzziness in the Hough space makes it possible to find the intersection, while noise is present.Furthermore, the dimension of the Hough space dependents on the formulation of the network.In principle it can have any dimension; in one dimensional case this is a simple histogram, but in higher dimensions this will translate into a line which radially decreases in weight.For this study we implement a three-dimensional Hough space, which for our toy example will look like as in figure 4b, though due to visualization limitations, the fuzzy borders are not visualized.
The advantage of a Hough search space is the resistance to multiple outliers.It builds support and is not reliant on the whole group of observations.Especially, when a second or third dimensional space is used, the chances of random (line)crossing decrease significantly in parameter space.Hence, such events will stand out when multiple measurements do align.Furthermore, random measurement errors can be incorporated through introducing a distribution function.In our implementation this is a Gaussian, but other functions are possible as well.The disadvantage of the fuzzy Hough transform is its limitation to implement a large and detailed search space, as the dimension and resolution depend on the available computing resources.
The fuzzy Hough transform functions as a selection process to find observations which are to a certain extent in agreement.
With this selection of inliers the velocity can finally be estimated through ordinary least-squares estimation.The model is the same as used to construct the network.However, the observations without consensus (i.e.outliers) are not used.The remaining observations can, nevertheless, still be misfits, such as from shadow casting, as no ice flow behavior is prescribed in the design matrix of equation 1.

Smoothing
Because the voting and least-squares adjustment in our implementation has no neighborhood constraints but is rather strictly per matching grid point, the velocity estimates contain systematic, gross and random errors, though reduced with respect to the initial data set.This least squares adjustment with voted displacements results in a spatio-temporal stack of velocity estimates that have a regular temporal spacing.However, due to undersampling as a result of cloud cover, the stack might have holes or ill-constrained estimates.We apply a spatial-temporal smoothing taking both spatial and temporal information into account using the Whitacker approach that tries to minimize the following function (S), ( This formulation is the one dimensional case, where for every location (x, where i denotes the index of the grid) an estimate (denoted by •) is searched for that minimizes S.Here ∆ denotes the difference operator, thus ∆x i = x i+1 − x i .Similarly, ∆ 2 is the double difference, describing the curvature of a signal (∆ . For the implementation of this method we use the procedure presented by Garcia (2010).This routine has an automatic procedure to estimate the smoothing parameter (λ) and has robust adaptive weighting (w).Its implementation is conducted through a Discrete Cosine Transform (DCT), which eases the computational load.Furthermore, a discrete cosine transform operates both globally and locally, and in multiple dimensions.Lastly, in order to include all data at once, the vector field is configured as a complex number field.In order to make this isotropic, the vertical axis needs to extent, so the spread is similar in both directions.
The smoothing parameter (λ) is operating over both space (2 dimensions) and time (1 dimension), but the smoothing parameter is a single scalar.Hence in this form it would be dependent on the choice of grid resolutions in time and space.Therefore, in order to get rid of this dependency and fulfill the isotropy property, the spatial and temporal dimensions are scaled.For this scaling estimation we construct an experimental variogram and look at its distribution (Wackernagel, 2013).Along the spatial axis, the variogram in Figure 5a shows spatial correlation up to about 10 kilometers.This sampling interval is then used to look at the spatio-temporal dependencies, as illustrated in Figure 5b.Around a year temporal distance, one can see a clear correlation, which corresponds to the seasonal cycle of glacier velocity.From this variogram a rough scaling was estimated, and the anistropy was set towards a factor of four.In our case the pixel spacing is 300 meters and the time separation is 32 days.

Method performance
Two different temporal networks (combinations of time intervals) can be formulated in order to calculate a velocity estimate, as is described in section 3.1.The "open" configuration includes a greater number of velocity estimates from image pairs, but this has consequences.It results in a more complete dataset, with coherent velocity fields, but when short-term glacier dynamics occur, temporal resolution of the event may be smeared.
Our methodology thrives when several displacement estimates are present, otherwise testing is not possible.This prerequisite is especially apparent when a "closed" configuration is used, then the collection of displacements are reduced.In order White regions correspond to data without an estimate.
to show this dependency a slice is shown in Figure 6, in this time interval the Western side has several overlapping displacement estimates in space and time, while this is limited at the Eastern side.The exact distribution of the available data for this time interval is also highlighted in Figure A1, within the appendix.The lack of data can be seen at the western border of both both estimates.On the eastern side the glacier velocity structure is more clearly visible.Especially in the "open" configuration .
In Figure 7 more details of the two different configurations are shown.However, by including more imagery as with an open configuration, the velocity estimates are more complete, as can be seen along the outlets of Guyot Glacier.In the same subsection, not only the completeness increase but also the consistency, which is mostly appearant by coherent low velocities over stable ground.With more displacement vectors in the configuration, smaller-scale details, such as tributaries flowing into the large Kaskawulsh glacier, become more evident.
The spatio-temporal least-square estimates still have to some extent variation in direction and amplitude as well as outliers.
Therefore spatial-temporal smoothing is applied, in order to extract a better overview from the data, as is described in section 3.3.The results of this smoothing for the same time interval as in Figure 6 is shown in Figure 7.
In the smoothing procedure the surroundings of glaciers, which are stable-or slow moving terrain, are included.Consequently, high speed-ups such as on the surge bulge on the Steele glacier are dampened, as in this case it has a confined snout within a valley.They do not disappear, as the signal is strong and persistent over time, but damping does occur.An aspect of concern is the velocity bulge retreat of many outlet glaciers; their fronts with large velocities seem to retreat in the smoothed version, while this is not the case for the original least-squares estimate.This is an effect caused by surrounding zero-valued water bodies.Damping also occurs at turns such as at Hubbard glacier, the mask does reduce the effect of stable terrain, but has no specific glacial properties.The isotropy function included into the smoother might work for the direct neighborhood, but breaks down for fast moving outlets.For this smoother, weights are given in relation to neighborhood.However for glacier flow, the magnitude might be more similar in the direction of the flow lines, while in cross-flow direction the orientation might be more similar.This relation is not included in the smoother and therefor damping occurs.Nonetheless, the damping effect of the smoothing outweighs the advantage of having a clear image and overview over the mountain ranges.
Because the surrounding terrain, which has no movement, impacts the smoothed glacier velocity estimates in particular for surge and calving fronts (i.e. for strong spatial velocity gradients), the smoothing can be supported by a glacier mask.In our case, this mask is a rasterization of the Randolph Glacier Inventory (Pfeffer et al., 2014), with an additional dilation operation, to take potential advance or errors in the inventory into account.The difference in result for this masking procedure is shown in Figure 7, with some highlights.In general, the mask does compensate a little bit for the damping, but because the regions are mostly covered with ice its effect is small.

Validation of post-processing procedures
The voting used in our procedure is assessed through validation with an independent velocity estimate.Terrestrial measurements are limited in the study area, hence we use satellite imagery from RapidEye satellites over a similar timespan.Data from this constellation has a resolution of 5 meters and through processing in a pyramid fashion, a detailed flow field can be extracted.This velocity field functions as a baseline dataset to compare the GoLIVE and the synthesized data.Here we will look at a section of Klutlan Glacier, which flows from west to east, thus aligned with one of the map axis.The velocity of this glacier is, due to its surge, of significant magnitude, and therefor will have a wide spread in the voting space.
The two RapidEye images used over Klutlan Glacier were taken on 7 th of September, and on 7 th of October 2016.To retrieve the most complete displacement field of the glacier, we used a coarse-to-fine image matching scheme.The search window decreased stepwise (Kolaas, 2016) and the matching itself was done through Orientation Correlation (Heid and Kääb, 2012a).
At every step a local post-processing step (Westerweel and Scarano, 2005) was implemented, to filter outliers.The resulting displacement field over one axis (that is x, the general direction of flow) for this period is illustrated in the upperleft inset of measures are frequently used to produce and analyze multi-temporal datasets (Dehecq et al., 2015).
When looking at this time period for the GoLIVE data, a clear displacement field is shown, as both images (11th Sept., 13th Oct.) from Landsat 8 were cloud free.The pattern is in close agreement with the RapidEye version.When looking at the voted estimate a similar pattern is observable but more corrupted.In some respect the median estimate grasps the direction of flow, but over estimates the velocity, probably due to the surge that occurs.The spread might confirm this, as shown by the MAD, as this is considerable and will not help to justify which displacement is correct.Furthermore, the voted estimate is an estimate over a short interval, while the median estimate is calculated over the full stack.
To better assess these results, the distribution of both two-dimenstional displacement fields are illustrated in Figure 9a and 9b.Two groups of displacement regimes are clearly visible, a cluster showing little movement, or a flock of displacements with a dominant movement eastwards.The voted distribution has more spread, and outliers are present, but in general the mapping has the correct direction and magnitude.When the x-component of these displacements are compared against the RapidEye displacements.The median of this difference is 0.45 m/day for the voting and 0.27 m/day for the good GoLIVE pair.A similar trend can be seen in Figure 9c, which again shows the distributions are similar.Hence, the illustrated validation do show the voting scheme is able to grasp the general trend of the short term glacier flow through a large stack of corrupted velocity fields.While the voted estimate is worse than the clean GoLIVE estimate, we stress that the chosen GoLIVE dataset is one clean example within a large collection of partly corrupted displacement fields.Hence it is a step towards efficient information extraction, though the implemented voting has many potential aspects of improvements.

Validation over stable terrain
A second component for validation is an analysis of the stable ground, and the effect of the smoothing of the voted estimates.

5
The non-glaciated terrain are the locations stemming from a mask.A similar mask, also based on the Randolph glacier inventory, is used within the GoLIVE pipeline.Here, displacements over land and non-glaciated terrain are used to co-align the imagery Fahnestock et al. (2016), as geo-location errors might be present in the individual Landsat imagery.The fitting is done through a polynomial fit, still, in general these should be random, and its mean should be zero.The distribution of these stable terrain measurements, more than 65 million in total, are illustrated in figure 10.Similar to the visual inspection already illustrated in Figure 6b, the distributions also show a clear improvement.This is a welcome property as the voted estimates still seem to be noisy with significant outliers.

Glaciological observations
When looking at the spatio-temporal dataset some patterns that are observed by others also appear in our dataset.For example, the full extent of Bering glacier slows down, as highlighted by Burgess et al. (2012), however our time series cover a period where the full deceleration towards a quiescent state can be seen.This observation of a slow down can also be made for Donjek Glacier (Abe et al., 2016) and Logan Glacier (Abe and Furuya, 2015), see suplementary video.In the time period covered by our study some surges appear to initiate.For example, our dataset comprises a surge traveling along the main trunk of Klutlan Glacier, see also Figure B1&B2 in the appendix.Glacier (for location see Figure 7).The used GoLIVE data configuration is shown in Figure A1, and the data is from the smoothed dataset.
When looking at the surge occurring at Klutlan Glacier, the dataset does capture the evolution of its dynamics, as can be seen in Figure 12.The surge initiation seems to happen in the central trunk of the glacier, as the surge front progresses downwards (with steady bulk velocities around four meters per day).The surge also propagates upwards mainly into the most westwards basin.The eastward basin does also increase in speed, but to a lesser extent, while the middle basin of this glacier system does not seem to be affected significantly.Though, its extent is mostly pronounced in the lower part of the glacier, as the upward creep of velocity increase is limited and does not reach the headwalls of any basin.In Landsat imagery of late 2017, there is no indication of any heavily crevassed terrain in the upper parts of these basins, which supports the hypothesis of a partially developed surge.
When looking at the velocities over the flowline of the glacier, as in Figure 12a, both the extension downstream as well as the upstream progression of the surge can be seen.Most clearly, the surge front seems to propagate downwards with a steady velocity, but seems to slow down around the 50 km mark (see dashed line in Figure 12a), as shown by the break in slope.Here, the glacier widens but the surge does continue.This can suggest ice thickness is homogeneous here or ice thickness does not seem to play an important role in surge propagation.
At the end of the summer of 2016 the tributary just north of the 20 kilometer mark of Klutlan Glacier seems to increase in speed.This can be confirmed by tracing the extent of the looped moraines, as in Figure 12b.In the same imagery the medial moraines of the meeting point of all basins are mapped as well.Here, the moraine bands before and after the event align well in the junction, indicating a steady or similar contribution over the full period.Or an insignificant effect, as the surge has not been developing into very fast flow.In contrast, the lower part of this glacial trunk has moraine bands that do not align.
The surge behavior we observe for Klutlan Glacier is not a special mechanism and similar propagation behaviour can be observed at other glaciers within the mountain ranges.For Fisher Glacier, a similar increase in speed is observed within the main trunk that later propagates downstream as well as upstream.Similarly, this seems to be the case for Walsh Glacier, where a speed increase at the eastern trunk triggers a surge on the northern trunk leading to glacier-wide acceleration.On its way the fast flowing ice initiates surges in tributaries downflow, but the surge extent also creeps upslope and tributaries that were more up-glacier from the initial surge start to speed up.This is also true for Steele Glacier, that develops into a surge and Hodgson Glacier is entrained into this fast flow as well.
These analyses were best observed with the help of an animation (see supplement) but the identification was done through a simple visualization of the spread of flow speed, see Figure 13.Here the surging glaciers stand out, but also most of the tidewater glaciers, which have a highly dynamic nature at their fronts.Not only large glaciers can be identified but also the dynamics in smaller tributaries.For example, a tributary of the Chitina Glacier seems to have pushed itself into the main trunk 500 000 600 000 700 000 800 000 900 000 within a two-year time-period, see figure C1.

Discussions
Synthesized velocity estimates from our post-processing chain of GoLIVE image-pair velocity determinations are dependent on the number and distribution of measured displacements (see Appendix A).Surface features in the same season have similar appearances, thus good displacement fields can be estimates from imagery which are a year apart, as is typical working practice (Heid and Kääb, 2012b;Dehecq et al., 2015).Such annual displacements fields could be helpful when areas are cloud covered for long periods, such estimates can then function as gap fillers in the least squares estimation.Furthermore, the adjustment model assigns equal weighting to individual displacements if no other information is available.Hence some velocity changes might be missed or blurred in time.Such a drawback might be overcome with spatial constraints, such as an advection pattern imposed on the data, although this would increase the amount of post-processing.
Another limitation of our method concerns the glacier kinematics that are constrained by our model.In the current implementation the deviation (σ) is dependent on the time interval.From a measurement perspective this makes sense, but the model does thus not inherently account for speed change.Hence, for long time intervals the fuzzy function forces the deviation to become slim.This reduces the ability to get a correct match, especially when glacier-dynamical changes are occurring.It might therefore be worthwhile to explore the improvements occurring when a fixed deviation is set instead.In addition, the low score over glaciated terrain, might indicate the deviation of the displacement is set too tight.When this deviation is given higher bounds, the score increases and such parameter can then be used as a meaningful measure.
The smoothing parameter used is a single global parameter that assumes isotropy.In order to fulfill this property the spatiotemporal data has been scaled accordingly.However, when severe data gaps are present, the velocity dataset still seems to jump.This will improve when more data is available, for example by including Sentinel-2 data or incorporate across-track matching (Altena and Kääb, 2017a).An increase in votes will result in a better population of the vote space, as can be seen in Figure 6.In addition, the voting score, that is the consensus score in the Hough space, can be used for the initial weighting for the smoothing procedure (w in equation 2).This might reduce the amount of iterations used by the robust smoother.Nonetheless, much improvement can be made by the smoother, as it has a simple neighborhood function and has no knowledge of glacier specific properties.

Conclusions
In the past couple of years the increase in the number of high quality optical satellite systems have made it possible to extract detailed and frequent velocity fields over glaciers, ice caps and ice sheets.The GoLIVE dataset is a repository of such velocity fields derived from Landsat 8, and available at low latency for analysis by the community.Discovery and exploration of this resource can be complicated due to its vast and growing volume, and the complexity of spatio-temporal changes of glacier flow fields.Hence, in this study we introduce an efficient post-processing scheme to combine ice velocity data from different, but overlapping time-spans.The presented methodology is resistant to multiple outliers, as voting is used instead of testing.
However, since cloud cover can still hamper velocity estimation and spatial flow relations are not incorporated, the resulting synthesized time-series still have gaps or outliers.We use a data-driven spatio-temporal smoother to address this issue and enhance the visualization of real glacier flow changes.
Our synthesized time-series has a monthly (32 days) temporal interval and 300 meter spatial resolution.The time-series spans 2013 to 2018 and covers the Saint Elias Mountains and vicinity.Within this study area, we identify several surges in different glaciers at different times and their development over time can be observed.Such details can even be extracted for small tributary glaciers.Even velocities for the snow-covered upper glacier areas are in general estimated accurately.Thus our synthesized time-series can provide an overview of where and when interesting glacier dynamics are occurring.
This study is a demonstration of the capabilities of the new GoLIVE-type remote sensing products combined with an advanced data filtering and interpolation scheme.We demonstrate that our method can be implemented with ease for a large region, covering several mountain ranges.The derived smoothed time-series data contains many subtle additional changes that could be investigated.If this time-series is combined with digital elevation model (DEM) time-series (Wang and Kääb, 2015), it becomes possible to look at changes in ice flux in great detail.
The presented velocity time-series has a high temporal dimension, especially in respect to the sensor 16-day orbit repeat cycle.Though temporal or spatial data-gaps are still present (due to the short temporal interval, cloud cover or visual coherence loss) this might partly be addressed by enlarging the temporal resolution or through additional data, such as from Sentinel-2 (Altena and Kääb, 2017a).Fortunately, harmonization with other velocity datasets can be easily implemented, because our procedure uses only geometric information and are not dependent on sensor type.With our framework it is thus possible to make a consistent time-series composed of a patchwork of optical or SAR remote sensing products.Competing interests.All authors declare no competing interests.

Figure 1 .
Figure 1.Nominal Landsat 8 footprints used over the studied region.The purple text colors annotates the different satellite paths of LAND-SAT, while the black text indicates the relative overpass time in days in respect to path 63.

Figure 2 .
Figure 2. Graphical and matrix representation of a network.Here acquisition pairs within a network are illustrated and written down in an adjacency matrix (AG).The dark gray squares indicate acquisitions within a period to be estimated.The connecting colors symbolize an open (red) or closed (blue) selection of displacements to be used for the velocity estimation over this period (v).

Figure 3 .
Figure 3. On the left is a graphical representation of Landsat 8 acquisitions at different times (t) illustrated as nodes and matching solutions with displacements (d) are shown as edges, within a network.The subscript for the displacement denotes the time interval.The velocities (v) are estimated through the collection of all these displacements.On the right the mathematical representation or its parameter (Hough) space for voting of displacements over two time instances is shown, in this case the parameter space of only v0,1 & v1,2 are shown.
three dimensional Hough space.

Figure 4 .
Figure 4. Two modifications used in this study, which deviate from the standard Hough-transform as given by the toy example in figure 3. (a) the change from a clear ideal line, towards fuzzy weights.(b) the extension towards a higher dimension, in this case the lines transform to planes, that intersect with each other.

Figure 5 .
Figure5.Experimental variograms over a slice of the stack and over a subset of the spatio-temporal stack.The colorbar along the axis of figure5a, is used for the coloring of figure5b.The aspect ratio of figure5bis at the resolution of the produced data-cube (32 days by 300 meters).In order to make this isotropic, the vertical axis needs to extent, so the spread is similar in both directions.

Figure 6 .
Figure 6.Least square estimates of velocities with different network configurations, see Figure 3 for a toy example of the terminology.The study region spread over several UTM zones, hence the dataset is in Albers equal-area projection (ALB) with North American Datum 83.

Figure 9 .
Figure 9. Figure 9a & 9b show the distribution of velocities for a section of Klutlan glacier, their map view are shown in Figure 8.In figure 9c the same x-component data as in Figure 8 is shown, but now the distributions are shown.

Figure 10 .
Figure10.Distribution of the speed over stable terrain, for displacements extracted from the voting process, or after spatial temporal smoothing.The mask used is within the inset.

Figure 11 .
Figure 11.Snapshots of ice speeds at different time instances from a data compilation for the summer 2016 surge occurring on Klutlan

Figure 12 .
Figure 12.The speed over the central flowline of Klutlan Glacier.The markings of this flowline are shown in figure 11a.In Figure 12b the convergence of different basins of the Klutlan glacier is shown, data is from a RapidEye acquisition on the 5 th of September 2013 and at the 23 rd of September 2017.For comparison the 2013 image is overlain with the two morraine positions.

Figure 13 .
Figure 13.Spread of variation in flow speed (using the difference between the 20 th and 80 th percentile) over the observed period (2013-2018).Different dynamic glaciers are encircled, and the square indicates the tributary glacier shown in Figure C1.

Figure A1 .Figure B1 .
Figure A1.Node network of the velocity fields used in this study, a total of 2736 velocity fields.

Figure B2 .
Figure B2.The temporal data shown in Figure B1 is combined and the least squares consensus is plotted against the smoothed estimate.The slanted line in all figures corresponds to the 1:1.

Figure
Figure C1.A tributary of Chitina Glacier surged in the period 2015-2016.Images are both acquired by Landsat 8, its location is indicated by a square in figure 13.The location of the time-series in figure C1b, is indicated by a red dashed circle in figure C1a.