Detecting seasonal ice dynamics in satellite images

Fully understanding how glaciers respond to environmental change will require new methods to help us identify the onset of ice acceleration events and observe how dynamic signals propagate within glaciers. In particular, observations of ice dynamics on seasonal timescales may offer insights into how a glacier interacts with various forcing mechanisms throughout the year. The task of generating continuous ice velocity time series that resolve seasonal variability is made difficult by a spotty satellite record that contains no optical observations during dark, polar winters. Furthermore, velocities obtained by 5 feature tracking are marked by high noise when image pairs are separated by short time intervals, and contain no direct insights into variability that occurs between images separated by long time intervals. In this paper, we describe a method of analyzing opticalor SAR-derived feature-tracked velocities to characterize the magnitude and timing of seasonal ice dynamic variability. Our method is agnostic to data gaps and is able to recover decadal average winter velocities regardless of the availability of direct observations during winter. Using characteristic image acquisition times and error distributions from Antarctic image 10 pairs in the ITS_LIVE dataset, we generate synthetic ice velocity time series, then apply our method to recover imposed magnitudes of seasonal variability within ±1.4 m yr−1. We then validate the techniques by comparing our results to GPS data collected on Russell Glacier in Greenland. The methods presented here may be applied to better understand how ice dynamic signals propagate on seasonal timescales, and what mechanisms control the flow of the world’s ice.


15
Earth-observing satellites have been in orbit for over half a century, but it was only in 2011 that a sufficient quantity of data had been collected to complete the first pan-Antarctic map of ice velocity (Rignot et al., 2011). Since then, new satellites have led to follow-on mappings that identified regions of changing ice flow (Gardner et al., 2018), and today data are being collected at such a rate that velocity mosaics can be generated every year for nearly all of the world's land ice (Joughin, 2017;Mouginot et al., 2017;Gardner et al., 2019). So for a field of study that was born in an age of in situ stake networks and dead reckoning, 20 the data revolution of the past decade has completely upended glaciology. Where once our challenge was to squeeze as much information as possible from a few sparse field measurements, our biggest challenge now lies in processing massive, often unwieldy datasets, and finding all the meaningful signals that lie hidden in this new abundance of data.
One of the most direct and insightful ways to understand how ice moves, what controls its flow, and how it responds to changes in its environment, is to observe dynamic variability under a wide range of periodic forcings. Long-term trends and 25 interannual variability (e.g., Moon et al., 2012;Christianson et al., 2016;Greene et al., 2017a;Dehecq et al., 2019) can provide a sense of how glaciers respond to sustained forcing, but on much shorter timescales we are able to see how velocity variations initiate and what physical processes control a glacier's movement. Several targeted studies have shown that glaciers can exhibit observable dynamic responses to ocean tides, and that tidal signals can propagate well inland of the grounding line on daily to fortnightly timescales (e.g., Anandakrishnan et al., 2003;Walter et al., 2012;Rosier and Gudmundsson, 2016;Minchew 30 et al., 2017;Robel et al., 2017), but in many glaciers around the world, a significant gap exists in our understanding of how ice dynamic changes develop between tidal and interannual timescales.
It is certainly understood that many mountain glaciers speed up and slow down throughout the year (Burgess et al., 2013;Armstrong et al., 2017;Yasuda and Furuya, 2015;Kraaijenbrink et al., 2016), and that some of Greenland's glaciers respond to seasonal cycles of subglacial hydrology or calving dynamics (Joughin et al., 2008;Howat et al., 2010;Bartholomew et al., 35 2010; Sole et al., 2013;Moon et al., 2015;King et al., 2018). Seasonal variability has even been reported in a few studies of Antarctic glaciers (Nakamura et al., 2010;Zhou et al., 2014;Greene et al., 2018); but to date, no global-scale mapping of seasonal dynamics of the world's ice has been completed, due in part to the technical challenge of working with optical data in polar regions, where the surface is not touched by sunlight for months-long periods each winter. A comprehensive mapping of the world's seasonal ice dynamics would permit direct inter-comparison of seasonal evolution in regions with different driving 40 processes; provide a basis for analysis of long-term changes in seasonal behavior; and supply models with a zeroth-order understanding of global ice dynamics.
In this paper, we describe a precise method of measuring the average seasonal cycle of ice flow dynamics, with the ultimate goal that our method may be used to map the typical magnitude and timing of the seasonal glacier dynamics worldwide.
Our study is primarily focused on Antarctica, where seasonal variability is poorly understood, and where data limitations 45 currently present the greatest challenges to making such measurements. We test the sensitivity of our method on several thousand synthetic ice velocity time series, then validate it by applying the method to satellite data covering Russell Glacier in Greenland, where we compare our results to GPS observations that show persistent seasonality.

Feature-tracked velocity data
The method we present applies to ice velocity datasets such as GoLIVE  or ITS_LIVE (Gardner et al.,50 2019), which have been derived by feature tracking techniques applied to satellite image pairs. For a detailed review of the principles of feature tracking we refer readers to Scambos et al. (1992) or Fahnestock et al. (2016), but for the work presented here it is essential to know only that feature-tracked velocities are measured as the integrated surface displacement that occurs between the acquisition times of two satellite images of a given location. That is, each measurement represents an average velocity between image acquisition dates, and the required passage of time between images precludes direct measurement 55 of instantaneous velocity at any given time. As a result, any high-frequency variability that occurs between images is not represented, and seasonality may appear missing or deprecated in velocity measurements obtained by feature tracking over long timespans (Fig. 1). Nonetheless, by fitting a cyclic function to the time series of displacements rather than average velocities, we show that it is possible to accurately recover the magnitude and phase of seasonal velocity variability. The upper panel shows an ice velocity time series in blue, which integrates to form the cumulative displacement time series shown in the lower panel. We use satellite images to measure ice velocity as the cumulative displacement of crevasses and other glacier features that occurs between acquisition times of any two satellite images. Here, four images taken at times t1 through t4 provide six unique combinations of image pairs that yield the measurements depicted as horizontal gray lines connecting the acquisition times of each image pair. Vertical gray lines show measurement uncertainty and a black dot is placed at the center times of each image pair for visual clarity. A dashed sinusoid is fit to velocity measurements at the center times of each image pair to highlight the inadequacy of fitting directly to velocity data for determination of seasonal amplitude or phase. This paper describes an alternate, exact approach, wherein sinusoids are fit to accumulated displacements, which are then converted to velocity to produce the light red velocity sinusoid shown in the upper panel.
This paper focuses primarily on the Landsat image pairs that populate the ITS_LIVE dataset, in part because their record 60 extends back as far as 1985 in some locations. The long Landsat record may help ascertain a climatological seasonal cycle of ice dynamics at any given location; however, we face the limitation that optical satellites like Landsat do not collect data throughout the dark winters in polar regions, where land ice is most prevalent (see Fig. 2). Despite the lack of direct observation in winter, we will show that the magnitude and timing of seasonal variability can be accurately retrieved from Landsat data, regardless of when in the year the maximum velocity occurs.

65
For any given 240×240 m pixel in Antarctica, the ITS_LIVE dataset may contain from a few dozen to more than 10,000 velocity measurements (Fig. 3), which are taken as the ice displacement observed between two satellite images collected on different days. Each satellite image may serve as the first or second image in multiple image pairs, resulting in many overlapping measurements of ice velocity, as shown in Fig. 4. Georegistration error of each image leads to some visible disagreement between the overlapping velocity measurements, but despite the noise, a coherent pattern of interannual variability is apparent 70 as the clusters of velocity measurements move up or down from year to year. 3 Method of analyzing image pair velocity data

Assess and remove interannual variability
The first step toward quantifying seasonal variability for any pixel is to remove any interannual variability from the time series. Interannual variability can be determined by smoothing the velocity data using any of several common methods. A 75 polynomial fit is robust and computationally efficient, but requires choosing a polynomial order, which is subjective and can lead to overfitting or underfitting the data. Alternatively, a moving average makes no prior assumption about the shape of the velocity curve, and can adapt to any arbitrary interannual variability. After exploring several approaches, we find the best results by first detrending the time series with a low-order polynomial, then using a hybrid of a moving average and a spline fit, which we describe below.

80
When assessing interannual variability, we temporarily ignore the duration over which each image pair measured ice displacement, and simply assign the average velocity to the center date t m of each image pair. Using the range of times t m in the time series, we require at least two years of velocity data and detrend using a polynomial whose order is chosen as one quarter of the range of t m in years, rounded up to the nearest integer. The result is that for up to four years of data, the time series is linearly detrended; from four to eight years of data, the time series is quadratically detrended, and so on. We assign the velocity weights w v in the polynomial fit using the formal error estimates σ v from the ITS_LIVE data product such that w v = σ −2 v . An example of a fourth-order polynomial fit to the velocity data is shown in Fig. 4.
After detrending the time series with a weighted polynomial fit, we characterize any residual interannual variability with a spline fit to the mean velocities of each year. We take the weighted mean velocity of all measurements whose t m lies within 183 days of winter solstice of that year. We assign the weighted mean velocity of each year to the weighted mean date of those 90 velocities, then interpolate with a shape-preserving piecewise cubic hermite polynomial to obtain a measure of interannual variability corresponding to each image pair's center date t m .
In the method described thus far, most summer image pairs whose ∆t is small contribute much less to the weighted mean annual velocities than do long ∆t image pairs that span winter. This is because velocity error by any feature-tracking algorithm stems strictly from displacement error σ d , while timing error is essentially zero. So although summer offers more image 95 pair measurements, their shorter ∆t values are associated with greater velocity error σ v and therefore they are significantly Center dates tm are shown as dark gray dots for visual clarity. downweighted in the calculation of annual mean velocities. In other words, it is possible that either due to the higher weights of winter velocities or the higher quantity of measurements available during summer, the weighted mean annual velocities could be biased toward one season or the other. We account for this possibility by iteratively solving for interannual and seasonal variability, as we describe later in Section 3.3.

Assess seasonal variability
We characterize the magnitude and timing of seasonal variability as a simple sinusoid that can be applied in the x and y directions separately to build a two dimensional understanding of how ice moves throughout the year. Limitations of the sinusoidal approximation are discussed later in Sect. 6, but we justify our approach as it is the simplest means of describing cyclic behavior, and by nature it is constrained to capture the sum total of ice displacement that occurs throughout the year.

105
After characterizing interannual variability, we subtract it from the velocity time series at the center date of each image pair.
The residuals v r after removing interannual variability can then be assumed to contain only seasonal variability and noise. To address blunders, we remove outliers whose absolute value exceeds 2.5 robust standard deviations of v r (see Appendix A).
The remaining task is to fit a sinusoid to the v r time series, but given that each velocity measurement is a single value that represents several weeks to more than a year, we cannot fit a sinusoid directly to the velocity time series or assume that values 110 of v r correspond to the image pair center dates t m . Instead, we operate on the displacements associated with each image pair, taken as the integrals of velocities v r .
We seek to define seasonal velocity variability v s in the form where A and φ are the amplitude and phase of the seasonal velocity cycle, respectively, t represents time in decimal years, 115 and T = 1 year is the period of oscillation. By our method, the constant C 0 is ideally zero after detrending and removing interannual variability, but we include it here in case some residual offset is present in the v r time series.
For a more robust least-squares solution, we employ a trigonometric identity to rewrite Eq. 1 as which is related to Equation 1 by A = C 2 1 + C 2 2 and φ = atan2(C 2 , C 1 ).

120
Again, we cannot solve Eq. 1 or Eq. 2 directly with image-pair velocities, because they do not represent instantaneous velocity at any known times t. Instead, image pair data track displacements, which equate to the integral of velocities between times t 1 and t 2 when the two images were acquired. Accordingly, we can take the definite integral of Eq. 2 to solve for the seasonal displacement cycle d s as, (3) 125 We solve for the coefficients of Eq. 3 using a least-squares approach with weights given by w d = (σ v · ∆t) −2 . This type of approach has previously been applied to study Earth deformation (e.g., Hetland et al., 2012), and is similar to an approach that has been used to analyze ice dynamic responses to tidal forcing Minchew et al., 2017).
By employing Eq. 3 to solve for A and φ, we obtain a first approximation of the seasonal variability of ice velocity. However, the possibility remains that the initial estimate of interannual variability may have been partly aliased by seasonal variability 130 due to uneven temporal sampling. Thus, we refine our estimates of interannual and seasonal variability by iterative means.

Iteratively refine interannual and seasonal variability estimates
If our first estimates of A and φ are correct, then seasonal variability must have aliased each initial estimate of interannual variability by a certain amount v a . The amount of velocity aliasing equates the seasonal displacement over time, which we can obtain by dividing the definite integral of Eq. 1 by ∆t, or After obtaining initial estimates of A and φ, we subtract v a from the original detrended velocity measurements and repeat the process of calculating interannual and seasonal variability. In most cases, we find that initial estimates of A and φ are reasonably close to their true values, and that just a few iterations are sufficient to yield accurate final estimates of seasonal amplitude and phase. We explore the number of requisite iterations in Sect. 4.3 and show the effects of iterating in Fig. 5.
140 Figure 4 shows an example of our method applied to a time series of 14,208 ITS_LIVE image pairs acquired in a single pixel near the grounding line of Byrd Glacier, Antarctica. Because no strong long-term velocity trends are present, the green fourth-order polynomial fit is relatively flat, albeit with a slight downturn at the unconstrained end of the time series. The blue curve of interannual variability follows the multi-year rises and falls of velocity much more closely, and is uncontaminated by seasonal variability on this tenth iteration of the solution. The red curve adds seasonal variability with an amplitude of 23 145 m yr −1 to the interannual variability, and given that it is a least-squares solution, it represents the only solution that minimizes the misfit between the curve and all available observations. Indeed, while at first glance it may appear visually that the timing of image acquisitions themselves are the only seasonal pattern in the ITS_LIVE data plotted in Fig. 4, close inspection shows a persistent pattern of summer slowdown in the image pair data that becomes especially clear after the 2013 launch of Landsat 8.

Sensitivity analysis 150
To assess uncertainty of our method, we generate synthetic random continuous velocity time series, then artificially sample them using random subsets of image acquisition times from Landsat image pairs that cover Antarctica. We then apply the methods described in Section 3 to determine how accurately seasonal ice dynamics can be measured, and under what conditions.

Synthetic time series generation
Any synthetic velocity time series used for testing should resemble the true nature of ice dynamics in its variability on all 155 timescales. Accordingly, we create realistic interannual variability by matching the distribution of the standard deviations of velocities in each grid cell in the ITS_LIVE annual velocity mosaics from 2013 to 2018. Figure 3 shows a map of interannual variability for the grid cells that contain all six years of data, and a histogram of those values for grid cells that contain at least 100 total image pairs and a minimum mean velocity of 15 m yr −1 . Within this subset, the median interannual variability is characterized by a standard deviation of 4.2 m yr −1 .

160
To create each synthetic time series of interannual variability, we generate uniformly distributed random values centered about zero, at daily temporal resolution. We apply a first-order low-pass Butterworth filter with a cutoff period of 548 days to each random time series to ensure that no annual cycles are present, and then we discard the first and last 548 values of each time series to eliminate edge effects of the filter. We then scale the remaining time series such that its standard deviation matches a prescribed level of interannual variability.

165
Our handling of interannual variability is an attempt to mimic the observable ways in which Antarctic glaciers speed up and slow down from year to year. Ideally, we would also carry forth in a similar manner for seasonal variability, imposing cyclic behavior that matches the true character and distribution of the types of seasonal variability that exist in nature. However, as the intent of this paper is to develop the methods that will be necessary to understand where and how ice velocities vary on seasonal timescales, we cannot at present create synthetic seasonal variability distributions to match what truly exists in nature.
170 Instead, to each synthetic time series we add a sinusoid with a period of 365.25 days, a random phase, and a random amplitude between 0 and 100 m yr −1 .

Synthetic time series sampling
Each synthetic time series is artificially sampled using characteristic acquisition times and error distributions from Antarctic Landsat-derived ITS_LIVE image pairs. In Fig. 3 we see that among grid cells containing at least 100 image pairs, and whose 175 mean velocity is at least 15 m yr −1 , we can expect a median of 1153 image pairs per grid cell. Accordingly, we use the acquisition times (Fig. 2) and corresponding error estimates (Fig. 3) from random subsets of 1153 image pairs (sampled from all 1.8 million Antarctic image pairs) to artificially sample each synthetic time series. Each synthetic velocity measurement is calculated from the cumulative sum of daily velocities that occurred between the first and second image in a given pair, but before dividing the total displacement by the time ∆t between images, a random gaussian value of displacement is added 180 according to the formal error estimates associated with that image pair in the ITS_LIVE dataset. The result is a time series of synthetic measurements that resemble the acquisition times and error characteristics of ITS_LIVE image pairs, but whose underlying continuous velocity signal it represents is known.

Seasonal amplitude and phase recovery
We conducted several tests to determine the accuracy with which we can recover the amplitudes and phases of seasonal cycles 185 in synthetic velocity time series. The parameters of each test are detailed in Table 1.  Table 1 indicate that a-b: amplitude and phase errors level off after the first few iterations of removing interannual variability and solving for seasonal variability; c-d: at least 500 to 1000 image pairs are necessary for the most accurate, stable solutions; e-f: strong interannual variability can drastically affect seasonal amplitude and phase detection; g: seasonal amplitude measurement uncertainty is independent of the seasonal ice velocity amplitude itself, but h: seasonal phase is measured most accurately when the seasonal ice velocity amplitude is strong; and i-j: although some faint residual effects of the seasonal bias in sampling persist after 10 iterations, amplitude and phase accuracy are effectively independent of the timing of the seasonal ice velocity variability.
In the first test, we sought to understand how many iterations are necessary to achieve a stable solution. In this test, we generated 10,000 synthetic velocity time series, each having interannual variability with a standard deviation of 4.2 m yr −1 .
Sinusoidal seasonal variability was added to each time series, characterized by a random phase and a uniform distribution of amplitudes between 0 and 100 m yr −1 . Each time series was then synthetically sampled using the acquisition times and error 190 characteristics of 1153 random image pairs in the ITS_LIVE dataset. Following the method described in Sect. 3, we analyzed each time series by assessing interannual variability in the synthetic measurements, removing it, removing outliers, and then obtaining the amplitude and phase of each seasonal cycle. After accounting for any aliasing that would have been caused by the seasonal amplitude and phase (Sect. 3.3) that were measured in the first iteration, we conducted a second iteration of assessing interannual and seasonal variability. We then followed with a third iteration, and so on, up to 15 iterations. Results in Fig. 5a, b 195 show that solutions tend to converge after the first few iterations, but it is possible that in some situations more iterations could be necessary, depending on sampling and the characteristics of the velocity time series. Accordingly, we use 10 iterations in all of the tests that follow.
We conducted a second test to determine how many image pairs are necessary to detect seasonal cycles in a velocity time series. We found that with just 32 image pairs, we were able to recover phase with sufficient accuracy to correctly identify the 200 season of maximum velocity (Fig. 5c,d). We also found that performance improves dramatically with increasing image pairs, until errors in phase and amplitude approach their asymptotes between 500 and 1000 image pairs. Beyond that, the number of image pairs used in the analysis has a negligible impact on the accuracy of seasonal signal detection.
In a third test, we found that one factor more than any other threatens the accuracy of seasonal variability detection. Given a velocity time series sampled by 1153 image pairs, seasonal amplitude error increases approximately linearly with the level of 205 background interannual variability (Fig. 5e,f). This is because despite our attempts to account for interannual variability, it is inevitable in a time series of finite length that some residual variability will influence the least-squares fit of the seasonal cycle.
Nonetheless, if we consider ±45 days of phase uncertainty to be a threshold that indicates accurate detection of the season of maximum velocity, our method performs adequately in the presence of interannual variability with a standard deviation exceeding 100 m yr −1 . Further analysis suggests that for any level of interannual variability, the amplitude of the seasonal 210 cycle must be at least one third of the standard deviation of interannual variability to reliably be detected (see Sect. A1). We note, however, that because we have used the simple metric of the standard deviation of velocity to describe all forms of interannual variability, it is likely that our method will perform better than these error estimates suggest wherever interannual variability is dominated by a long-term trend.
The final four panels (g-j) of Fig. 5 provide valuable insights into the capabilities and sensitivities of our seasonal detection 215 method. Most notably, that for a constant level of background interannual variability, seasonal amplitude errors are independent of the amplitude of the seasonal cycle. However, phase detection benefits with increasing seasonal amplitude. This should not be surprising, as a signal must not only exist, but also be sufficiently strong for its phase to be accurately measured. The last two panels (i-j) of Fig. 5 show that for all practical purposes our seasonal amplitude and phase estimates can be considered insensitive to the timing of the underlying seasonal signal, despite having no images during winter.

220
Thus far we have determined that 10 iterations are more than sufficient for the model of seasonal variability to converge. We have also found that our measurements can be considered agnostic to the timing of ice velocity variability and to the timing of satellite image acquisitions, when applied to at least 500 to 1000 image pairs. With this understanding, we now apply our method to 100,000 synthetic time series that typify the interannual ice velocity variability and satellite image acquisitions that have been measured across the Antarctic continent as a whole. Interannual variability in the synthetic time series was randomly 225 sampled from the distribution shown in Fig. 3d and a uniform distribution of seasonal amplitudes from 0 to 100 m yr −1 were added to the interannual variability. Each synthetic ice velocity time series was observed with a number of synthetic image pairs taken randomly from the distribution of values in Fig. 3c that exceed 500 image pairs, and displacement errors randomly sampled from the distribution in Fig. 3e were added to each synthetic measurement before dividing by the times ∆t that separate each image pair. The results of this test, shown in Fig. 6, indicate that with the Landsat images that have been acquired 230 to date, we should be able to measure seasonal amplitudes within ±1.4 m yr −1 and phase within about ±2 days.

Validation with in situ GPS observations
To confirm that we can extract seasonal ice velocities from feature-tracking data, we compare an in situ GPS position time series against results obtained by applying our method to ITS_LIVE velocities that were generated from 15 m resolution  Table 2. data, so we directly apply the methods discussed above without any adjustments. Likewise, our methods could just as easily be applied to ITS_LIVE data covering Alaska, the Canadian Arctic, Svalbard, the Russian Arctic, Patagonia, or High Mountain Asia. We compare the GPS velocity time series to velocities obtained from 5189 ITS_LIVE image pairs in the pixel closest to the fixed median location of the GPS receiver throughout the course of its data collection. Key findings are listed in Table 2.
250 Secular mean velocities in the x and y directions measured by the two independent datasets agree on the order of 1 m yr −1 .
Seasonal velocity amplitudes also show agreement on the order of 1 m yr −1 , as we expect on this glacier where interannual variability has a standard deviation of 3.6 m yr −1 in the x direction and 3.4 m yr −1 in the y direction (see Appendix A2). Table 2. ITS_LIVE analysis compared to GPS. The secular mean velocities and x and y components of the seasonal characteristics of the glacier velocity time series shown in Fig. 7. The acquisition dates of the two independent measurements differ slightly, yet agreement is within a few percent by every measure. Uncertainty estimates are described in Appendix A2. The largest disagreement in Table 2 lies in the phase of the x direction, which at 14 days is only about 4% of a year. We note further that disagreements listed in Table 2 do not necessarily reflect errors in the ITS_LIVE data or the methods we have 255 employed to process it. Rather, disagreements may simply reflect that the GPS receiver did not record data from late 2010 to mid 2012, and no GPS data were recorded during two of the winters that followed. Meanwhile, we use ITS_LIVE data from images that were collected more than a year before the start of the GPS record, but image pair data have not yet been processed corresponding to the final months of the GPS record. In addition to these differences in timing of data collection, we also note that the GPS record represents a Lagrangian measurement of ice velocity, whereas the ITS_LIVE data record an Eulerian 260 measurement at a pixel that the GPS receiver passed through only once in its decade-long life. Nonetheless, agreement between the GPS solutions and our method of image pair data analysis is notable, and lends credence to this method as a precise means of measuring seasonal variability in ice velocity.

Discussion
The method we present in this paper requires a multi-year record with at least several hundred image pairs to confidently 265 identify the amplitude and phase of seasonal variability. Some key regions of interest, such as Totten Glacier in Antarctica, do not yet offer enough cloud-free images to meet this threshold, so a few more years of data collection may be required before our methods can be successfully implemented there. In addition, it may be difficult or impossible by our methods to detect seasonality in places of interest such as Pine Island and Thwaites glaciers, which are currently undergoing dramatic interannual change (Mouginot et al., 2014) that could confound our measurements of seasonal variability.

270
Nonetheless, Fig. 7 illustrates the sensitivity of our method to tiny variations in glacier flow when conditions are favorable.
The vertical scale of the upper panels spans just ±15 m; yet, by our method of analyzing 15 m resolution Landsat Band 8 images, we are able to detect minor nuances in position that occur entirely within this range. We know of no other sensor or dataset that can offer such insights into these kinds of subtle variabilities of ice flow that have occurred over the past few decades.

275
The most significant limitation of the method we present may lie in our approximation of seasonal variability as a sinusoid.
Where seasonal dynamic variability has previously been documented, it has been found that in some cases a sawtooth function or other higher-order fits might match seasonal variability more closely than a sinusoid (Moon et al., 2014;Van de Wal et al., 2015;Vijay et al., 2019). In particular, although events such as springtime calving or summer melt may occur on an annual cycle, a glacier's complete response to them may only last for only a few days (Schoof, 2010). To approximate such events as 280 smoothly varying sinusoids will underestimate the magnitude of any brief glacier acceleration, while potentially giving a false sense that the ice responds to an impulse event continually throughout the entire year.
Despite the tendency of a sinusoid to oversimplify complex time series, we contend that no other description of seasonal variability is as elegant or robust over decadal timescales, and that understanding seasonal ice dynamics must begin with a zeroth order description of the amplitude and phase of ice velocity variability. While it is true that a glacier can accelerate in 285 response to a transient event and return to an equilibrium velocity within just a few days (Stevens et al., 2015;Andrews et al., 2014), in the climatological sense, nature does not consistently time such events as calving or increases in basal water pressure with any greater precision than the method we have presented to detect them. Short-lived speedup events do not consistently occur on the same day each year, so the decadal average tends to be smoother and more sinusoidal than the velocity profile of any single year.

290
The approach we present captures the fundamental mode of subannual variability in ice velocity and conserves all ice displacement that occurs throughout the year. The simple two-term explanation of amplitude and phase provides a robust description of seasonality that is less prone to error than higher-order fits such as three-term sawtooth functions. By providing a simple measure of amplitude and phase, we offer a straightforward method to compare how neighboring glaciers respond to a common seasonal forcing or investigate how dynamic signals propagate upstream or downstream in a given glacier.

295
In Sect. 5 we applied our method to ITS_LIVE velocity measurements and compared against GPS position data, with the intention that the in situ GPS station would provide a reliable ground-truth reference. However, although the two independent measurements show close agreement whenever GPS data were logged, the receiver's harsh polar environment has led to several long gaps during which no GPS position data were acquired. This suggests that as a means of measuring ice dynamic climatology, our method might not only meet, but exceed the performance of the in situ GPS receiver while providing insights 300 into dynamic behaviors that occurred years before the station was installed. The particular ITS_LIVE grid cell in Greenland that we analyzed in this paper is hardly unique in its potential to provide historical context, as decades-long velocity records exist in most of the 240×240 m pixels which cover nearly all of the world's ice.
The methods presented in this paper have focused primarily on optical satellite data because no other type of sensor provides such a long record of ice velocity. As more radar data become available, particularly since the launch of Sentinel 1a/b, the 305 problem of missing winter data will be eliminated, but the methods presented in this paper will still hold. When an abundance of feature-tracked velocities from radar become available, Eq. 3 may then be easily modified to include additional terms to simultaneously solve for cyclic velocity variability on tidal timescales as well as seasonal variability.

15
Given a relatively continuous time series of at least 500 to 1000 image pairs, our method can extract the amplitude of seasonal 310 variability with a precision on the order of about 1 m yr −1 by describing ice flow as an oscillation, provided the level of background interannual variability does not overwhelm the overall signal. We find that if the amplitude of seasonal variability is at least a third of the standard deviation of interannual variability, the method we describe can reliably detect the season of maximum ice velocity. Ability to detect seasonal amplitudes is independent of the amplitude and phase of the seasonal velocity variability, but phase accuracy benefits with increasing amplitude of seasonal variability.

315
With the method we describe, we may begin to map seasonal ice dynamic variability on a global scale, in a consistent and meaningful manner; and by providing a method that can be employed independently in the two dimensions of Cartesian coordinates, we hope to gain a more complete understanding of how dynamic signals propagate through the world's ice.
Code and data availability. Data analysis in this paper relied upon Antarctic Mapping Tools for MATLAB (Greene et al., 2017b) and the Climate Data Toolbox for MATLAB (Greene et al., 2019). All data analyzed in this paper and code necessary to generate the figures are 320 included as a supplement to the paper, or are available online at http://www.github.com/chadagreene. PROMICE data is freely available at http://www.promice.dk. ITS_LIVE global image-pair velocity data is freely available at its-live.jpl.nasa.gov.

Appendix A: Additional uncertainty analysis
Phase and amplitude uncertainty in Sect. 4.3 were calculated from the differences between imposed and recovered values.
Differences in phase were wrapped such that their absolute values did not exceed 182.62 days. During this work, we found 325 that in some cases, simple standard deviations of differences could be strongly influenced by a few extreme outliers, and as a result, standard deviations of differences (and similarly, root-mean-square errors) did not reflect the true gaussian distributions of errors. So for a more robust and meaningful measure of error distributions, we define σ as 1.4826 times the median of the absolute value of differences. uncertainty contour marks an approximately linear relationship between interannual variability and seasonal amplitude. Indeed, the dashed blue line that nearly coincides with the 45 day phase uncertainty contour shows the simple slope whereby 340 seasonal amplitude is one third of the standard deviation of interannual variability, with zero offset from the origin. This tells us that there is a quite simple relationship between the amplitude of a seasonal cycle, the level of background interannual variability, and our ability to detect phase-As long as the seasonal amplitude is above the noise floor of about 1 m yr −1 and interannual variability does not exceed three times the amplitude of the seasonal cycle, we can accurately detect the seasonal cycle.

A2 Uncertainty estimates for ITS_LIVE/GPS comparison
The values of seasonal amplitude and phase uncertainty listed in Table 2 were each calculated as the robust standard deviations of amplitude and phase errors from 5000 synthetic time series sampled by 5189 image pairs. In the x direction, ITS_LIVE image pairs indicate an interannual variability of 3.85 m yr −1 in our pixel of interest, with a seasonal amplitude of 28.8 m yr −1 and a maximum velocity occurring around December 13th of each year. We generated 5000 synthetic time series matching 350 these criteria, and were able to recover the 28.8 m yr −1 seasonal amplitude within σ=0.87 m yr −1 and phase within σ=4.7 days. Similarly for the y direction, from 5000 synthetic time series with an interannual variability of 3.38 m yr −1 , seasonal amplitude of 14.4 m yr −1 , and maximum velocity occurring on January 1 of each year, we recover seasonal amplitudes within σ=0.91 m yr −1 and phase within σ=10.7 days.