Seven-year variation in glacier surface velocity at Narsap Sermia in Southwest Greenland

ABSTRACT

Glacier velocity is a critical factor closely related to ice dynamics of the Greenland ice sheet (GrIS). We present variations of ice velocity using space-based synthetic aperture radar (SAR) observations from July 2013 to September 2020. We focus on a tidewater outlet glacier of Narsap Sermia (NS), located in southwest Greenland, which is experiencing significant ice loss. A time series of high spatiotemporal resolution maps of ice velocity was traced by using an offset tracking approach with seven-year SAR acquisitions to understand glacier dynamics in time. We observed seasonal variations showing a dramatic increase in late spring (May) and early summer (June) followed by a sharp decrease in late summer (September). The seasonal variation was confirmed through Seasonal and Trend decomposition using Loess (STL) analysis. Overall, the NS glacier experienced a gradual speed-up of ice velocities. The ice velocity increased by 0.1–0.3 km/year from 2015 to 2017, then after a temporary decrease in 2018, the velocity moved faster by 0.9 km/year in 2020. The speed of glacier movement accelerated as it approached the tidal terminus of the glacier, showing approximately 6.0 km/year in 2020. The ice speed in 2020 increased by 0.9 km/year at the terminus region compared with 2014. We also discuss the relationship between the glacier velocity and temperature data such as air, land surface, and sea surface temperatures.

KEYWORDS:

• Ice velocity
• Narsap sermia
• offset tracking
• SAR
• glacier
• Time series analysis

1. Introduction

The Greenland ice sheet (GrIS) is among the largest ice bodies in the world and plays a significant role in the global climate system (Robinson, Calov, and Ganopolski 2012; Vizcaino et al. 2015; Pattyn et al. 2018). Cumulative ice loss of 3,900 billion tonnes from 1992 to 2017 in Greenland has caused an approximate 11 mm rise in global sea level (Bamber et al. 2018). The Intergovernmental Panel on Climate Change declared that the melting of the GrIS during the 21^st century is set to contribute to the global sea-level rise by 4–27 cm, which threatens populated coastal cities around the world (Shepherd et al. 2020; Briner et al. 2020; Pörtner et al. 2019). The ice surface velocity is a hypothetical flow velocity estimated at the given physical phase of the glacier. Ice velocity is a crucial factor in studying ice dynamics, ice sheet mass balance, and numerical ice sheet models (Teng et al. 2018; Cheng and Lötstedt 2020; Barnes and Gudmundsson 2022) because the ice velocity is mainly dependent on ice thickness, the slope of the surface, properties of basal sediments, and ice temperature as well as weather and seasonal conditions. A time series of glacier velocity maps with the high spatial and temporal resolution is beneficial for estimating glacier mass loss, calculating glacier volume, and understanding the relationship between glacier systems and climate environments (Vijay et al. 2019; McNabb et al. 2012). Various researchers have reported that the glacier velocity in Greenland varies seasonally using remote sensing observations (Joughin et al. 2014; Ahlstrøm et al. 2013; Bevan, Luckman, and Murray 2012). However, as most previous studies have focused on spring and summer velocities to examine driving forces of its speed, seasonal variation in glacier velocity has not yet been fully clarified. Ice velocities with a high spatial and temporal resolution can be utilized to estimate discharge and mass flux of glaciers and to analyze glacier behavioral patterns that might depend on weather conditions and seasonal effects. A time-series of glacier velocity maps with high temporal and spatial resolution can be helpful to understand the local glacier dynamics and seasonal ice flow variation. Furthermore, it may be possible to also use the high temporal resolution maps to evaluate the impact of local temperature changes. Frequent or continuous observations are a key to constructing a high temporal glacier velocity map. A terrestrial measurement, such as the Global Navigation Satellite System, ground-based laser, or radar survey, can estimate the most accurate ice velocity with a very high temporal resolution. However, they cannot produce a two-dimensional high spatial resolution map, which is helpful to trace the history of glacier movement over a wide area. Optical observations by airborne and unmanned aerial vehicles have been attractive methods to capture the earth’s surface with a high spatial resolution map. However, their practical applications are often limited by operational costs and logistics. Polar orbit earth observation satellites can revisit high latitude regions more frequently; thus, they have been widely utilized to calculate the ice velocity (Han, Hee Kim, and Kim 2022; Sahu and Gupta 2021; Gómez et al. 2020; Yellala, Kumar, and Arild Høgda 2019; Vijay et al. 2019; Mouginot, Rignot, and Scheuchl 2019; Lemos et al. 2018; Bhushan et al. 2018; Rignot and Mouginot 2012; Joughin et al. 2010). Active radar imaging using microwave frequency like synthetic aperture radar (SAR) is preferred because optical observations are limited by changing weather conditions.

Differential SAR interferometry (DInSAR) using active-imaging synthetic aperture radar has been widely used to monitor various surface displacements such as earthquakes, volcanoes, ground subsidence, water level change, including glacier displacements (Massonnet et al. 1993; Rignot 1996; Massonnet and Feigl 1998; Amelung et al. 1999, 2000; Hong et al. 2010). Although the DInSAR application can detect relative surface displacements from cm to mm accuracy on the solid surface, severe decorrelation in fast-moving bodies like a glacier or lava flow prohibits displacement calculations (Hanssen 2001; Strozzi et al. 2006; Riveros et al. 2013). Offset tracking is an alternative approach for estimating transitional displacement correlating with the same feature between two observations (Waechter, Copland, and Herdes 2015; Schellenberger et al. 2016; Wychen et al. 2018). The accuracy of the offset tracking method depending on the pixel spacing rather than the wavelength of the radar signal is slightly degraded compared to the few cm or mm accuracy of the DInSAR technique (Strozzi et al. 2002; Nagler et al. 2015). The degradation of correlation between two SAR images could often occur in the summer because wet surfaces due to melting might affect the backscattered amplitude. Furthermore, the window patch size used to evaluate the correlation in tracking offset is closely related to the spatial resolution of the velocity map. Additionally, intensive computation power would be required to calculate the offset with a large window patch size for a vast study area. Despite its disadvantage, offset tracking has been intensively utilized to estimate glacier velocity for cryosphere research due to its robustness.

Tidewater glaciers that terminate in the sea produce icebergs through calving by fast ice flow in Greenland. Twila et al. (2014) proposed three types of tidewater glaciers with their speed change depending on the amount of seasonal meltwater (Twila et al. 2014). Type 1 glacier is characterized by sustained speedup even at the end of the melt season and acceleration from late spring (May) to early summer (June). Type 2 glacier shows speedup in late spring (May) and early summer (June), and slowing in the middle of summer (July–August), and returns to pre-acceleration speed during the winter season (December–March). Type 3 glacier slowdown in speed in the middle of summer (July–August), reaching a minimum speed by late summer (September). Research on variations in the tidewater glacier velocity related to subglacial drainage and channel formation has also been undertaken (Motyka et al. 2017; Davison et al. 2020).

A few researchers have reported that Narsap Semia (NS) experienced 3.3 km of a large calving retreat into an over-deepened basin from 2010 to 2014, then stabilized at a retreating rate of 150 m/year between 2013 and 2015 (Motyka et al. 2017). The glacier speed of 1.4 km/year was one of the fastest in the southwest of Greenland from the International Polar Year 2008–2009 dataset (Rignot and Mouginot 2012). Davison et al. (2020) proposed that in three tidal glaciers adjacent to NS, the development of hydrologically efficient glacial channels by surface-derived meltwater controls ice velocities and variations. They indicated that rapid glacier retreat depends on seasonal variation considering hydraulic conditions, as observed through Sentinel-1 SAR and Landsat-8 optic satellite data (Davison et al. 2020). Although they generated seasonal glacier velocity maps using multi-temporal satellite observations, the temporal resolution is insufficient to show a detailed variation of the ice flow. In this study, we analyzed the relationship between the ice velocity and temperature of NS and examined trends and seasonality through multi-temporal ice velocity maps using very high spatial and temporal resolution velocity maps.

Here, we generated a time series of glacier velocity maps over the course of seven years (2013–2020) on NS, with a high temporal resolution using vast 325 SAR observations from TerraSAR-X (TSX)/TanDEM-X (TDM), 6 COSMO-SkyMed (CSK), and Sentinel-1 satellites. Seasonal and Trend decomposition using Loess (STL) method was applied to evaluate seasonal variation. We also examined the correlation between the variation of ice velocity and temperature with daily air temperature data, sea surface temperature (SST), and land surface temperature (LST) products from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor. We first describe the dataset and methodology used for the glacier velocity map and then present variations in glaciers using time-series maps. We present a correlation between the glacier velocity and various datasets of temperatures. Finally, we discuss the advantages and limitations of the offset tracking method for glacier research and summarize our study.

2. Study area

Narsap Sermia is one of the several large tidewater outlet glaciers in the southwestern margin of the GrIS that is connected to Godthåbsfjorden, a 190-km-long sill fjord along with Kangiata Nunaata Sermia and Akullersuup Sermia, as shown in Figure 1(a) (Mortensen et al. 2014). The NS is the closest marine-terminating glacier to Nuuk, Greenland’s capital city, and often calves icebergs into Nuuk Fjord (de Vries, James, and David 2022; Meier and Post 1987). Terrain elevation above sea level of NS is about 444 m, and NS can be affected by the groundwater from warm and saline Irminger Current (Holland et al. 2008). The Godthåbsfjorden water forms a relatively warm intermediate water mass between the sills near Nuuk, extending until approximately 111 km from NS (Mortensen et al. 2013). The grounding line depth of NS is approximately 160 m and the tidewater glacier sends its meltwater and glacier ice to the fjord (Lorenz et al. 2017). As the water depth of the tidewater glacier terminus increases, the calving speed increases, and the calving glacier speed increases as it approaches the ice front (Brown, Meier, and Post 1982; Pelto and Warren 1991). The terminus of the glacier in 2020 retreated by approximately 1.1 km compared to the location of the terminus in 2013 on the centerline of glacier, as shown in Figure 1(d). We used terminus positions derived from optical images from 2013 to 2020 from the NASA MEaSUREs dataset (https://nsidc.org/data/nsidc-0642/versions/2: (Joughin et al. 2021). The terminus position data were obtained between February 3^rd–5^th of each year. Nominal ground resolution of the mosaics and imagery ranges from 15 m to 50 m and nominal uncertainty in terminus position is 50 m. To derive the terminus position, Sentinel-1A/B, RADARSAT-1, Landsat 1–5, Landsat 7, and Landsat 8 were used (Moon and Joughin 2008; Joughin et al. 2021).

Figure 1. a) Location map with polar stereographic coordinate (EPSG 3413) of Godthåbsfjorden, NS and Nuuk in Greenland. b) Geocoded average amplitude image of COSMO-SkyMed SAR observations showing study area. Selected points (P1-P6) for analysis of multiyear ice velocity were marked with various color schemes. The dashed line is a centerline of the NS along the selected points. c) Variation of the annual terminus location of the NS using a product of the NASA MEaSUREs program (Joughin et al. 2021). d) a relative variation of the terminus retreat of NS with respect to the terminus position in 2013.

3. Data and methods

3.1 Data

We utilized three datasets from different SAR payloads to evaluate consecutive glacier velocity for multi-year. Three datasets comprise 95 TerraSAR-X (TSX)/TanDEM-X (TDM) single look slant range complex (SSC) images collected from 22 July 2013 to 9 February 2017, six COSMO-SkyMed (CSK) single look complex slant data acquired from 9 April 2015 to 6 July 2015, and 224 Sentinel-1 interferometric wide mode single look complex images from 4 January 2017 to 15 September 2020 (Table 1). The pixel spacings of the TSX/TDM, CSK, and Sentinel-1 images are 0.9 m × 1.9 m, 1.1 m × 1.9 m, and 2.3 m × 13.9 m in the range and azimuth directions, respectively. Two X-band satellites of TSX/TDM and CSK provide higher spatial resolution, which is useful to estimate the detailed variation of glacier velocity compared with the C-band Sentinel-1 imagery. Slightly different incidence angles ranging from 21.2° to 39.2° were used to observe the study area. The repeat cycle of the CSK and Sentinel is 16-days and 12-days, respectively. However, the revisit time is significantly reduced to 1–8-days and 6-days, respectively, due to their constellation missions. Although the TSX/TDM operate as constellation missions, they usually collect data using either bistatic or pursuit monostatic mode, with an insufficient temporal period to detect glacier velocity (0–10 secs). Frequent observations with short revisit times have the benefit of examining the temporal variation of the glacier velocity, which might be affected by seasonal condition.

Table 1. Characteristics of the SAR observations used in this study.

	TerraSAR-X/ TanDEM-X	COSMO-SkyMed	Sentinel-1A/B
Azimuth pixel spacing	1.94 m	1.93 m	13.89 m
Range pixel spacing	0.91 m	1.14 m	2.32 m
Flight direction	Ascending	Descending	Ascending
Period	22 July 2013- 9 February 2017	9 April 2015- 6 July 2015	4 January 2017- 15 September 2020
Revisit time	11 days	16 days (1, 3, 4, 8 days in constellation case)	12 days (6 days in constellation case)
Number of images	95	6	224
Incidence angle	21.4 °	32.2 °	39.2 °
Carrier frequency (wavelength)	X-band (3.1 cm)	X-band (3.1 cm)	C-band (5.6 cm)

We used the Greenland Mapping Project (GIMP) digital elevation model (DEM) for topographic sensitivity analysis with glacier velocity. The GIMP DEM was generated with a combination of the Advanced Spaceborne Thermal Emission and Reflection Radiometer and Satellite Pour l’Observation de la Terre-5 DEM for ice sheet margin and Advanced Very High-Resolution Radiometer photoclinometry in the interior ice sheet, and the average of ICESat Laser Altimeter System (GLAS) collected from 2003 to 2009 was used to calibrate vertically and horizontally (Howat, Negrete, and Smith 2014). The spatial resolution of the DEM is 30 m and 90 m, and we utilized the 30 m resolution DEM, which was provided as a complete set of 36 tiles in Greenland. Validation using ICESat shows that the mean square root error is ±10 m with ±1 m error on the majority of the glacial surfaces and ±30 m error in large undulating relief areas.

We also examined the correlation between the variation of ice velocity and temperature with daily air temperature data, SST, and LST products from MODIS sensors. The air temperature data were collected at Nuuk, the nearest city to NS (https://en.tutiempo.net/climate). The 8-day products of SST and LST data with 4 km of spatial resolution from 20 July 2013 to 20 September 2020 were collected. The 1-day products could not be used because of a large number of missing values due to the cloud effect in the Arctic. The MODIS temperature products corresponding to the area of interest were spatially averaged to compare with the ice velocity.

3.2 Methodologies

3.2.1 Glacier velocity with offset tracking

We utilized a commercial radar interferometric Gamma software to apply an offset tracking approach to construct a time series of glacier velocity (Werner et al. 2001). The offset tracking method calculates the cross-correlation between two SAR intensity images to track the offset showing the same feature in images collected at different times (Figure 2). The coherence tracking approach to finding the maximum interferometric-phase coherence value was not considered because the rapid movement of glacier surfaces leads to temporal decorrelation, which then degrades interferometric coherence (Strozzi et al. 2002). The SAR observations of each dataset were coregistered to a reference image before applying the offset tracking. Prior to coregistration, unlike the TSX/TDM or the CSK dataset, the Sentinel-1 SAR observations have been debursted to remove an artifact in the azimuth direction and deramped to suppress azimuth-phase ramp, created by its Terrain Observation with Progressive Scans (TOPS) burst acquisition mode (Wegnüller et al. 2016). For offset tracking applications to estimate glacier velocity, a suitable window patch size should be selected by considering glacier seasonal speed. Larger window patch sizes might often be required for faster moving surfaces to track same features in both SAR observations. The window patch size might also be dependent on the time span between two observations because the glacier moves much further away over a long time period. We evaluated the window patch size from 32 × 32 to 512 × 512 for three kinds of SAR datasets to choose the optimum size to estimate the glacier velocity considering spatial and temporal resolutions. Hence, we applied 256 × 256 pixels for both TSX/TDM and CSK datasets and 200 × 40 pixels for Sentinel-1 in the range and azimuth directions. However, 512 × 512 window patch size was also applied to track a large glacier movement for a few longer temporal baseline cases. The estimated offset in pixels could be converted into displacement in meters, which could be translated to the glacier velocity. We utilized consecutive temporal baselines (mostly from 6 to 33 days) in each SAR dataset to apply the offset tracking method. A multi-temporal dataset with a much higher temporal baseline can provide more detailed temporal and seasonal variations of the glacier velocity. Although insufficiently short or long temporal baselines between the two SAR observations are not good enough to track the offset feature, the minimum of 6-day temporal baseline in the dataset could show fast moving glacier velocity in the NS. We adopted 32 × 32 pixels of window steps for TSX/TDM and CSK datasets and 20 × 4 pixels for Sentinel-1 in the range and azimuth, respectively. Since the offset tracking displacement map has different temporal baselines, we converted the estimated offset into velocity per year by considering the time span between two observations, which was used to calculate the annual glacier speed. Finally, three time series were geocoded to have the same geometric coordinates using the GIMP DEM product. Six points were selected to examine the variation of the ice velocity according to location in time along the glacier centerline (Figure 1(b)). We selected those points randomly along the centerline of the glacier. Three points are chosen to show profiles of the faster glacier movement and the others for the slower ice movement. We decomposed the time series of the selected points into three components: trend, seasonality, and residual, using the STL method (Cleveland et al. 1990). The STL is one of the most widely used additive decomposition models and is a robust approach for handling abnormal outliers. The STL works for different seasonalities, such as monthly or quarterly periodicities (Theodosiou 2011; Hyndman and Athanasopoulos 2018). A multi-year time series of glacier velocity clearly showing annual seasonality has been decomposed to evaluate the trend.

Figure 2. Flowchart to estimate glacier velocity using the offset tracking approach with three datasets of SAR observations composed of TSX/TDM, CSK, and Sentinel-1. Considering the ice speed dependent on the period between two observations, adjustable window patch size according to spatial resolution of each SAR satellite should be selected.

3.2.2 Temperatures data preparation

Nuuk daily air temperature, SST, and LST data were utilized to examine the correlation with the variation of ice velocity. The LST temperature was generated with the geographic grid of the 7200 × 3600 Climate Modeling Grid (CMG), which represents the entire globe and composes of the 17 layers, such as temperature, observation times, and view zenith angles from the MOD11C2 products. The LST was provided in Scientific Data Sets (SDS) format and could be converted into Kelvin temperature multiplying by a scaling factor (~0.02). Then, the centigrade temperature could be calculated by subtraction of 273.15°. The SST temperature was extracted from the MODIS Global Level 3 Mapped mid-IR SST product and was provided as the centigrade temperature. We used SST4 products obtained using the mid-infrared 3.95 and 4.05 µm wavelength channels in the SST products. The LST and SST datasets were spatially resampled considering the nearest neighbor distance to the NS, because they have a lower resolution (4–5 km) and occasionally have missing values. The LST and SST temperature datasets as well as the ice velocity were interpolated with 1-day time interval to compare each other, because they have various temporal resolutions from different data resources. The LST and SST temperature data were interpolated using modified Akima spline interpolation, a type of non-smoothing spline method, and shows better fits in oscillation curves to estimate the correlation coefficients with ice velocity.

3.2.3 Uncertainty analysis

The uncertainty of glacier velocity estimated using the offset tracking method depends on the characteristics of utilized SAR data, the glacier type, and the parameters of the offset tracking (Gómez et al. 2020; Wendleder, Friedl, and Mayer 2018). Fitting errors to calculate a polynomial equation to resample the secondary images were estimated with sub-pixel accuracy in both range and azimuth directions. We evaluated the surface velocity in stable bedrock, which is expected to move more slowly, than glacier surfaces (Waechter, Copland, and Herdes 2015; Schellenberger et al. 2016; Wychen et al. 2018; Yellala, Kumar, and Arild Høgda 2019). We assumed that the surface displacement should be almost zero in the stable bedrock area and that the uncertainty could be estimated if there is some displacement even in the stable area. We selected approximately 17,500 pixels in the bedrock area near NS to estimate the uncertainty error. The estimated surface velocity in the bedrock was 0.053 m/day in the TSX/TDX dataset, 0.032 m/day in the CSK dataset, and 0.134 m/day in Sentinel-1 observations. Sentinel-1 images may have a slightly faster surface velocity because they have a larger spatial resolution (5×20 m) than that of the CSK (3×3 m) or the TSX/TDX satellites (1.2×3.3 m).

4. Results

Multi-temporal ice velocity maps and terminus positions show that the terminus of the glacier in 2020 retreated by approximately 1.1 km along the centerline of the glacier in Figure 1(c) compared to the terminus in 2013. Most of the retreat was discovered in the southern part of the glacier, whereas the northern part of the terminus line was not much changed in Figure 1(c). Relative terminus position has been gradually retreating from 2017 to 2020 with respect to the terminus in 2013. However, the glacier terminus has slightly expanded to the seafront from 2015 to 2017 as shown in Figure 1(d). To evaluate the two-dimensional spatial variation of the glacier movement over time, average ice velocity maps were created during April–June of each year (Figure 3). As the data was partially collected in the 2013 and 2020 datasets (Table S1–1), we only used the April–June dataset, which showed the highest ice movement over the course of a year. Although the velocity maps show the displacement of ice between two SAR images with different time spans, we assumed that the average ice velocity using multiple time-series of the dataset could represent the annual trend of the ice movement.

Figure 3. (a-g) Average ice velocity map of NS during April – June from 2014 to 2020. The velocity map indicated that the glacier might have a typical flow pattern reflecting ice surface elevation (around the red dot circle area in region A).

We found that the overall pattern of glacier movement was very similar over a seven-year period. The slowest ice movement could be detected in 2014 in the entire area of NS (Figure 3(a)) for a selected time window. From 2015 to 2017, an overall slight increase of 0.1–0.3 km/year in the ice velocity of the glacier was detected near the glacier terminus (red dot circles in Figure 3(c-e)). However, the ice flow was relatively slower temporally in 2018 than in previous years (Figure 3(e)). The ice velocity’s speed at the terminus in 2020 abruptly increased by approximately 0.9 km/year compared with 2014, as shown in Figure 3(g). The ice velocity maps indicate that NS had a typical flow pattern showing faster movement in red circle areas, which might reflect an ice surface elevation feature.

We plotted the mean ice velocity from 2014 to 2020 along the centerline of the glacier to investigate its variation in time and in distance from the terminus (Figure 4). Although the overall velocity pattern trends for seven years were pretty similar, the ice speed accelerated with the lapse of time. While the ice velocity far from the terminus in the area of B had not changed drastically, the profiles showed a sharp increase toward the terminus over 9 km (region A) with dramatic acceleration. To estimate the increment of velocity, we compared the average of the first two years (2014–2015) and the last two years (2019–2020). The overall velocities in region B were approximately 1.6 km/year with a similar flow pattern, however, the ice speed in 2019–2020 increased by approximately 10% compared with 2014–2015. The velocity at the terminus ranged from 4.0 to 5.3 km/year before 2019, whereas the ice speed in 2020 increased significantly by up to ~ 6.0 km/year. Though the ice velocity in region A, at the terminus, increased by 210–240%, the velocity in region B did not change much. Glacier surface elevation variations along the centerline in NS were extracted from the GIMP DEM at a distance of 12 km from the terminus between P1 and P6 to examine the relationship with the ice velocity in Figure 4. The elevation profile might be divided around a small turning point with elevation change (red arrow in Figure 4). Region B, where the ice speed did not change considerably, had a low slope of the glacier surface height. However, a slightly steeper slope in region A might cause the ice speed to accelerate toward the terminus. We will discuss the topographic effect in more detail in the discussion section.

Figure 4. Average ice velocity from 2014 to 2020 during April – June along the centerline which is marked in Figure 1(b). A traverse line of glacier surface height extracted from the GIMP DEM and its derivative (slope) was plotted to examine the relationship between the elevation change and the ice speed of the glacier. The ice velocity increased dramatically at approximately 9 km from the terminus.

The time series of ice velocity over six selected points are displayed in Figure 5 to examine the seasonal variation of the NS. Three LST, SST, and air temperature datasets were also plotted to examine the dependency with the ice velocity. Although seasonal variation could be detected from the ice velocity and all of the different types of temperatures, they did not show a strong correlation. We noticed that more dynamic glacier movement occurred in region A, whereas much slower glacier flows were estimated in region B. Fine temporal resolution of the SAR acquisitions promises more detailed glacier speed, which is useful to evaluate the seasonal effect in the glacier movement. Especially, the 6-day temporal resolution of Sentinel-1 from 2017 shows more elaborate profiles of the glacier speed. Multi-temporal ice velocities indicate that the glacier speed increased dramatically in early summer (June) and decelerated abruptly before the end of summer (September) overall. After deceleration, the rate returned to the pre-acceleration level regardless of the temperature. The fastest and slowest ice velocities of each year increased gradually at the P1–P3 points. Both the maximum ice velocity (~6.5 km/year) and the sharpest decline were detected in 2020. The velocity lines indicate that the ice velocity flowed faster at the glacier terminus. The maximum and minimum ice speed in each year indicates that the ice velocity of the NS has been gradually increasing. The high temporal resolution of the constructed time series could provide a detailed variation of the ice velocity. However, the seasonal variation at the P4-P6 points did not show significant changes.

Figure 5. A time series of the ice velocity from points P1-P6, air temperature, SST, and LST for seven years. P1-P3 belongs to region A, and P4-P6 belongs to region B. The dashed line indicates a trendline showing the variation at each point. The variation from 2018 to 2020 was a relatively complicated pattern due to the more acquisitions with a short revisit cycle of the Sentinel-1 SAR satellite.

We validated the constructed time series of the ice velocity map using the Program for Monitoring of the Greenland Ice Sheet (PROMICE) dataset based on the results of the Sentinel-1 SAR offset tracking (Solgaard and Kusk 2022). The overall pattern between the two observations shows good agreement; however, there are a few inconsistencies, as shown in Figure 6. We suspect that this inconsistency might be a result of the temporal baseline, window patch size, and window step size. Note that our time series of the velocity map shows more detailed information, due to the higher temporal resolution of the collected SAR observations.

Figure 6. Plots to validate constructed time series of the ice velocity with the PROMICE dataset with cross mark derived from Sentinel-1 SAR offset tracking approach. The PROMICE dataset is available from October 2016.

The STL decomposition results for each point (P1–P6) were displayed using seasonal and trend analysis to evaluate temporal variation in more detail (Figure 7). The STL method was applied assuming that the seasonal component of the ice velocity has a one-year cycle. A significant variation of the trend components could not be shown at the P4–P6 points. However, a gradual increase in the trend component from 2018 could be detected at the P1–P3 point. As the sinusoidal seasonal patterns of the P1–P3 points are irregular, relatively noisy residual components may remain, in comparison with those of the P4–P6 points. All seasonal components confirmed that the ice velocity decreased sharply in early summer (June) and then increased gradually.

Figure 7. The plots showing the STL decomposition analysis. P1-P3 belongs to region A, and P4-P6 belongs to region B. The seasonal (red), trend (blue), residual (green) component and the ice velocity (orange) were displayed. The seasonal component was determined on a one-year cycle basis.

Figure 8. (Continued).

Figure 8. (a-f) Plots between ice velocity and temperature (air, land surface, and sea surface) on P1-P6. The color schemes of the plot are red, green, blue, cyan, magenta, and black at the P1 to P6, respectively.

5. Discussion

Practical applications to track surface velocity using offset tracking have limitations with regard to spatial resolution and computation performance. Basically, the offset tracking method implies moving windows to calculate correlations as for a window patch in the entire image. Thus, heavy computation is required regardless of the window patch size. Furthermore, because a larger window patch size is often required to estimate the offset between two SAR observations, much more computation time is spent. Even though we used multi-core parallel computing power to calculate the correlation, it often took several days. To reduce the computation time, it is possible to enlarge the moving window step size for coarse offset estimation. However, the larger window step size degrades the spatial resolution of the correlation maps. The various window patches and step sizes might be useful for computing efficiency, as Baek et al. reported (Baek, Jung, and Chae 2018; Chae et al. 2017).

The spatial resolution of the offset tracking fully corresponds to the pixel spacing of the satellite images. Thus, super-high-resolution SAR images, such as spotlight mode, promise better glacier velocity maps, which can be helpful in examining local ice dynamics. Wrong estimation of the correlation peak with a large window size in a fast-moving glacier with a long temporal baseline can increase uncertainty in the velocity map. A short revisit cycle of the constellation SAR mission, which is enough to detect ice velocity, enables us to use a smaller window size. Hence, enhanced spatial resolution of the offset tracking map could be acquired.

The grounding line provides a better understanding of the glacier behavior in the terminus. A fast-moving surface that exceeds the wavelength of the SAR signal prevents the generation of the interferometric phase at a relatively long temporal baseline condition. A ground-based radar interferometer with a very high temporal baseline of a few minutes can produce a very coherent double differential interferogram to estimate the grounding line. The combination of phase change from the ground-based radar interferometer and amplitude-based change from the satellite can be attractive for investigating the glacier environment.

Multi-temporal variation of the glacier velocity maps and terminus positions show that the glacier terminus between winter and spring every year retreated continuously. Many studies have reported accelerated retreat rates in Greenland (King et al. 2018; Carr, Stokes, and Vieli 2017; Bevan, Luckman, and Murray 2012). The overall pattern of the ice velocity shows that it increased up to a rapid rate in early summer (June) every year before showing a tendency to drop down suddenly; by summer, the pace had gradually recovered. The velocity profiles along the centerline of the NS indicate that the glacier moved faster around the terminus (region A). At the same time, it flowed relatively slowly after 9 km from the terminus (region B). We suspect that the glacier surface height change could affect the ice velocity. To further evaluate the slope effect of the ice surface elevation, we calculated the derivative of the elevation, which was plotted with a solid black line in Figure 4. The slope changes gently (~0) from 9 km to 12 km in region B, whereas the slope in region A is much steeper. Note that there is a speedup zone where the large slope change can be detected 2 km, 6 km, and 9 km from the terminus. To evaluate the impact of the ice surface elevation, it is very useful to look at the elevation change over time. However, an interferometric SAR observation with a very short temporal baseline is required to maintain coherence over the moving ice surface (Hong et al. 2018). If a multi-temporal constellation mission like TDM is available over the NS, it merits further investigation.

Our results show an acceleration pattern in every late spring (May) and early summer (June), following rapid deceleration and gradual increment from fall (October–November) to spring (April – May). This glacier moving pattern might be consistent with “type 3 glacier” which is often observed at marine-terminating outlet glaciers (Davison et al. 2020; Vijay et al. 2019; Moon, Joughin, and Smith 2015; Twila et al. 2014). The seasonal component through the STL analysis indicated that the glacier velocity pattern is similar to “type 3 glaciers.” Although the STL method provided reasonable results to understand the ice flow pattern of the NS, the undesirable residual component remained at the P1-P3 points, located near the terminus. This may be related to irregular and abrupt variations of the ice velocity at the terminus.

The meltwater that flows between the base beneath the glacier and the ice could play a role in accelerating the ice velocity (Zwally et al. 2002; Parizek and Alley 2004). In contrast, it drastically suppresses the ice speed by forming a narrow channel beneath the glacier, enabling meltwater to flow faster (Sundal et al. 2011; Andrew et al. 2013; Vijay et al. 2019). The seasonal components also showed that the ice velocity increased in late spring (May) and early summer (June), then decreased rapidly before the end of summer (September). The meltwater effect could explain the abrupt variation of the ice velocity. A field survey evaluating meltwater condition might be helpful to investigate its relationship with glacier movement from very high temporal resolution maps of ice velocity.

Although both the ice velocity and the temperature represent seasonal variations, it appears that they do not show a strong correlation. To examine the correlation between them, we plotted the glacier velocity versus air, land surface, and sea surface temperatures for each month from 2014 to 2019 (Figure 8). Overall, similar plot patterns in both air and land surface temperature could be estimated. The air temperature from the plots may have a delayed impact on the sea surface temperature. Between 2014 and 2015, the glacier velocity did not change much from January – May and from September – December, whereas a relatively large variation in ice velocity was detected during June – August. Although they did not show a strong correlation, we did notice a slightly negative correlation in June and July, and a small positive correlation in August. These correlation patterns could be found in the terminus area (P1–P3), because significant variation could not be estimated inRegion B (P4–P6). In 2016, the variation in ice velocity could be noticed with a slight positive correlation in May and a negative correlation in June and July. Significant variations could not be monitored since 2017; however, very interesting plot patterns are displayed during 2018–2019. The linear pattern of the plots was displayed from late fall (October) and early spring (March) season from 2014 to 2017. This suggests that the glacier flows with almost constant speed regardless of temperature. Yet, interestingly, the plots in 2019 were scattered with relatively large variation (Table S2, Figure S1). The standard deviation of the monthly glacier velocity is relatively small during the winter season indicating a relatively stable glacier movement. However, in January and December of 2019, the standard deviation of glacier velocity moderately increased even during the winter season. We suspect that the ice velocity from 2019 could flow faster even in the winter season. Glacier dynamics can be unpredictable because they may fluctuate not only in summer but also in winter. For a more comprehensive investigation of the relationship with temperature, a time-series of high temporal glacier velocity maps with additional SAR observations could be a challenging topic for future research.

Figure 8. (Continued).

6. Conclusion

In this study, we analyzed the multi-temporal variation of the ice velocity in NS using the offset tracking method with seven years of SAR observations from July 2013 to September 2020. As the glacier got closer to the ice front, the ice velocity increased, and significant temporal variation was observed. The fastest ice speed of the selected point (P1) at the terminus of NS, reached approximately 6.5 km/year in 2020, as shown in Figure 5. We noticed that the ice velocity in 2020 increased by approximately 0.9 km/year compared with 2014. We decomposed the ice velocities into seasonal, trend, and residual components in the STL analysis. A gradual increment of the ice velocity could be shown in the overall trend component at the P1–P3 points, whereas there was no significant change at the P4–P6 points. Seasonal variation could be observed in both glacier velocity and the temperatures (air, land surface, and sea); however, we do not find positive feedback among them.

The time series of ice velocity estimated in this study has a relatively good temporal resolution to monitor the variation of the glacier. More frequent observations, which rely on revisiting the cycle of the satellite, promise much more detailed variations of the glacier movement. Despite heavy computation and coarse spatial resolution of the offset tracking method, the small SAR constellation missions, composed of several satellites, such as ICEYE, Capella X-SAR, or various near future SAR missions to be developed, could be beneficial resources for persistent monitoring and understanding of the glacier behaviors in polar regions.

Acknowledgments

The authors thank the German Space Agency (DLR), the Italian Space Agency (ASI), and the uropean Space Agency (ESA) for providing the TerraSAR-X/TanDEM-X (PID: GEO3061), COSMO-SkyMed and Sentinel-1 data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available from the corresponding author, SHH, upon reasonable request.

Supplemental data

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15481603.2023.2207301.

Additional information

Funding

This work was supported by National Research Foundation of Korea (NRF) grant funded by the Korean Government (MSIT) under Grant [NRF-2023R1A2C1003609] and was funded by the Meteorological/Earthquake See-At Technology Development Research [Grant KMI2018-02710].