Measuring the state and temporal evolution of glaciers using SAR-derived 3D time series of glacier surface flow

SAR-derived 3D time series of glacier surface flow Sergey Samsonov1, Kristy Tiampo2, and Ryan Cassotto2 1Canada Centre for Mapping and Earth Observation, Natural Resources Canada, 560 Rochester Street, Ottawa, ON K1S5K2 Canada 2Earth Science & Observation Center, Cooperative Institute for Research in Environmental Sciences, University of Colorado, Boulder, CO 80309 USA Correspondence: Sergey Samsonov (sergey.samsonov@canada.ca)

Glacier is the world's largest piedmont glacier covering approximately 2200 km 2 on the flat coastal foreland (Sharp, 1958;Muskett et al., 2003;Sauber et al., 2005) and is partially fed by Seward Glacier, a surge-type glacier that originates in the upper reaches of Mt. Logan (Sharp, 1951;Ford et al., 2003). A mass budget deficit in the Malaspina-Seward complex has long been recognized (Sharp, 1951). Agassiz Glacier is another surge-type that flows in an adjacent sinuous valley northwest of the Malaspina-Seward complex (Muskett et al., 2003;Sauber et al., 2005). The Klutlan Glacier is an 82-km long surge-type valley glacier located at elevations between 1300 and 2100 m; it has surged repeatedly over the last few hundred years (Wright, 1980;Driscoll, 1980). Walsh Glacier is a 90-km long surge-type valley glacier located at a higher elevation of about 1500-3000. It is fed by two major branches, one from the north and one from the east and converges with the Logan Glacier downstream (Fu and Zhou, 2020). Our results demonstrate that this new technique can be used to analyze 3D flow velocities of glacier surfaces over large regional scales using nearly 30 years of archived SAR data and for near real-time monitoring of these glaciers us-75 ing rapid revisit SAR data from satellites, such as Sentinel-1 (6 days revisit period) and forthcoming NISAR (12 days revisit period).

Data and Model
We downloaded 198 ascending (track 123) and 214 descending (track 116) Sentinel-1 single-look complex (SLC) images from the NASA Distributed Active Archive Center (DAAC) operated by the Alaska Satellite Facility (ASF) ( Table 1). Two 80 ascending and two descending frames along the azimuth directions were concatenated for each, resulting in 99 and 107 swaths, respectfully. Ascending and descending sets were processed individually using GAMMA software (Wegmuller and Werner, 1997) that produced range and azimuth offsets for consecutive pairs (Figure 3). To compute offsets, we used a 64×16 pixels sampling interval and 256×256 pixels correlation window. Such a large window was required to obtain a distinct peak of the 2D cross-correlation function. Offsets were spatially filtered using a median filter with a 2 km (6 sigma) filter-width, geocoded 85 using TerraSAR-x 90 m DEM and resampled to a common grid with a ground resolution of 200 m.
The 3D displacement time series were computed by inverting a set of linear equations, first solving for the north, east and vertical flow velocities V n,e,v for each acquisition epoch and then for cumulative 3D flow displacements D n,e,v Here, in a matrix form, RO are the range and AO are the azimuth offsets computed from Sentinel-1 data; L is the Tikhonov 90 regularization matrix multiplied by the regularization parameter, λ, and A is the transform matrix constructed from the time intervals between consecutive SAR acquisitions and the range S r and azimuth S a directional cosines with north, east, and vertical components where φ is the azimuth and θ is the incidence angles.
This work is now possible due to the improved availability over large areas of high-quality, high-resolution, temporarily-dense ascending and descending SAR data and the increase in computational power that allows computing a large amount of range 110 and azimuth offset maps and inverting large matrices. Since this method does not make any assumptions about the direction of motion, it provides the optimal solution applicable to any phenomenon (e.g. glacier flow, tectonic and anthropogenic deformation, etc). The typical size of the transform matrix A exceeds 100s and often 1000s of columns and rows for each coherent pixel. It is 609×1014 in our case. The singular value decomposition (SVD) algorithm from the Linear Algebra PACKage (LA-PACK) library called from C++ code is used for inverting this matrix for each coherent pixel. Depending on the number of 115 cores in the processing unit and the number of pixels, this process can take from several hours to several days.
First-order Tikhonov regularization was applied during the inversion, resulting in temporal smoothing. The magnitude of smoothing is controlled by the regularization parameter λ that can be selected, for example, using the L-curve method (Hansen and O'Leary, 1993;Samsonov and d'Oreye, 2017). The choice of a regularization order is not critical in our case; for example, the second-order Tikhonov regularization produces visually indistinguishable results.

120
An estimate of the 3D mean linear flow velocity was computed from cumulative flow displacements. A straight line was fit to the 3D flow displacement time series and then divided by the total length of time of the record, here three years. For aesthetic purposes, horizontal flow velocity was resampled to a coarser resolution and values less than 5 m/year were removed. Areas of interest were reduced to three small sub-regions (AOI1, AOI2 and AOI3 in Figure 1) that cover only five glaciers.

125
The SAR intensity images (   erence frames is well established in the glaciological community, it is less known in the solid-earth geophysical community (Samsonov and Tiampo, 2006;Samsonov et al., 2007;Gourmelen et al., 2010;Shen and Liu, 2020). For the latter, GNSS and SAR measurements can be considered nearly identical only when flow velocity at a fixed geophysical location (i.e. SAR pixel) is equal to the flow velocity at the GNSS site; that is when GNSS-derived displacements are contained within a single SAR pixel. This occurs when the flow velocity is very small and the material is rigid, as in tectonic deformation studies. For rapidly 195 deforming glacier surfaces, the differentiation is far more critical. Although the time-series in Figure 7 resemble GNSS-derived displacements, it is important to remember that these Eulerian measurements represent the cumulative displacement at any one pixel over time. Hence to emphasize this difference, we use the flow displacement terminology. The downward flow is expected in the upper reaches of accumulation zones because of firn compaction, and in areas with steeply dipping surfaces due to sloping bed topography; however, downward flow in the lower ablation zone is more concerning.
In general, the accumulation of snow and ice in high elevations produces a net mass gain that replenishes ice lost through ablation processes along the lower glacier. In steady-state, these processes balance each other and lead to submergent flow Velocities along Klutlan Glacier vary in more complex and interesting ways with multiple zones of upward and downward flow observed. This surge-type glacier (RGIConsortium, 2017) has a 30 (Meier and Post, 1969) to 60-year surge cycle (Wright, 225 1980;Driscoll, 1980  can only be derived from side-looking SAR measurements that capture horizontal and vertical components of motion. SAR measures glacier motion at a certain depth rather than at the surface. Previous studies for this region suggest that the C-band SAR penetrates about four meters into the glacier's firn layer in dry conditions (Rignot et al., 2001). The penetration depth is affected by water content and thus varies by season; it is likely that some fraction of the signals observed in Figures   245 4-7 is due to the seasonal fluctuation in the penetration depth. A combination of SAR and optical data in a complementary manner could provide useful information about the variation of flow velocity with depth throughout different seasons.

Conclusions
We presented a novel flow displacement technique to observe variations in glacier surface flow in 3D using ascending and descending SAR scenes. The 3D flow displacement (and/or velocity) time series computed allowed us to map in unprecedented