TermPicks: a century of Greenland glacier terminus data for use in scientific and machine learning applications

Abstract

Marine-terminating outlet glacier terminus traces, mapped from satellite and aerial imagery, have been used extensively in understanding how outlet glaciers adjust to climate change variability over a range of timescales. Numerous studies have digitized termini manually, but this process is labor intensive, and no consistent approach exists. A lack of coordination leads to duplication of efforts, particularly for Greenland, which is a major scientific research focus. At the same time, machine learning techniques are rapidly making progress in their ability to automate accurate extraction of glacier termini, with promising developments across a number of optical and synthetic aperture radar (SAR) satellite sensors. These techniques rely on high-quality, manually digitized terminus traces to be used as training data for robust automatic traces. Here we present a database of manually digitized terminus traces for machine learning and scientific applications. These data have been collected, cleaned, assigned with appropriate metadata including image scenes, and compiled so they can be easily accessed by scientists. The TermPicks data set includes 39 060 individual terminus traces for 278 glaciers with a mean of 136 $\pm$ 190 and median of 93 of traces per glacier. Across all glaciers, 32 567 dates have been digitized, of which 4467 have traces from more than one author, and there is a duplication rate of 17 %. We find a median error of $\sim$ 100 m among manually traced termini. Most traces are obtained after 1999, when Landsat 7 was launched. We also provide an overview of an updated version of the Google Earth Engine Digitization Tool (GEEDiT), which has been developed specifically for future manual picking of the Greenland Ice Sheet.

1 Introduction

Since the 1980s, the Greenland Ice Sheet (GrIS) has been in negative mass balance due to increased surface melt and ice discharge with projected increases in sea level of 5 to 33 cm by 2100 from Greenland alone . Long-term historical trends in ice sheet mass loss show that approximately 50 % of the total mass loss since the $\sim$ 1990s is from ice dynamics alone, via fast-moving outlet glaciers that drain into the ocean . In part, this acceleration in dynamic loss may have been triggered by a warming climate (atmosphere and ocean) that induces sudden rapid retreat of outlet glacier termini . Observations of glacier retreat, however, show a high degree of heterogeneity in the magnitude, timing, and temporal patterns of this retreat across the ice sheet , which complicates our understanding of future mass change from outlet glaciers. This suggests that knowledge of past terminus change, and the potential for future terminus change, is critical for accurate forecasting of the GrIS contribution to sea level rise (e.g. ).

Glacier termini have long been an indicator of climate change, and terminus change data have been used to understand a range of processes over multiple timescales (e.g. ). In the long term ($>$ annual), terminus records are used to inform the timing of, regional patterns within, and climate controls on marine-terminating glacier retreat . Outlet glaciers can also change at sub-annual timescales, and the examination of terminus change on shorter timescales ($\sim$ seasonal) aids interpretation of the specific environmental and glaciological processes that influence glaciers . Such studies are valuable because glacier termini respond to a diverse set of mechanisms related to the geometry of the glacier–fjord system, inland ice dynamics, and the strength of climate forcing . However, determining the variables controlling seasonal variations can be difficult because changes in the climate system occur simultaneously (e.g. ). Recent work suggests that the shape of the terminus trace and how it evolves over time may provide additional information about the nature of processes dominating any given glacier . Such studies demonstrate the need for detailed tracing of the full terminus width (in map-view) at as high a temporal resolution as possible.

Numerous studies have digitized termini manually (Table ) for use in interpreting glacier dynamics in response to climate variability; however, the lack of coordination across these studies has resulted in duplicated data and heterogeneity in terms of format, quality, method, location, temporal coverage, and availability. Such factors limit the utility of terminus data to future researchers. In addition, manually picking glacier termini is a laborious process. For example, the data set from used the entire Landsat record to digitize 15 glaciers in central west Greenland, and the authors estimate that it took three undergraduate researchers nearly two summers working 15 h a week each to download imagery and digitize the full width of the terminus, or approximately 48 h per glacier. Rapidly replacing manual-picking are machine learning techniques, which have recently been developed for automated extraction of glacier termini across a number of satellite sensors (e.g. ). Manually digitized data are still needed for validation of machine learning methods and as training data. For example, methods using over 1500 training data inputs result in classification in $\sim$ 94 % of detectable images, under ideal conditions . Further, machine learning methods fail in images where ice conditions do not permit easy delineation of the terminus (e.g. mélange-choked fjords, shadowed termini), and therefore manually digitized termini will still be needed until machine learning algorithms improve. Importantly, future satellite missions imaging the polar regions are expected to continue for the foreseeable future, suggesting an ongoing need to coordinate terminus data in addition to other important glaciological observations that are highly coordinated (e.g. velocity and elevation). Here we present the most complete set of manually digitized terminus data for Greenland's outlet glaciers, re-processed for use in machine learning methods and scientific analysis. Data have been cleaned, associated with appropriate metadata where possible, and the metadata normalized so they can be easily accessed by scientists.

Table 1

Original sources for terminus traces for the TermPicks data set. Spatial coverage describes the number of glaciers and name/region(s) of the traces. Date range are the years covered by the data set. Resolution is the temporal resolution; annual is approximately one trace per year, sub-annual is more than one trace per year, decadal is approximately one trace every 10 years, and sub-decadal is more than one trace every 10 years but not each year. Method is the tracing method used by the author to digitize the terminus. The author key is the label given to that data set in the TermPicks data set.

Published source	Spatial coverage	Date range	Resolution	Method	Author key
	GrIS wide; $n=$ 47	1999–2018	Annual	Full width	PROMICE
	GrIS wide; $n=$ 14	1985–2011	Sub-annual	Full width	Bevan
	Kangerlussuaq; $n=$ 1	1985–2018	Sub-annual	Full width	Bevan
	SE GrIS; $n=$ 132	1931–2010	Decadal–sub-decadal	Full width	Bjørk
	NW GrIS; $n=$ 87	1972–2021	Annual	Box	Black
	Kangerlussuaq; $n=$ 1	2013–2018	Sub-annual	Box	Brough
	NW and SE; $n=$ 276	2000–2015	Annual	Box	Bunce
	NW GrIS; $n=$ 10	1976–2012	Decadal to monthly	Box	Carr
	GrIS Wide; $n=$ 273	1992–2010	Decadal	Box	Carr
	Humboldt; $n=$ 1	1975–2012	Decadal–sub-decadal	Full width	Carr
	CW GrIS; $n=$ 15	1965–2018	Sub-annual	Full width	Catania
	GrIS wide; $n=$ 65	1972–2019	Sub-annual	Full width	Cheng
	E GrIS; $n=$ 10	1993–2012	Sub-annual	Box	Sole
	GrIS wide; $n=$ 224	1984–2017	Annual	Full width	Fahrner
	N GrIS; $n=$ 21	1916–2015	Annual	Box	Hill
	N GrIS; $n=$ 18	1948–2015	Annual	Box	Hill
	GrIS wide; $n=$ 452	1978–1987	Annual	Full width	Korsgaard
	GrIS wide; $n=$ 203	1992–2007	Sub-decadal	Box	Moon
	GrIS wide; $n=$ 199	2000–2010	Sub-annual	Full width	Murray
	GrIS wide; $n=$ 28	1990–2016	Sub-annual	Full width	ESA
TermPicks	E and W GrIS; $n=$ 13	1985–2019	Sub-annual	Full width	TermPicks
	GrIS wide; $n=$ 226	1992–2017	Annual	Full width	Wood
	Helheim, Jakob.,	2009–2015	Sub-annual	Full width	Zhang
	and Kanger.; $n=$ 3

2 Methods

2.1 Input data

Terminus traces were collected through email requests to authors who had published papers that made use of such data, or they were taken from publicly available online databases (Table ). Since there was no open call for data submission, there may be other sources of terminus trace data that are available and/or unpublished. Authors used a range of image sources (Table ), but the bulk ($\sim$ 70 %) of terminus traces originate from Landsat images. Collectively, we refer to these collected data as input data to differentiate these data from the output (cleaned, reformatted) training data generated.

All data were provided in ESRI shapefile format (Fig. ) with the bulk of data provided as polylines and a smaller volume of data provided as polygons or polygon boxes. In these latter cases, the polygons were cropped at the terminus and converted into polylines. All glacier terminus traces were exported into a single ESRI line shapefile format consistent with file formats typically used in machine learning techniques. All shapefiles were re-projected into NSIDC (National Snow and Ice Data Center) Sea Ice Polar Stereographic North (EPSG:3413).

Figure 1

Flow chart showing processing pipeline for producing consistent terminus trace training data.

[Figure omitted. See PDF]

Glacier termini were commonly traced by importing geographically rectified images into GIS software (e.g. ArcGIS, ENVI, and QGIS) and manually digitizing the ice–ocean boundary (terminus). Authors used a range of methods for tracing termini including picking the full width or variations on the box methods. Box methods consist of using a fixed-width rectilinear or curvilinear box along the length of a fjord tracing the terminus within those bounds (for a description of these methods see ). For consistency in data format, we excluded termini that were identified with only a center point e.g. because these data do not cover the entire width of termini. Individual terminus trace files are largely indistinguishable between authors, with the exception of those who used the box method for picking the terminus since this method often produces terminus traces that are truncated before they reach the fjord wall. Across all authors, terminus traces have an average of 23 vertices per kilometer with a median of 10 vertices per kilometer.

2.2 Glacier identification

As the GrIS has several hundred marine-terminating glaciers, proper identification of glaciers is important for data management. Several prior authors have produced identification files (ID files) for GrIS glaciers including (Moon IDs) who created a glacier ID file by identifying all non-stagnant glaciers that terminate in the ocean with terminus widths of roughly 1.5 km or greater. The Moon IDs identify 239 glaciers that are assigned a numerical ID, including 6 ice cap glaciers that are marine-terminating. We received terminus traces for 278 glaciers but subsequently identified 282 glaciers by including all glaciers with a Moon and Joughin (2008) ID and additional glaciers with the following criteria: (1) surface speeds $>$ 50 m yr${}^{-\mathrm{1}}$, (2) grounding lines below sea level as determined from the BedMachineV3 bed topographic product , and (3) termini greater than or equal to 1 km in width. We excluded terminus traces where only one pick was available for the glacier over all authors, as well as land-terminating glaciers . Using this new ID system, here termed TermPicks ID, we assigned glacier IDs to each glacier in our database (Fig. ).

Our TermPicks ID file maintains consistency with the Moon IDs by including the corresponding Moon ID with the TermPicks ID file. We also include other information in the TermPicks ID file that is relevant for wide community use, including outlet glacier flux gates identified by and glacier naming schemes catalogued by in an ESRI multipoint shapefile so the data can be easily referenced with other data sets.

2.3 Data cleaning

The number of terminus traces included in an input shapefile varied across the input data. Some authors represented multiple dates per glacier within each shapefile, while others included single dates per glacier for each shapefile. Our output data merged all terminus traces for all dates together into one shapefile, and so input data were re-processed to fit into this format. Some authors included multiple glaciers per date for a shapefile, particularly when glaciers were adjacent to one another. Where possible, these shapefiles were manually split into traces representing separate glaciers, consistent with our output data format (Fig. c). This was accomplished using the MEaSUREs Greenland Ice Mapping Project (GrIMP) 2000 Image Mosaic for glaciers to be properly sorted along fjord wall boundaries or ice stream where appropriate. Traces were also clipped using the GrIMP ice mask in order to remove fjord wall traces . The mask was extended where it did not intersect earlier traces.

Figure 2

Common issues addressed in data cleaning and labeling. (a) Box method glacier traces are contained within a box that is smaller than the full terminus width at Glacier 224, (b) Landsat 7 ETM+ Scan Line Corrector (SLC)-off image line artifacts at Glacier 291, and (c) a single shapefile containing several different glaciers (IDs 27–30) that need to be split manually into separate glaciers to be consistent with the ID scheme. Additionally, all three images show varied levels of obstruction of the terminus in the fjord due to ice mélange. Landsat-7 and Landsat-8 images courtesy of the U.S. Geological Survey.

[Figure omitted. See PDF]

Traces that were digitized using the box methods were not interpolated to the fjord wall. In many cases, the box spans nearly the entire width of the fjord, but several data sets use boxes that are much smaller than the width of the fjord (Fig. a). The lack of data at the edges of glacier termini may lead to differences in total retreat using these data compared to other data . Thus, terminus traces digitized using the box method are flagged in the metadata (Table ).

Table 2

Image sources used in this compilation of manually traced glacier terminus trace data set. SAR signifies synthetic aperture radar.

Source	Start	End	Spatial	Temporal	Sensor
name	date	date	res. (m)	res. (days)	type
ASTER	Jan 2000	Nov 2020	15–19	16	Multispectral
Landsat 1	Jul 1972	Jan 1978	80	18	Multispectral
Landsat 2	Jan 1975	Aug 1983	80	18	Multispectral
Landsat 3	Mar 1978	Sep 1983	80	18	Multispectral
Landsat 4	Jul 1982	Dec 1993	30	16	Multispectral
Landsat 5	Mar 1984	Jan 2013	30	16	Multispectral
Landsat 7	Apr 1999	Ongoing	30	16	Multispectral
Landsat 8	Feb 2013	Ongoing	30	16	Multispectral
Sentinel 1	Apr 2014	Ongoing	20	6–12	SAR
Sentinel 2	Jun 2015	Ongoing	10	12	Multispectral
SPOT-1	Feb 1986	Dec 1990	20	26	Multispectral
Corona	Jul 1959	May 1972	7.5	Irregular	Photograph
Seventh Thule Expedition	1933	1933		Single	Photograph
aerial oblique photos
British Arctic Air Route	1931	1931		Single	Photograph
Expedition (BAARE)
Danish aerial photos	1978	1987		Single	Photograph
US Navy/US Army Air Forces	1943	1943		Single	Photograph
ALOS-PALSAR	Jan 2006	Apr 2011	10–20	14	SAR
ENVISAT	Mar 2002	Apr 2012	30	35	SAR
ERS-1	Jul 1991	Mar 2000	30	3, 35, and 168	SAR
ERS-2	Apr 1995	Sep 2011	30	3, 35, and 168	SAR
JERS-1/ Fuyo-1	Feb 1992	Oct 1998	18	44	SAR
TerraSAR-X	Jan 2008	Dec 2020	40	11	SAR
RADARSAT 1	Nov 1995	Mar 2013	100	11	SAR

Table 3

Flags assigned to output terminus trace data, created in conjunction with CALFIN . All data in the TermPicks data set have the prefix of X = 0.

Flag code	Issue
X $=$ 0	Manually digitized trace
X $=$ 1	Machine-generated trace
X0	No issues
X1	Trace uncertainty due to environment or image
	issues (clouds, shadows, missing data, etc.)
X2	Supplemented trace
X3	Landsat 7 SLC-off
X4	Incomplete/box method
X5	Automatically assigned image ID

2.4 Metadata creation

Consistent and uniform metadata are critical for the use of training data in machine learning and scientific studies. Feature extraction using image segmentation techniques rely on accurate attribution of training data to the correct time, location, and satellite image used for terminus tracing. Input data used for TermPicks suffered from a lack of consistency in the metadata, such as date format, author and satellite identification, image ID, and digitization techniques. Here we describe the metadata format for the output TermPicks data set (Fig. ). The TermPicks metadata format was chosen to be consistent with the largest archive of machine-digitized terminus traces from , known as CALFIN. For example, CALFIN includes the date, quality flags, satellite sensor, and image ID, all of which are important for machine learning. Figure  shows examples of the metadata structure for the data.

• Date columns. The date column represents the acquisition time for the image used to digitize the terminus for that trace. There are four additional columns for year, month, day, and decimal date. The date column is a string, and the format is ”YYYY-MM-DD”. Year, month, and day are integers. If a trace included only year information, the date column format is “YYYY-00-00”.

• Satellite. Satellite refers to the original sensor or satellite that produced an image used to digitize the terminus. This information was taken from existing attribute tables or file names from the input data and was used to determine the image ID where possible. The names used are listed in Table .

• Author. All people contributing traces have been listed as authors in this paper. Included in the metadata is the author identifier connected to a specific citation using the data provided. We also provide a code block in the code repository to produce citations for the authors of terminus traces that are used in data downloads. This allows for proper attribution to the correct author depending on the location and time span of data downloaded. In the data set, the author “TermPicks” refers to terminus traces produced with TermPicks GEEDiT (Google Earth Engine Digitization Tool) but not published elsewhere (Supplement C).

• Image ID. Image ID refers to the image scene identifiers for the original image used to digitize the individual glacier trace. This corresponds directly to the sensor. For example, a Landsat product ID is an example of an image ID. Certain images (e.g. some aerial images) were used to digitize multiple traces. The image ID includes information on the date and location for the original image. This may be listed as a file name that the original author used and may be stored locally (Fig. ; Glacier 291) or an image ID from a different satellite (e.g. Sentinel-1 product folder name). If an author included an image ID, the text was kept the same in case users need to contact the original author for image access.

• Glacier IDs. The glacier ID refers to the TermPicks glacier ID scheme that was created for this project (described in Sect. 2.2).

• Center X and Y. A centroid point was created for each trace in WGS 84 (EPSG:4326) so that the TermPicks data can be easily referenced with other data sets.

• Quality flag. Quality flagging is used to identify and classify traces that may have issues leading to sources of error. This quality flagging scheme was created in conjunction with to enable data synthesis between our data and machine-generated terminus traces. We assign a prefix “X” for all data defining whether the trace was created automatically or manually, with X $=$ 0 for TermPicks data and X $=$ 1 for CALFIN data, or any machine-generated terminus traces that may be included in the future. In addition, traces can have multiple quality flags. We follow the quality flag scheme in Table . In this scheme, flags are assigned if there are no issues with the terminus trace (X0), if there is uncertainty in the trace due to environmental or image issues, for example clouds partially obscuring the terminus (X1), if the trace was supplemented (two images were used to digitize the terminus) (X2), if the trace was digitized with the Landsat 7 sensor when the Scan Line Corrector (SLC) was off (X3), if the trace was digitized using the box method and is thus incomplete (X4), and if the image ID was automatically assigned because of lack of information provided in the input metadata (X5). The X1 and X2 flags are only used if the trace author indicated this information, and so many traces will not include these flags. If there are multiple flags, they are separated by commas (Fig. ; Glacier 278).

Figure 3

Example metadata for the TermPicks data set. Each column corresponds to the description in Sect. 2.4 (“Metadata creation”).

[Figure omitted. See PDF]

2.5 Landsat image scene identifiers

Satellite image scene identifiers (image IDs) are useful to find the original image from which a glacier terminus was digitized, which is a requirement for these data to be useful for machine learning. Including image IDs is also useful in cases where scientists want to explore other features in the scene at the time of a terminus trace (e.g. iceberg distribution, sediment plume occurrence). These were provided in very few of the input data sets. Where no image ID was available, Landsat scene identification is assigned to terminus traces that were originally digitized using Landsat data. Scenes were assigned by geolocating a path–row from the Worldwide Reference System (WRS-1 for Landsat 1-3; WRS-2 for Landsat 4 onward) that is closest to the terminus trace and then searching by date using Google Cloud Services. As Landsat scenes are freely available for level-1 data on Google Cloud Services and most ($\sim$ 70 %) of the data are derived from Landsat images, only terminus traces that were known to be digitized with Landsat data are assigned IDs (Fig. ). Some glaciers share multiple overlapping Landsat path–row combinations resulting in some terminus traces having two scenes assigned. In these cases, both image IDs are appended to the metadata. Glaciers with automatically assigned image IDs have the quality flag of 05 (Fig. ; Glacier 3). Further, some terminus traces did not have dates that corresponded to an image ID from Google Cloud Services and were not assigned an image ID.

2.6 Calculation of terminus change and variability

In addition to providing manually digitized terminus traces for glaciers in Greenland, we also computed terminus position change. As many previous studies have already published on terminus change over time, we provide these estimates largely as a check on our data set. We compute terminus position in two ways. First, we calculate terminus position using a method developed in in which equally spaced points along each terminus trace are projected to the nearest location along the glacier centerline. The average position of all projected points on the centerline thus becomes the average position of the glacier terminus for that date of the terminus trace. We call this the “interpolation method”. The interpolation method is most accurate when the glacier traces are all approximately the same length (i.e. not a mixture of full-width and box-method termini). Second, we calculate the fluctuation in terminus position simply by taking the point where the terminus intersects the centerline of each glacier following , here named the “centerline method”. Traces that were missing day and month information were assumed to have a timing of mid-year. Retreat rates were then calculated by taking the distance between each of these terminus positions over time. We use centerlines from where available for the glaciers in our database. Remaining centerlines were manually mapped from the MEaSUREs GrIMP 2000 Image Mosaic through the center of the glacier and the terminus traces.

We also computed the terminus seasonality as a measure of the total variation in the terminus position over the annual cycle. This is quantified using the standard deviation of the difference between raw terminus position data and smoothed terminus position data from the centerline following . We estimated seasonality for glaciers in years when there are terminus traces in at least three unique months.

Finally, we calculated the terminus sinuosity as a way to characterize the shape of the terminus as the sinuosity quantifies how much the terminus deviates from a straight line. Sinuosity is classically used in river morphology to describe map-view morphological changes in river channel patterns and is the ratio of along-channel length to valley length . Here, terminus sinuosity is measured as the length of the terminus divided by the straight line distance between the terminus endpoints. The sinuosity of rivers depends on river valley geology with typical values between 1 and 3 ; however, we do not expect glacier termini to exceed a sinuosity of 2 (i.e. the terminus will be less than twice the length of the distance across the fjord) because calving will likely occur for the parts of the terminus that are extremely anomalous. Increased sinuosity of glacier termini may be associated with crenulated terminus morphology that is thought to result from localized terminus melt as a result of buoyancy-driven plumes ; however, a smooth but highly concave terminus may also have a high sinuosity. Low-sinuosity termini may be associated with glaciers that calve via full-thickness calving events, causing fjord-width step changes in the terminus position with each calving event . While additional metrics of the geometry (e.g., curvature) may be necessary to completely describe the morphology of glacier termini, the change in sinuosity in time may reveal differences in processes affecting a single glacier.

2.7 Error estimation

Terminus traces from different authors on the same date do not necessarily align with each other, and so we quantified the difference between these traces. As a metric of error between data sets, we calculated the Hausdorff distance (commonly used in pattern recognition), the greatest minimum distance between two lines . A larger Hausdorff distance indicates two lines are less similar to each other; however, large Hausdorff distances could also indicate that two otherwise identical lines have different endpoints (different lengths). To avoid this latter issue, we trimmed each terminus trace to a glacier reference box, modified from those used by , before computing Hausdorff distances. We also excluded traces that did not span the width of these glacier boxes. Excluding short traces reduced the data set to 25 355 (65 % of the original TermPicks data set). Then, we calculated the Hausdorff distance between every pair of traces for traces that were digitized at the same glacier and on the same date by multiple authors. We identify 2671 individual instances where multiple authors digitized a glacier on the same date (sometimes more than two authors). This resulted in a total of 5748 duplicated traces.

3 Results

The TermPicks data set includes 39 060 individual terminus traces for 278 glaciers with a mean and median number of traces per glacier of 136 $\pm$ 190 and 93, respectively. However, the trace count varies depending on an author's interest in a specific glacier or region of glaciers (Fig. ). Across all glaciers, 32 567 dates have been digitized, of which 4467 have traces from more than one author. This represents duplicated efforts of $\sim$ 17 % of the input data. Traces extend back to 1916 for a small number of glaciers, but the greatest number of traces were obtained between 2000 and 2017 (Fig. ). See the Supplement for information on individual glacier coverage and statistics (Figs. S9, S10, S11), as well as access to a kmz file that can be viewed in Google Earth that produces a quick look at location and coverage for each glacier.

Figure 4

(a) Terminus trace count for glaciers in Greenland. Each circle is centered on a location of a glacier in the TermPicks ID file. The size of the circle reflects the total number of terminus traces available for that glacier. (b) The same data organized by drainage basin. Circle size reflects the total number of traces for that basin. The numbers inside of or adjacent to the circle represent the number of individual glaciers in each basin with terminus traces. Each basin is defined by the ESA/NASA ice sheet mass balance inter-comparison exercise 2016 (IMBIE; ) which includes basins from and . They are labeled by their geographic location. Region labels are NO $=$ north, NE $=$ northeast, SE $=$ southeast, SW $=$ southwest, CW $=$ center west, and NW $=$ northwest.

[Figure omitted. See PDF]

3.1 Terminus change and variability

The retreat time series using the interpolation method reveals small errors that are present as anomalous spikes in the retreat record possibly due to traces that have different endpoints (e.g., Fig. ). Centerline retreat as an average over each decade of the observational record (1940–2010 when sufficient data permit) shows regional patterns of retreat before 1990 and more ubiquitous retreat after 1990 (Fig. ). Glacier terminus seasonality varies over time and space. Out of the 19 authors in our data set, 10 are able to resolve a seasonal signal for at least one glacier for at least 1 year (Fig. ). The Catania data are able to resolve seasonal signals across the longest time period (1985–2019); however, this is only for 15 glaciers. The Murray data set resolves seasonality for 199 glaciers but only between 2000 and 2009. In contrast, the TermPicks data set resolves seasonality for the most glaciers ($n=$ 221) at different levels of completeness over the longest period of time (1985–2019). For example, Glacier 116 has traces from seven authors (Fig. ), allowing us to examine changes in seasonality from over $\sim$ 35 years between 1986 and 2017. In contrast, the data from Murray only resolve seasonality for Glacier 116 for 8 years between 2000 and 2008. Finally, we find increases in the amplitude of terminus seasonality during periods of terminus retreat for all three of our example glaciers (Fig. ).

Figure 5

Heatmap of glacier traces in each regional basin from ESA/NASA ice sheet mass balance inter-comparison exercise 2016 (IMBIE; ) in this study. Total number of traces per region can be found in Fig. . The $x$ axis is year, and the $y$ axis is the basin ID. The color corresponds to the number of traces for that basin's glacier per year; 0 traces are grey.

[Figure omitted. See PDF]

Figure 6

Terminus positions for Glacier 152 (Kangerlussuaq Gletsjer) from 11 August 2006 for three authors. Bunce (pink) and Cheng (blue) traces end before the northern fjord wall, while the ESA (yellow) trace ends at the northern wall. The table shows each calculated retreat amount since the 1978 position using the interpolation method and the centerline method. Landsat-7 image courtesy of the U.S. Geological Survey. Base image has reduced saturation to increase contrast with traces.

[Figure omitted. See PDF]

Figure 7

Decadal retreat patterns for available TermPicks data using the centerline method. For each panel, the entire decade of traces were averaged to produce an average position for that decade. The 1940/1950s are an average over both decades as there are fewer traces available in the 1950s. Then the average position for the decade is differenced from the average position of the previous decade. The size correlates to the magnitude of terminus change, while red (negative) indicates retreat and blue (positive) indicates advance.

[Figure omitted. See PDF]

Figure 8

Locations of glaciers that include terminus delineations for at least three unique months, which is the minimum number of traces required to resolve seasonality, for the entire TermPicks data set and a subset of authors. The size of the blue circle indicates how many years that there are enough traces to resolve seasonality, ranging between a single year to up to 40 years.

[Figure omitted. See PDF]

Figure 9

Example seasonality plots for three glaciers, F. Graae Gletscher (116), Heinkel Gletscher (109), and Humboldt Gletsjer (85). The location of each of these glaciers is noted in Fig. . Each color corresponds to either the entire TermPicks data set (purple) or an individual author. Glacier 85 has no individual author data set that can resolve the seasonality.

[Figure omitted. See PDF]

We calculate the sinuosity of Kangerdlugssup Sermerssua (Glacier 291) and Sermeq Silarleq (Glacier 288) between 1990 and 2020 as there is the highest density of traces after 1990 (Fig. ). Terminus sinuosity is found to vary generally between values of 1 (straight across) and 2 (highly sinuous). We examine two examples with different retreat histories. Glacier 291 is a stable glacier over the observational time period and has a similarly stable sinuosity with a mean of 1.43 $\pm$ 0.12 between 1990 and 2020 (Fig. ). In contrast, Glacier 288 undergoes a two-stage retreat beginning in 1998 with a slower-paced stage of retreat until $\sim$ 2010 when retreat accelerated through to today. This glacier has a mean sinuosity of 1.35 $\pm$ 0.17; however, we observe that the slower period of retreat is tied to a period of increased terminus sinuosity of 1.41 $\pm$ 0.18 (Fig. ), while the period of more rapid retreat experiences a decrease in terminus sinuosity to values of 1.29 $\pm$ 0.13 (Fig. ).

Figure 10

(a) Terminus change between 1990 and 2020 colored by sinuosity for Glacier 288 (Sermeq Silarleq). The dashed grey line is the start of progressive retreat as defined in . (b) Corresponding map-view terminus traces for Glacier 288 with every fifth trace colored by sinuosity. (c) Terminus change between 1990 and 2020 colored by sinuosity for Glacier 291. (d) Corresponding map-view terminus traces for Glacier 291 (Kangerdlugssup Sermerssua) with every fifth trace colored by sinuosity. The base map in (b) and (d) is the bed from BedMachine . The black pixels in (b) are errors; however, they do not impact the overall interpretation of the bed. The bed scale bar applies to both (b) and (d). The white arrows indicate glacier flow direction. The red star in the inset map is the location of the glacier on the Greenland Ice Sheet.

[Figure omitted. See PDF]

3.2 Spatial and temporal bias

Heatmaps of the output data demonstrate the temporal coverage and frequency of the data. We present heatmaps for both regional groups of glaciers (Fig. ) and individually for each glacier (Figs. S9, S10, S11). These figures demonstrate that terminus data availability is intimately tied to Landsat image acquisition. A combination of US-centric acquisition strategies, ground station coverage, and limitations on data transmission and duty cycles meant that much of the world did not have regular repeat Landsat coverage until 2013 with the launch of Landsat 8, which follows a continental acquisition strategy . Further, the failure of Landsat 6 upon launch in October of 1993 meant that imagery was only obtained in a limited capacity (via extension of the Landsat 5 satellite) until the successful launch of Landsat 7 in 1999, when we observe an increase in terminus trace data (Fig. ). We further compute the percentage of terminus traces for a given glacier compared to all available Landsat images that cover any particular glacier (see Fig.  for four examples) in order to examine the completeness of the terminus data for all glaciers. All glaciers have an individual coverage figure that is contained in our Google Earth file (Supplement). From this analysis we find that Sermeq Silarleq (ID 288) has traces from 33.1 % of all available Landsat images (including cloudy images), the most of any glacier in our data set. However, on average only 5.8 % of available Landsat images have been manually traced per glacier.

Figure 11

Examples of Landsat image availability (gray) versus termini traced (orange) for (a) Kangilliup Sermia (Rink Isbræ; 1), a relatively well-traced glacier, (b) Qeqertaarsuusarsuup Sermia (Tracy Gletscher; 73), a glacier representative of the average number of total traces for this data set, (c) Sermeq Silarleq (288), the glacier with the highest percentage of available Landsat images that have been traced in this data set, and (d) an unnamed glacier (251), representative of the average percentage of available Landsat images that have been traced in this data set.

[Figure omitted. See PDF]

Figure 12

Median error between pairs of authors, for instances where those authors have duplicated a glacier trace on a given date. No color indicates two authors have no duplicated traces between them.

[Figure omitted. See PDF]

Regional differences in data availability also exist (Figs.  and ). Higher-latitude glaciers experience more frequent coverage by satellite image sensors than lower-latitude glaciers due to increased scene overlap at high latitudes (e.g. Fig. b after 2013). In southwest Greenland, there are fewer traces simply due to the lack of marine-terminating glaciers in this region, which is primarily drained through land-terminating ice. There are also fewer overall traces in north and northeast Greenland than central west Greenland, a region with a similar number of glaciers, potentially due to less interest in tracing in north and northeast Greenland (Fig. ). The densest coverage is in central west and northwest Greenland (IDs 279 to 3) where nearly every available image from Landsat and other sensors was traced to create as complete a record as possible of regional glacier change. Other glaciers of interest include Helheim, Kangerlussuaq, and Sermeq Kujalleq (Jakobshavn; IDs 181, 152, and 278), which also have dense coverage.

3.3 Error in manual digitization

The overall median error between pairs in this reduced data set is 107 m, which is comparable to that obtained in most machine learning studies when comparing machine-traced termini to manually traced termini . The median error between any given pair of authors varies with the greatest median error (7350 m) between Cheng and Hill, and the lowest median error (58.6 m) between Fahrner and TermPicks (Fig. ). The magnitude of errors are not necessarily due to inaccurate digitization by authors but can be explained by Hill and other authors focusing on northern glaciers (which can be difficult to trace due to the presence of near-terminus crevasses), as well as Fahrner focusing on late summer observations when the glacier margin is often most clear. The mean and median of the median errors for each author are presented in Table , and there was no clear distinction in error based on methodology used (box versus full-width tracing). Traces with $>$ 500 m error between traces were manually checked for errors (220 traces). If two traces were on the same date but the trace was not equivalent (e.g. the trace did not appear to be from the same front), then the trace with more complete metadata (e.g. includes the original image ID) was kept. If a trace had three authors and one was not equivalent, it was removed. Only 0.4 % of total traces were removed from the data set through this manual checking. In some cases, there are glaciers that have higher errors than other glaciers (e.g. IDs 39, 73, 86, 99, 100, and 101) due to the fact that they appear to have highly fractured ice tongues, and they develop long, linear cracks that authors may or may not trace in their entirety.

Table 4

Mean vertices per kilometer of trace, as well as mean and median of the median errors of each author compared to other authors. n/a – not applicable

Author	Vertices	Mean median	Median median
	per km	error (m)	error (m)
Bevan	2.5	227.5	145.8
Bjørk	14.2	113.6	113.6
Black	5.7	181.9	111.2
Brough	n/a	n/a	n/a
Bunce	14.1	109.0	88.3
Carr	7.1	201.0	98.0
Catania	18.3	112.7	100.9
Cheng	211.1	720.5	171.8
ESA	10.4	321.9	317.8
Fahrner	5.9	139.3	122.5
Hill	10.0	1458.8	309.1
Korsgaard	9.7	n/a	n/a
Moon	5.5	148.0	148.0
Murray	6.3	106.7	96.5
PROMICE	16.5	355.5	133.2
Sole	5.4	228.1	144.5
TermPicks	11.8	113.7	78.7
Wood	23.1	114.5	96.7
Zhang	55.7	421.8	452.0

Termini traced with different methods or widths of the glacier may have some systemic differences in terminus retreat over time . For example, Fig.  shows Glacier 152 (Kangerlussuaq Gletsjer) on 11 August 2006. This date was digitized by three separate authors (Bunce, Cheng, and ESA) at different extents of the glacier front. When the interpolation method is used, there is a 0.5 km difference in terminus position change because the endpoints for each trace are different. Bunce and Cheng will show a higher retreat compared to ESA because the interpolation method accounts for the entire width of the glacier. Therefore the mean positions of the Bunce and Cheng traces will be further up-glacier as they do not include the lateral tails seen in the ESA trace. While there is no large-scale difference between retreats calculated from the box method versus full-width traces, users of these data should be aware of this potential misfit between traces based on endpoints. For example, Bunce traces use the box method, while Cheng traces use the full-width method; however, they both end before the fjord wall. Glacier 152 has dead ice on its northern margin, and, as shown in the image, the scan line errors in the Landsat 7 imagery block some of the ice, and so some authors may or may not digitize the entire front for numerous reasons.

4 Discussion

This is the first published study of manually traced Greenland-specific marine-terminating glacier traces with consistent metadata and formatting across multiple data sets from different authors. Glacier terminus traces have been a staple indicator of glacier change for decades (e.g. ). From this paper alone, 22 sources have digitized and interpreted terminus positions in Greenland, with many more using these data to aid interpretation of GrIS change. However, all of these efforts have happened independently, with duplicate efforts and lack of consistency across data format and accessibility. For example, Fig.  shows a time series of Glacier 116 (F. Graae Gletscher) with author labels for each trace. This figure demonstrates the utility of combining data sources, which enables a more complete view of terminus change at this glacier than any previously published individual study. We find similar ice-sheet-wide retreat patterns as previously published sources. For example, total retreat for 2000–2010 is $\sim$ 252 km in 225 glaciers (Fig. ), which is comparable to who found $\sim$ 267 km in 199 glaciers. We find the greatest retreats occur from 1990 to 2010 (Fig. ), similar to and . Finally, we find a rapid increase in retreat beginning in the 1990s–2000s (Fig. ), similar to , , and . While we recognize that not every glacier has a complete time series or the ability to resolve seasonal changes in terminus position over all years and that there remain limitations in drawing large-scale conclusions on retreat patterns with these data alone, we find increases in the amplitude of terminus seasonality during periods of terminus retreat (Fig. ). This may be related to the changes in fjord geometry that glaciers experience as the terminus retreats through overdeepenings.

Figure 13

Example terminus change for Glacier 116 (F. Graae Gletscher). Color and symbol correspond to different authors for each pick.

[Figure omitted. See PDF]

An additional value of the TermPicks data set is that it provides map-view trace data, not just centerline data, thus informing on morphological changes to the terminus over time. We explore the value of this through the examination of the terminus sinuosity, but other measures (e.g., terminus curvature) may also be valuable in contextualizing terminus morphology. While the mean sinuosities for Glaciers 288 and 291 (Fig. ) are similar, we find variations in sinuosity for the glacier that experienced large-scale retreat (Glacier 288) compared to the one that has remained stable over the observational period (Glacier 291). Glacier 291 is known to have a terminus that is dominated by plume-driven melting , and so we might anticipate increased sinuosity related to local melting associated with these plumes . In contrast to this, the terminus of Glacier 288 begins with a relatively low sinuosity and then during the period of slow retreat (1998–2010) experiences an increase in sinuosity (Fig. ), suggesting that this glacier may also have experienced enhanced terminus melting due to subglacial discharge plumes during this time. Subsequently, Glacier 288 experiences a period of more rapid retreat as the glacier terminus moves into an overdeepened portion of the bed. Here, sinuosity decreases, and terminus change is dominated by full-thickness calving .

Although machine-enabled terminus tracing has made great strides in the past few years, there will be a continued need for manually tracing glacier termini. This is because certain environmental conditions, such as heavy shadows, cloud cover, ice mélange, and low solar illumination, make it difficult for current machine learning algorithms to accurately trace all available images. The data provided here will aid improvements in machine learning that will ultimately reduce the need for future manual tracing. Ideally, machine- and manual-tracing efforts would work in concert, with data gaps or large errors reported by machine learning quickly identifying where the need is the greatest for the manual-tracing team. For example, both the data presented here and the data in CALFIN are not extended beyond 2020, and there is no funding in place to provide continued coordinated (between machine- and manual-tracing authors) updates to terminus positions in the future. Coordinated effort between machine- and manual-tracing teams is warranted to ensure regular delivery of future data, given the importance to the wider scientific community.

Until fully automated, frequently updated, and publicly available terminus traces are available for Greenland and elsewhere, we anticipate that authors will continue to manually trace in studies that are spatially or temporally limited. Ideally, future efforts would occur in conjunction with our work, producing data with similar format, metadata, and visibility. To that end, we recommend the use of a bespoke version of the Google Earth Engine Digitisation Tool (GEEDiT; ) within Google Earth Engine’s (GEE) Application Programming Interface (API) . This GEEDiT–TermPicks version builds substantially on the original GEEDiT, with improvements made to the digitization interface, metadata options, sensor availability, and image accessibility. A user guide is provided in the Supplement to this paper. A major advantage of GEEDiT–TermPicks over traditional repository download and visualization approaches is that it accesses the archive of Landsat, Sentinel-1, Sentinel-2, and ASTER images on the Google Cloud servers within a standard web browser. It therefore allows for much faster access to imagery compared to the alternative of downloading, extracting, and processing each individual image. This is combined with an interface for the easy digitization of margins that now uses GEE’s DrawingTools functions to improve both speed and flexibility of digitization for users.

To ensure that future data generated using this tool will be consistent with our data set, the GEEDiT–TermPicks interface visualizes the TermPicks ID locations, allowing the user to easily identify the glaciers present and access relevant imagery. Once a glacier is chosen, GEEDiT–TermPicks provides rapid access to all available satellite images of that glacier, which can be pre-filtered by date and satellite. If the image is clear, the termini can be extracted by simply clicking on the screen along the glacier margin. Images with glacier termini that are low in quality can be compared with previous or subsequent images that are near in date to help better determine the location of the terminus for a specific date and/or time. If this is done, it will automatically be flagged in the image metadata, though this (and other) image quality flag options can be manually selected, including options to provide a written note as to why the image is inadequate. Data exported from GEEDiT–TermPicks will therefore include as standard all metadata required for easy inclusion into future TermPicks data releases.

Finally, we recommend a minimum of 11 vertices per kilometer of trace for quality, which is consistent with this database. We also recommend tracing across the entire width of the glacier terminus as previous studies have shown that information about mass loss processes can be obtained from studying the map-view change in trace morphology at high levels of detail .

5 Conclusions

We present a new compilation of outlet glacier terminus traces for the GrIS spanning a time period from 1916 to 2020 obtained through manually tracing the ice–ocean boundary. Data were cleaned, reformatted, assigned to image IDs, and quality controlled for use in machine learning algorithms that will enable semi-automated terminus tracing. Termini are provided in the same format and with similar metadata to ongoing machine-learning-based terminus tracing. We have combined TermPicks data with those from CALFIN in our data repository. We find errors in TermPicks on the order of $\sim$ 100 m, similar to machine-identified termini. We find biases in terms of data coverage with well-studied glaciers with high coverage of terminus trace data, as well as other glaciers devoid of consistent coverage, showcasing the need for further manual and machine learning efforts to provide terminus data. We provide tools for future tracing efforts and include software to enable the use of these data for the broader scientific community.

Code availability

This work was performed using freely available software, primarily Google Earth, Google Earth Explorer, Python, and QGIS. Code to generate a text file that includes the digital object identifier of citations for users is available on the GitHub site (https://github.com/sgoliber/TermPicks, last access: 20 July 2021) and Zenodo (10.5281/zenodo.6954113, ). GEEDiT–TermPicks can be accessed through Google Earth Engine Code Editor (10.5281/zenodo.6962120, ).

Data availability

Terminus trace data will be made available at NSIDC, a NASA DACC. Until the data submission is approved, data are currently available on Zenodo (10.5281/zenodo.6557981, ). A shapefile of combined CALFIN and TermPicks data is included in this repository. We also provide terminus retreat data, the TermPicks ID shapefile, and kmz file that can be viewed in Google Earth and provides a quick look at temporal coverage (compared to imagery availability) for all glaciers.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-16-3215-2022-supplement.

Author contributions

SG collected, cleaned, and formatted the data and metadata and performed most of the analyses. TB performed some data analysis and synthesis. SG, TB, and GC wrote the manuscript. JML created GEEDiT–TermPicks. DC provided expertise on machine learning. All authors contributed to the editing and refining of the manuscript and data contribution.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We acknowledge support from a NASA Earth and Space Sciences fellowship to Sophie Goliber (18-EARTH18F-323) and terminus tracers everywhere. At the University of Texas this includes Lucero Casteneda, Gene Hsu, David Peters, Andrea Jimenez, and Mason Fried. Niels J. Korsgaard was supported by the Programme for Monitoring of the Greenland Ice Sheet (PROMICE). Michael Wood was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, California Institute of Technology, administered by the Universities Space Research Association under contract with NASA. James M. Lea is supported by a UKRI Future Leaders Fellowship (grant no. MR/S017232/1). Dominik Fahrner acknowledges support for this study through the EPSRC and ESRC Centre for Doctoral Training on Quantification and Management of Risk and Uncertainty in Complex Systems Environments (grant no. EP/L015927/1). Tavi Murray is funded by the Leverhulme Trust Research Leadership scheme F/00391/J and the UK NERC NE/G010366/1. Landsat images used in figures were downloaded from USGS Earth Explorer https://earthexplorer.usgs.gov/ (last access: May 2021).

Financial support

Sophie Goliber has been supported by the NASA Earth and Space Sciences fellowship (18-EARTH18F25 323).

Review statement

This paper was edited by Kristin Poinar and reviewed by three anonymous referees.