Observing traveling waves in glaciers with remote sensing: New flexible time series methods and application to Sermeq Kujalleq (Jakobshavn Isbræ), Greenland

The recent influx of remote sensing data provides new opportunities for quantifying spatiotemporal variations in glacier surface velocity and elevation fields. Here, we introduce a flexible time series reconstruction and decomposition technique for forming continuous, time-dependent surface velocity and elevation fields from discontinuous data and partitioning these time series into shortand long-term variations. The time series reconstruction consists of a sparsity-regularized least squares 5 regression for modeling time series as a linear combination of generic basis functions of multiple temporal scales, allowing us to capture complex variations in the data using simple functions. We apply this method to the multitemporal evolution of Sermeq Kujalleq (Jakobshavn Isbræ), Greenland. Using 555 ice velocity maps generated by the Greenland Ice Mapping Project and covering the period 2009 – 2019, we show that the amplification in seasonal velocity variations in 2012 – 2016 was coincident with a longer-term speedup initiating in 2012. Similarly, the reduction in post-2017 seasonal velocity variations was coincident 10 with a longer-term slowdown initiating around 2017. To understand how these perturbations propagate through the glacier, we introduce an approach for quantifying the spatially varying and frequency-dependent phase velocities and attenuation length scales of the resulting traveling waves. We hypothesize that these traveling waves are predominantly kinematic waves based on their long periods, coincident changes in surface velocity and elevation, and connection with variations in the terminus position. This ability to quantify wave propagation enables an entirely new framework for studying glacier dynamics using 15 remote sensing data.


Introduction
Until recently, observations of glacier and ice stream motion were limited to velocity snapshots measuring motion over distinct time periods, most commonly averaged over multiple years or annually repeating (Rignot et al., 2011;Gardner et al., 2018;Moon et al., 2012). While the increase in spatial coverage of velocity measurements facilitated by the increasing availability 20 of satellite-based remote sensing observations has allowed for ice sheet-wide analysis, the complexity of glacier dynamics requires observations at multiple temporal scales. Rapid responses in ice velocity to changes in external forces, such as ocean melt rate or calving frequency, may be superimposed on longer-term responses to variations in surface melt, ice geometry, and the final term λ m 1 is an 1 -norm term that encourages a sparse number of non-zero coefficients. The function in curly brackets in Eq. 2 is a convex cost function, which provides a solution that is guaranteed to be globally optimal. Implementation and the procedure for solving Eq. 2 form is detailed by Riel et al. (2014).
The coefficient λ in Eq. 2 is a penalty parameter controlling the strength of the sparsity-inducing regularization. Schemes for choosing values of λ based on the amount of data available and the desired smoothness of the solution are discussed by Riel et al. (2014). In practice, the 1 -norm regularization can be applied to a subset of m, which is assumed to be sparse, and depending on the reference functions in the dictionary that correspond to this subset, fewer non-zero coefficients may result in a smoother time series reconstruction. This regularization approach results in a compact representation for transient variations, which can aid in determining the dominant timescales and onset times captured in the data while potentially improving detection of signals with a lower signal-to-noise ratio (SNR).

100
In this study, we model transient signals as linear combinations of third-order integrated B-splines (B i -splines), which exhibit one-sided behavior of a particular timescale (Hetland et al., 2012;Riel et al., 2014). By allowing for combinations of B i -splines of different timescales and onset times, we can reconstruct a wide variety of transient signals. Strong seasonal variations in ice surface velocity and elevation, like those observed on Jakobshavn Isbrae, exhibit large changes in amplitude from year-to-year (Joughin et al., 2010(Joughin et al., , 2018. Therefore, we use a linear combination of third-order B-splines to reconstruct seasonal signals 105 with time-varying amplitudes (Riel et al., 2018). To encourage seasonal coherency of the B-spline coefficients, we construct C m in Eq. 2 such that the B-splines co-vary with other B-splines that share the same centroid time within any given year. The covariance strengths are constructed to decay exponentially in time. The flexibility in representing potentially complex temporal variations afforded by this approach avoids the severe limitations of using a single or small subset of sinusoidal variations (e.g., Minchew et al., 2017) and allows for a framework of transient and periodic variations that readily admit physical interpretation. 110 The interpretability of the resulting posterior modelm (which in this study primarily represents surface flow speeds) in terms of external drivers and intrinsic dynamics of the glacier is a marked advantage of our approach over time series approaches based on singular value decomposition (e.g., Samsonov, 2019). This advantage is amplified when using the sparse regularization techniques to constrain the timing, duration, and amplitude of transient events that are superimposed on periodic variations, as we describe in this study. Importantly, the framework described by Eq. 2 also enables quantification of the formal uncertainty 115 estimates in the inferred time series (i.e., posterior modelm).
Uncertainties for the estimated model coefficients can be formally quantified by combining observational uncertainties contained in the data covariance matrix C d with the dictionary G and prior covariance matrix C m (Tarantola, 2005;Bishop, 2006): 120 whereC m is the posterior model covariance matrix. Lower coefficient uncertainties can thus be obtained by a combination of reduced data noise and lower prior uncertainties on those coefficients. Similarly, the posterior covariance matrix of the reconstructed time series can be formally computed as (Tarantola, 2005;Bishop, 2006): effect on the covariances between different model coefficients and time epochs via the off-diagonal values in the matricesC m andC d , respectively. Properly quantifying these uncertainties and covariances is critical for any subsequent interpretation or analyses using the modeled time series.
After using Eq. 2 to estimatem for a given data time series d, the seasonal and transient signals can be reconstructed aŝ whered S andd T are the reconstructed seasonal and transient signals, respectively, S S represents the set of indices corresponding to the seasonal B-splines, S T represents the set of indices corresponding to the transient B i -splines, and G j is the vector corresponding to the j-th column of G. Correspondingly, we can use these reconstructed signals to detrend the data, e.g.
Examples of the practical implementation of these multitemporal methods are provided in the Results section.
In some cases it may be desirable to enforce spatial coherency when solving form to reduce the influence of data noise onm (Hetland et al., 2012;Riel et al., 2014;Minchew et al., 2015). Examples of such cases involve low-amplitude signals with spatial wavelengths longer than a single data pixel and, more generally, cases where the data have low SNR values. The 140 software made available with this study allows the user to enforce spatial coherency (cf. Riel et al., 2018), but for the velocity and elevation time-dependent data sets used in this study, we do not enforce spatial coherence between neighboring pixels because doing so can be computationally expensive and is unnecessary in our case. Instead, we solve form for each pixel independently and justify this decision based on the fact that the SNR is generally high over Jakobshavn Isbrae, reducing the need for spatial coherency in the inversion process.

Data
Data used in this study provide information on the time-dependent surface velocity and elevation fields. In both cases, the raw remote sensing data were processed to individual time-stamped fields and made publicly available through different projects and publications. In this section, we briefly describe these data sets.
observations of glacier velocity variations with increasing temporal resolution as more data from more sensors became available. Over select glaciers like Jakobshavn Isbrae, 11-day and monthly repeat TerraSAR-X/TanDEM-X (TSX) SAR pairs for both ascending and descending orbits are available and provide much higher temporal resolution than is available on many 155 other glaciers. The GIMP velocity maps are formed using speckle-tracking techniques on each SAR pair (Joughin, 2002).
Speckle tracking is generally more robust to large-scale displacements than standard interferometric methods and allows for uncertainty quantification through estimates of the statistics of correlation measurements. These formal errors are provided with the GIMP velocity fields.
In this work, we use 555 GIMP horizontal velocity fields generated from TSX data covering the period 2009 -2019. The 160 velocity fields are provided with 100-m grid spacing, and the short repeat times allow for mostly complete spatial coverage of Jakobshavn Isbrae due to the high number of coherent surface features facilitating high SNR for the speckle tracking techniques.
During the course of this work, we discovered systematic discrepancies in velocity measurements between SAR pairs collected from ascending and descending orbits at the same geographic location and approximately the same time. These discrepancies are generally limited to the shear margins of the glacier and cluster in areas much smaller than the glacier. To improve the 165 overall quality of the time series, we masked out areas where the velocities collected from near-coincident ascending and descending orbits differed by more than 250 m/yr.
To extend the velocity time series, we include monthly GIMP velocity maps derived from optical Landsat imagery for the years 2004-2009(Jeong et al., 2017. This longer time series allows us to compare long-term trends in surface velocity to long-term changes in glacier terminus position. The optical velocity maps are provided at 100-m grid spacing on a grid that is 170 not aligned with the TSX velocity maps. Thus, we resampled the Landsat-derived velocity fields to the same grid as the TSX velocity fields.

Surface Elevation Data
Ice surface elevation will rise and fall in response to variations in mass balance and ice flow. While time-dependent measurements of surface elevation are generally available at a lower temporal resolution than measurements of velocity, the increased 175 availability of optical imaging satellites and robustness of photogrammetry techniques has allowed for the generation of timedependent digital elevation model (DEM) strips at sub-annual epochs. Here, we use publicly available DEMs from ArcticDEM and the Oceans Melting Greenland (OMG) mission.
The ArcticDEM initiative automatically produces 2-meter resolution DEM strips over all land area north of 60 • latitude using stereo auto-correlation techniques (Porter et al., 2018). We use DEMs generated from optical data collected over Jakobshavn 180 from 2010 to 2017. The geographical extent of the strips and their associated acquisition times are irregular, so we interpolate the strips onto a uniform spatial grid using inverse distance weighting while using nearest neighbor interpolation in time to sample the strips onto a uniform temporal grid.
Because ArcticDEM data are available from 2010 to 2017, a timescale shorter than the GIMP surface velocity data, we augment the surface elevation time series with DEMs produced annually by the Oceans Melting Greenland (OMG) mission 185 in March 2017, 2018, and 2019. The OMG DEMs were constructed with Ka-Band single pass SAR interferometry from the Glacier and Ice Surface Topography Interferometer (GLISTIN-A) instrument, providing high-precision (< 50 cm) elevation maps over 10-12 km-wide swaths over various marine-terminating glaciers in Greenland (OMG Mission, 2016). These interferometric products were processed to an intrinsic resolution of 3 meters and then georegistered to a ground spacing of approximately 3 meters. For comparison with velocity data, both ArcticDEM and OMG DEMs were resampled to a grid 190 spacing of 100 meters.

Calving front positions
To elucidate the connection between observed changes in the terminus position and observed variations in ice surface velocity and elevation, we use calving front positions from a variety of sources. The main source are calving front positions automatically determined from TSX images acquired from 2009 to 2015 . We supplement these data with calving is a generalization of the often-mentioned process of retreat of the calving front into an overdeepened basin (e.g. Nick et al., 2009;Joughin et al., 2012) and has been reproduced in three-dimensional models incorporating detailed bed topography (e.g. Morlighem et al., 2016). We therefore generate a time series of the intersection of the calving front with the glacier centerline, referenced to the downstream position of the aforementioned shallow bed location. To reduce time series noise and facilitate comparison with the velocity time series, we again apply our linear least squares method to represent the front time series as 210 a combination of B-splines for short-term, seasonal signals and integrated B-splines for long-term transient signals. Since the absolute position of the front relative to the shallow bed location is important, we do not isolate the front seasonal (short-term) signal from the time series in the subsequent analysis.

Results
Hereafter, we focus on applying the time series analysis methods presented in Section 2 to analyze and decompose the observed 215 time-dependent velocity magnitude (speed) and surface elevation fields summarized in Section 3 into sub-annual (primarily seasonal) and multi-annual transient variations. Our focus is on quantifying the rates and distances over which stress perturbations of various frequencies propagate through Jakobshavn. We use only the time-dependent flow speeds because Jakobshavn flows along a deep, narrow channel in the underlying bed, which leads to temporally consistent velocity unit vectors throughout the observation period. Future work will be aimed at adding multidimensional capability to the time series analysis methods 220 we introduce in this study. Before proceeding to the analysis of the time-dependent velocity fields, we note for completeness that we observe magnitudes and spatial patterns of mean flow speed that are consistent with previous studies, with the highest speeds near the terminus decaying quasi-exponentially with upstream distance ( Figure 1A,B; Movie S1).

Seasonal Variations in Surface Velocity
The reconstructed velocity time series demonstrates the ability of our flexible method to smoothly interpolate the velocity 225 data in time in a manner that preserves the seasonal variations. In particular, the use of temporally coherent B-splines to model  Figure 1A,B represents a classical view of spatiotemporal variations in surface velocity, with a map of secular velocity ( Figure 1A) and time series of select points on the glacier ( Figure 1B).
235 Figure 1C shows, for the same points on the glacier, the seasonal variations, which are the total signal shown in Figure 1B less the inferred multi-year trends discussed in the next section. In Figure  km/year. Thus, there is a clear correlation between mean flow speed for a given year and the amplitude of seasonal variations, which we will explore in later sections.
In addition to attenuation, we are interested in constraining the rate of propagation of surface velocity variations. We quantify 255 these variations in terms of phase velocities by constraining the relative timing of peak velocity for different variation temporal frequencies.
As with amplitude, we present the absolute and relative timing of the velocity variations in multiple ways to help build a more intuitive framework for the reader. This presentation follows the same structure as the amplitude variations discussed in detail above, with a classical view shown in Figure 1, the space-time diagram in Figure 2, the mean over all seasons along a centerline transect in Figure 3B, and map view of the relative timing in Figure 4C (with formal uncertainties 260 in timing given in Figure 4D).
For seasonal variations in surface velocity, the time of peak velocity varies slightly from year-to-year, with the earliest peaks occurring around mid-August and the latest peaks occurring around mid-September. However, the spatial pattern of relative timing in the upstream direction from the terminus is broadly consistent even among years with large differences in mean velocities, which indicates a common mechanism for the seasonal cycle ( Figure 2A). As expected for marine-terminating 265 glaciers like Jakobshavn, the timing of the peak seasonal signal indicates that seasonal variations originate at the terminus and propagate upstream ( Figures 1C and 2A).
To investigate the spatial characteristics of the timing of peak velocity, we fit a simple temporal model using a sum of sinusoids to the velocity data from 2011 to 2018 at each pixel with the estimated long-term signals removed, d S (e.g., Figure   1C) such that where ω i is the angular frequency for the i-th sinusoid, C i and S i are the coefficients of the cosine and sine components, respectively, and C 0 is a constant offset. After estimating the values of C i and S i , the amplitude and phase (i.e., relative timing) and bootstrapped standard deviation for the transients (blue and orange). C) Centerline ice thickness (red) and bed depth (blue) using bed data from BedMachine V3 (Morlighem et al., 2017). For all plots, the gray shaded region represents the upstream region encompassed by the southern bend.
of each sinusoid can be recovered as While this model cannot accurately reproduce the amplitude changes or nonsinudoidal variations (e.g. Joughin et al., 2008Joughin et al., , 2014, the seasonal phase can be estimated robustly with seven years of data. Furthermore, we can compute the formal phase uncertainties following the procedure outlined in Minchew et al. (2017). The estimated seasonal phase is thus equivalent to the mean time of peak seasonal velocity for the 2011 -2018 period while the phase uncertainty is proportional to the formal 280 variance of the mean.
The seasonal phase map shows upstream transmission of velocity perturbations originating at the terminus ( Figures 3B and   4C). This propagation occurs rapidly in the first 5 km upstream of the terminus and then slows to a near-constant phase velocity of 398 ± 20 m/day (approximately 146 ± 7 km/yr), which is more than an order of magnitude faster than the mean flow speed near the terminus. Our estimates of the phase velocity within the first 5 km upstream of the terminus are limited by the temporal 285 resolution of the data, but by considering the uncertainties in timing, we estimate that the phase velocity in this region must be at least 500 m/day (182.5 km/yr), or approximately 18 times the local mean glacier flow speed. In the across-flow direction, about 8 km upstream from the 2017 calving front, the center of the glacier reaches its peak velocity earlier than the margins by about 15 days, indicating a nonlinear relationship between time-dependent lateral shear strain rates and centerline velocity, meaning that the effective width of the shear margins (defined here as the centerline velocity divided by the maximum shear 290 strain rate) must change over the seasonal cycle. Finally, we note that the phase uncertainty is generally lowest in the center of the glacier where amplitudes are higher and increases with distance upstream as the amplitudes decay ( Figure 4D), as expected from the formal uncertainties .

Multi-Year Variations in Surface Velocity
After isolating the long-term signals from the short-term seasonal signals, we can observe clear variations in multi-annual 295 amplitudes at different points along the glacier ( Figures 1D and 2B). The temporal density of the velocity time series allows us to quantify spatial variations in the amplitude and timing of the positive and negative multi-annual trends, much like we did with seasonal velocity variations in the previous sections. We observe two events in the data: a speedup that begins in 2012 and a slowdown that begins in late 2016 near the terminus, which we refer to as the 2017 slowdown. We present the results for multi-annual variations using the same general structure as for the seasonal, with the classical view shown in Figure 1, the 300 space-time diagram in Figure 2, along a centerline transect in Figure 3, and map view of the amplitudes and phase values in   (1σ). Seasonal amplitudes are measured as the difference between the summer high and winter low velocities in the short-term time series as shown in Figure 1C. The highest amplitudes occur at the terminus and decay exponentially upstream. amplitudes for the velocity variations become too low to reliably estimate the timing of arrival of the transient signals ( Figure   5D). The phase velocity between 8 and 20 km upstream is approximately 231 ± 16 m/day (84 ± 6 km/yr), which is a little more than half of the phase velocity for the seasonal signal and roughly seven times the mean flow speed near the terminus. The upstream phase velocity also has features that suggest a sensitivity to ice thickness and bed topography ( Figure 3C), but more work is needed to establish concrete connections. As with the observed seasonal variations, the phase velocity in the first 5 km 315 upstream of the terminus is at the limit of the temporal resolution of the data with a lower bound on the phase velocity of at least 500 m/day (182.5 km/yr), or approximately 18 times the mean glacier flow speed in this region. While the slowdown in wave propagation in the southern bend is coincident with a local high in the bed topography, more work is needed to evaluate whether the topographic effect is the dominant control on wave propagation. The apparent slowdown may also be an artifact of numerical errors caused by tracking of peak acceleration/deceleration rather than a multi-year average of sinusoidal phase 320 13 https://doi.org/10.5194/tc-2020-193 Preprint. Discussion started: 7 August 2020 c Author(s) 2020. CC BY 4.0 License. as for the seasonal signal. Nevertheless, the consistency in the spatial distribution of peak timing and amplitude reinforces the notion that a common physical mechanism is responsible for multi-year and seasonal velocity variations. In addition to the spatial distribution of phase delay and amplitude being consistent for the two events, these transient phase delays show strong similarities with the seasonal phase delay. However, the transient amplitudes show farther upstream propagation than the seasonal amplitudes.

Multi-Year Variations in Surface Elevation
Ice surface elevation varies in response to changes in snow accumulation and melt (the sum of which constitutes the surface mass balance; SMB), firn compaction (Herron and Langway, 1980;Huss, 2013;Meyer et al., 2019a), and dynamic thinning 325 (thickening) in response to increases (decreases) in the flux divergence of the ice. The interplay between observed elevation and velocity changes at different temporal and spatial scales can thus yield insight into the mechanisms driving longer-term elevation and velocity changes.

Calving Front Forcing
The position of the calving front has been reported to be the dominant factor influencing ice dynamics for Jakobshavn Isbrae over observable, particularly seasonal, timescales (Nick et al., 2009;Joughin et al., 2012;Bondzio et al., 2017). The position of 345 the calving front is heavily influenced by bed topography (Xu et al., 2013;Morlighem et al., 2016), local ocean temperatures that influence subaqueous melting of the terminus (Holland et al., 2008;Khazendar et al., 2019), and changes in melange rigidity (Joughin et al., 2008;Cassotto et al., 2015;Robel, 2017;Xie et al., 2019;Joughin et al., 2020). Changes in the calving front driven by iceberg calving can reduce back stresses upstream of the new calving front, allowing ice near the terminus to accelerate and subsequently thin. Local thinning of the ice steepens the surface slope, thereby increasing gravitational driving 350 stress and causing further acceleration and thinning . For glaciers where the calving front is located on a retrograde bed, thinning-induced retreat of the terminus corresponds to increasing terminus ice thickness, which further increases the driving stress (Joughin et al., 2020). The redistribution of mass during dynamic thinning propagates upstream as traveling waves, with phase velocities and attenuation length scales that we constrain using the methodology presented in this study. The thinning of the ice surface throughout the glacier is hypothesized to cause a further, indirect enhancement of 355 velocities by causing a reduction in overburden pressure (weight of the ice column per unit area), which influences the response of the glacier to processes such as basal lubrication by drainage of surface meltwater in the summer melt season. While the role of basal lubrication has been shown to affect the seasonal cycle on certain glaciers (e.g. Howat et al., 2010;Bevan et al., 2015), the magnitude and timescale of the influence of basal lubrication on the velocities in Jakobshavn remains unresolved.

Discussion
Decomposition of the time-dependent velocity and surface elevation fields into distinct temporal scales reveals a repeating pattern on Jakobshavn where velocity and surface elevation variations originate at the terminus in response to changes in calving front position. The coincidence between speedup and slowdown of the glacier with thinning and thickening, respectively, suggests a dynamic origin to the physical mechanism generating these variations. In this section, we detail the observed wave 395 phenomena, introduce the proposition that the observed traveling waves are kinematic waves, and discuss possible paths for future development of observational methods that will enable progression toward robust and efficient techniques for fusing remote sensing data from multiple sources and using in situ observations as prior information to constrain the inversions.

Wave phenomena
Our results indicate that velocity variations initiating at the terminus of Jakobshavn propagate upstream as traveling waves with 400 frequency-dependent propagation speeds (phase velocities) and attenuation length scales. To our knowledge, ours are the first which we find to be 231 ± 16 m/day (84 ± 6 km/yr). Thus, whether associated with retreat (glacier speeds up) or advance (glacier slows down) of the terminus position, multi-annual signals propagate upstream roughly seven times faster than the mean flow speed near the terminus of Jakobshavn, or more than an order of magnitude faster than the local mean flow speeds.
Taking both of the observed multi-annual signals to have periods of 3 years, we estimate the wavelength of the traveling waves to be 2πc p /ω ≈ 252 ± 18 km, where c p is the phase velocity and ω the angular frequency. At ∼50 times the width, ∼200 times the thickness of the glacier near the terminus, and ∼5 times the length of the glacier, the wavelength of observed multi-annual waves is much longer than any spatial dimension of the glacier. For seasonal frequencies, we observe phase velocities in the upstream 8-20 km section of 398 ± 20 m/day (146 ± 7 km/yr), which is more than an order of magnitude faster than the glacier flow speed near the terminus. Taking the period of seasonal variations to be 1 year, we estimate the respective wavelength to be ≈ 146 ± 7 km. At ∼30 times the width and ∼100 times the thickness of the glacier near the terminus, the wavelength of 445 seasonal variations is much longer than typical characteristic length scales for glacier flow but shorter than the multi-annualperiod waves by a factor of approximately 1/ √ 3, as discussed below. The attenuation length scales we report for both seasonal and multi-annual-period traveling waves are more than an order of magnitude shorter than the wavelengths, consistent with previous studies of kinematic waves that found diffusive behavior with attenuation timescales shorter than their periods (Lick, 1970;Jóhannesson, 1992;Gudmundsson, 2003).

450
The wavelengths of the observed waves are perhaps clearer when considered in terms of the phase velocities. Using the same periods, we see that the ratio of phase velocities is approximately equal to the square root of the ratio of the frequencies, such that c multi−annual in agreement with the estimates of a wavelengths just discussed. This estimate places useful constraints on the form of the 455 dispersion relation, which we will explore in future work. We note that the observed attenuation length scales appear to deviate somewhat from the square-root-frequency relation as in Eq. 10, with the attenuation length scale for the seasonal signal being approximately half that for the multi-annual signal (Fig. 3). As with the phase velocities, this relationship between attenuation and frequency provides useful constraints on the dispersion relation.
Given the estimates of phase velocity at two different frequencies, it is desirable to estimate the group velocity, defined 460 as the rate at which the envelope of a wave packet travels, and thus the rate at which the energy of the wave packet travels.
Mathematically, group velocity is defined as c g = ∂ω/∂k r , where k r = ω/c p is the real component of the angular wavenumber.
Our estimates of group velocity are crude given that we only observe waves with two frequencies and the period of the multi-annual signal is not well-defined because we use an integrated B-spline as a smooth step function to fit the transient.
Nonetheless, assuming, as before, that seasonal signals have a period of one year and the observed multi-annual signals both 465 have periods of three years, we estimate a group velocity of c g ≈ 230 km/yr, which is faster than the respective phase velocities.
Therefore, in this relation between phase and group velocities, the waves we observe are analogous to capillary waves, though the analogy goes no further as the physical mechanisms governing capillary waves are not all applicable to the low-Reynolds number regime of glacier flow.
we correct for the differences in mean flow speed between the two glaciers. Traveling wave speeds on Mer de Glace following numerous perturbations in surface mass balance were synthesized by Lliboutry and Reynaud (1981) to be in the range 450-725 m/yr, or about 4-6 times the mean flow speed of the glacier. The traveling waves we observe on Jakobshavn are moving two orders of magnitude faster upstream than traveling waves on Mer de Glace. The ratio of phase velocity to local mean flow speed is roughly a factor of two higher on Jakobshavn than Mer de Glace. There are marked differences in these two glaciers 475 and the perturbing forces are quite different -retreat of the terminus on Jakobshavn and perturbation in surface mass balance on Mer de Glace -but the qualitative differences in observed wave propagation speeds are noteworthy.
When comparing the results presented here with observed wave propagation on other glaciers, we note dramatic differences between the phase velocities and attenuation length scales on Jakobshavn and wave propagation related to much higher frequency (14.76-day period) variations in ice-surface velocity that we previously reported for Rutford Ice Stream, West Antarc-480 tica . On Rutford -with a mean flow speed near the grounding line of approximately 375 m/year -downstream variations in buttressing stresses over the ice shelf, grounding line position, and/or water pressure at the bed (Gudmundsson, 2006(Gudmundsson, , 2007Rosier et al., 2014Rosier et al., , 2015Minchew et al., 2017;Robel et al., 2017;Rosier and Gudmundsson, 2020) drive variations in the flow speed of Rutford that propagate with a phase velocity of approximately 24 km/day for the first 40 km upstream of the grounding line, and then at a faster rate of 34.3 km/day further upstream. Thus, observed waves on 485 Rutford driven by ocean tides propagate two orders of magnitude faster than those we observe on Jakobshavn. The attenuation length scales are also markedly different: approximately 45 km on Rutford versus 7 km (seasonal) and 14 km (multi-annual) on Jakobshavn. The marked differences in forcing frequencies and propagation speeds and distances suggest fundamental differences in wave types and mechanisms.

Kinematic wave propagation 490
We hypothesize that the traveling waves we observe can be classified as kinematic waves because of their long periods (months to years), strong correlation with calving front position, and marked dynamical thinning of the glacier corresponding to increases in velocity. Kinematic waves represent a special case in a broader spectrum of wave behavior that includes various forms of dynamical waves. As the name suggests, dynamical waves are driven primarily by pressure and stress gradients and arise in a variety of contexts, such as seismic waves, flexural waves (Lipovsky, 2018), shallow-water waves, so-called "sea-495 sonal waves" on glaciers (Fowler, 1982;Hewitt and Fowler, 2008), and the response of outlet glaciers to ocean tides (Rosier et al., 2015;Rosier and Gudmundsson, 2016;Minchew et al., 2017). On the other hand, kinematic waves arise primarily from the redistribution of mass, with propagation characteristics dominated by mass conservation rather than momentum balance.
Kinematic waves have been used to model phenomena as varied as glacier surges, river floods, and traffic flow (Lighthill and Whitham, 1955;Nye, 1958Nye, , 1960. All waves on glaciers will be driven by some combination of dynamic (momentum bal-500 ance) and kinematic (mass balance) processes, and the dominance of one driving mechanism over another informs the broad characteristics of the full dispersion relation for any given glacier at any given time.
Our hypothesis that the waves we observe are kinematic waves is consistent with previous studies that explore various aspects of kinematic wave propagation and the response of tidewater glaciers to changes in terminus position, from observations (Howat et al., 2005;Joughin et al., 2008;Felikson et al., 2017) to flow models (Thomas, 2004;Pfeffer, 2007;Nick et al., 2009;Joughin 505 et al., 2012;Williams et al., 2012). Often, these arguments are based on the logical conclusion of progressive draw-down of the glacier surface. In essence, this process plays out such that changes in the terminus position through calving locally alter longitudinal (normal) stresses, allowing the glacier to locally accelerate. The local acceleration alters longitudinal stresses and increases flux divergence, causing the glacier to thin. Localized thinning steepens the surface slope, thereby increasing gravitational driving stress and generating upstream acceleration . The process repeats as the surface 510 is progressively drawn down and a wave propagates upstream. Similar mechanisms operate in reverse for re-advance of the terminus position. This explanation is broadly consistent with our observations showing that changes in surface elevation are coincident with changes in surface velocity (Fig. 6). We reserve for future work a detailed exploration of wave propagation in glaciers, including tests of the kinematic wave explanation of the observations we present in this study.
However, we note two important caveats to the applicability of the kinematic wave description to our observations. The first 515 is that the wave speeds we observe are an order of magnitude faster than the local mean flow speeds, making them much faster than kinematic wave speeds proposed previously (Nye, 1960;Lliboutry and Reynaud, 1981;Weertman and Birchfield, 1983;van de Wal and Oerlemans, 1995). This discrepancy likely arises from a combination of the assumptions made in previous estimates and the fact that our observed phase velocities are not corrected for diffusivity, and so are not direct measurements of the kinematic wave speed (Lliboutry and Reynaud, 1981). As an example of the assumptions made in earlier estimates of 520 kinematic wave speeds, Nye (1960) assumed a shallow-ice approximation for the flow regime (wherein gravitational driving stress is balanced locally at the bed) and a power-law relation between slip-rate and drag at the ice-bed interface (i.e., the sliding law) with values for the exponent defined by Weertman (1957) for regelation and viscous flow of ice around roughness features in a monochromatic bed. Such a sliding law does not account for the full range of possible sliding mechanisms (Schoof, 2005;Joughin et al., 2019;Zoet and Iverson, 2020;Minchew and Joughin, 2020), meaning that estimates from Nye (1960), and 525 others that follow the same basic formulation, do not capture the range of possible values for kinematic wave speeds. More work is needed to reconcile our observations with models of traveling wave propagation.
The second caveat to our conclusion that the waves we observe are kinematic waves is that our observations indicate that the observed waves are dispersive. Recall that we define dispersion in the traditional sense to mean "frequency dispersion" -the frequency dependence of wave speed -whereas some authors, notably Lighthill and Whitham (1955), use the term "amplitude 530 dispersion" to describe the amplitude dependence of wave speed. The latter is a general characteristic of kinematic waves that can cause waveforms to change shape as faster waves overtake slower waves (Fowler, 1982(Fowler, , 2011. However, the amplitude dependence of wave speed does not account for our observations for three reasons. First, the phase velocity of waves with seasonal periods are relatively consistent despite large differences in the amplitude of seasonal variations (Fig. 2). Second, the multi-annual variations in glacier flow speeds have comparable amplitudes to the seasonal variations and yet propagate slower 535 by a factor of √ 3, nearly equal to the square root of the ratio of frequencies (Eq. 10 and Fig. 3). Third, the observed phase velocities for the 2012 speedup and 2017 slowdown propagate at nearly the same speed (Fig. 3) despite different amplitudes of variations in surface flow speed (Figs. 1 and 2) and elevation (Fig. 6). Thus, we conclude that differences in amplitude cannot explain the disparities in phase velocities from waves with seasonal and multi-annual periods. Our observations, therefore, provide evidence of (frequency) dispersion in traveling waves.

540
Evidence for dispersion is noteworthy in part because kinematic waves on glaciers are governed by the continuity equation, which can be expressed in the form of the first-order wave equation; kinematic waves are, therefore, non-dispersive in the classical theory (Lighthill and Whitham, 1955;Nye, 1960). More recent work by, for example, Pfeffer (2007)  stresses would entail inclusion of dynamical processes operating with greater or equal importance to the kinematic wave processes. As mentioned previously, kinematic waves are special cases in a broader spectrum of wave phenomena that must include various types of dynamical waves, and we should expect there to be a range of frequencies in which both mass and 550 momentum balances play non-negligible roles in wave propagation.
The rheology of ice is another candidate for dispersion. We may intuitively expect waves with periods comparable to the viscoelastic relaxation time (ratio of dynamic viscosity to shear modulus) to be dispersive as the combination of viscous and elastic responses mean that waves of different frequencies will travel through ice with with different effective rheologies. In this case, the restoring force provided by elasticity is proportional to strain, which will vary with wavelength. However, the 555 periods of waves we observe are much longer than the viscoelastic relaxation time -which is typically of order hours to weeks (Gudmundsson, 2007) -meaning that ice is essentially a purely viscous fluid over the timescales relevant for our observations, and therefore elasticity is an unlikely source of dispersion. The non-Newtonian viscosity of ice is also unlikely to generate dispersion in kinematic waves due to the long periods of the wave relative to the viscoelastic relaxation time, which is supported by previous studies of kinematic waves (Nye, 1960;van der Veen, 2001;Pfeffer, 2007;Felikson et al., 2017).

560
Other obvious physical processes that could give rise to dispersion in traveling waves with seasonal and multi-annual periods are related to water pressure at the bed. In the interest of brevity, we divide the relevant processes into two broad categories: subglacial hydrology and till porosity. For subglacial hydrology, the well-known dependence of the hydrological pressure gradient on the surface slope of the glacier may be an important consideration (Flowers, 2015). After all, kinematic waves propagate upstream due to a progressive draw down of the surface, which entails a transient steepening of the surface slope.

565
The response of the subglacial hydrological system to the transient changes in glacier geometry are beyond the scope of this study and are worthy of further consideration in future work. As another example, dilation and compaction of subglacial till has recently received attention as a possible mechanism for triggering surges in glaciers with deformable beds (Minchew and Meyer, 2020) and has long been considered a potentially important mechanism in the centennial timescale dynamics of ice streams (Tulaczyk et al., 2000;Robel et al., 2016;Meyer et al., 2019b). In this mechanism, till dilation is influenced by the 570 overburden pressure at the bed (i.e., the weight of ice per unit area) and the change in slip rate at the bed. As the surface velocity increases, it is likely that a corresponding increase in slip rate manifests at the bed and can be expected to cause the deforming till layer to dilate. For realistic values of hydraulic permeability of the till, dilation will cause a temporary decrease in pore water pressure. Minchew and Meyer (2020) showed that this processes can (in an idealized model) lead to glacier surges by delaying the evolution of the till to a new steady-state. In a related sense, dilation and the subsequent evolution of 575 pore water pressure provide a possible mechanism for altering the mechanical properties of the glacier bed in such a way that might generate dispersion in waves with seasonal and multi-annual periods.

Future work in remote sensing
The increasing availability of surface velocity and elevation fields sampled at monthly-to-sub-monthly timescales will continue to provide opportunities to study the rapid evolution of fast-flowing glaciers to various environmental forcings. The operational 580 capabilities of several working groups that produce velocity fields over the Greenland and Antarctic Ice Sheets will consistently improve as new data are made available and techniques for generating velocity estimates are refined. In particular, the upcoming NASA-ISRO Synthetic Aperture Radar (NISAR) mission will generate unprecedented volumes of data that are useful for quantifying surface change for a number of scientific applications, including glacier dynamics (NISAR, 2018). The wide imaging swath (~240 km) coupled with a 12-day repeat cycle and global coverage will allow for systematic observations of 585 high-resolution velocity variations over interconnected glacier networks and coupled ice stream and ice shelf systems. Such observations will facilitate quantification of the spatiotemporal responses of glaciers and ice streams to any changes to the stress state, such as changes to the terminus position, loss of ice-shelf buttressing, changes in frictional properties of the bed, evolution of the subglacial hydrology. These processes will likely result in wave phenomena similar to those observed at Jabovshavn Isbrae (this study) and Rutford Ice Stream  and would be well-observed with platforms like 590 NISAR. Furthermore, the quantification of phase velocities and attenuation length scales at multiple forcing frequencies would provide valuable constraints on a general theory for wave propagation for fast-flowing glaciers because the characteristics of wave propagation are intrinsic properties of any given glacier system, which includes the boundary conditions.
The temporal resolution of surface DEMs is currently a limiting factor in quantifying sub-annual dynamical thinning. In this study, we noted that the thinning signal in the ice adjacent to the fast-flowing regions may be due to oversmoothing of the time 595 series due to the limited temporal resolution of the ArcticDEM dataset. Therefore, elevation or altimetry datasets that have increased temporal sampling, such as IceSat-2, may help isolate the shorter-term dynamic signals from any longer-term SMBbased variations. In particular, future analysis would benefit from the 91-day repeat time of IceSat-2 for capturing seasonal elevation variations for direct comparison and synthesis with the velocity seasonal variations. This type of analysis could lead to a full three-dimensional velocity time series, which has the potential to improve quantification of strain and stress fields, 600 constraints on ice rheology, and assimilation of velocity data into state-of-the-art ice flow models. For glaciers and ice streams where a persistent ice shelf or tongue exists, tidal forcing may become an important stress perturbation, in which case accurate reconstruction of vertical displacements would be necessary in order to constrain the dominant tidal constituents of motion .
The flexible time series framework described here introduces the potential for using in situ observations as prior information 605 (encoded in the prior model covariance matrix C m ) in forming time-dependent surface velocity fields. One example of this synergy between remote sensing and in situ data is the use of GPS/GNSS observations to constrain the form of the temporal basis functions, as we did in Minchew et al. (2017). Another, less obvious, example is the potential for employing catalogs of calving events gleaned from seismic observations Nettles, 2017, 2019;Olinger et al., 2019) to constrain the timing and duration of transient accelerations in ice flow. Such constraints on the temporal evolution of the fields observed 610 from remote sensing observations should afford novel opportunities to constrain phenomena such as the localization of strain rates (and, thereby, stresses) associated with fracture and calving. We expect the usefulness of the flexible methods we present here to grow as more remote sensing and in situ data become available.

Conclusions
We have presented a framework for forming continuous time-dependent surface velocity and elevation fields from publicly 615 available surface velocity and elevation data. This framework is based on a sparsity-regularized linear regression method that reconstructs time series as a linear combination of relevant basis functions. The flexibility and expressive power of the basis function representation allows for accurate reconstruction of time series in the presence of noisy and missing data while also allowing for a natural decomposition of the total signal into signals of multiple temporal scales. Over Jakobshavn Isbrae, this decomposition permitted a detailed investigation into the spatiotemporal characteristics of the evolving seasonal cycle of ice 620 speedup and slowdown in response to terminus variations. Analogously, longer-term changes in velocity were isolated and shown to also be primarily driven by longer-term terminus variations. This type of analysis is directly applicable to many outlet glaciers in Greenland, Antarctica, and other areas where multitemporal remote sensing data is available and could improve our understanding of the dynamic response of glaciers to various geographic and environmental forcings.
We demonstrated that the time series reconstruction permitted the quantification of traveling wave propagation resulting 625 from terminus forcing functions at different temporal frequencies. These results build upon an important new area of research that aims to achieve a mechanistic understanding of glacier flow from time-dependent velocity data. To our knowledge, our results are first to show from observations that waves on glaciers with seasonal to multi-annual periods are dispersive, with a ratio of observed phase velocities approximately equal to the square root of the ratio of frequencies. We hypothesize that the observed waves can be classified as kinematic waves based on their long periods (much longer than the viscoelastic relaxation 630 time), correlation with changes in the terminus position, and coincident variations in surface velocity and elevation. These observations of traveling waves are only possible due to the strong velocity response to changes in terminus position, as well as our ability to isolate short-and long-term signals in the velocity data. Looking forward, we aim to assimilate other velocity sources for Jakobshavn Isbrae (e.g., optical or Sentinel-1), as well as other elevation and altimetry data sets to improve temporal sampling and to obtain full 3D surface velocity time series. The resultant dataset will likely lead to a marked improvement in 635 incorporating velocity data in ice flow models for simulation and inversion of mechanical properties.
Code availability. The time series analysis and decomposition tools are available at https://github.com/bryanvriel/iceutils.