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

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 ± 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 ∼ 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 (Mouginot et al., 2019; Enderlin et al., 2014) with projected increases in sea level of 5 to 33 cm by 2100 from Greenland alone (Aschwanden et al., 2019; Goelzer et al., 2020). Long-term historical trends in ice sheet mass loss show that approximately 50 % of the total mass loss since the ∼ 1990s is from ice dynamics alone, via fast-moving outlet glaciers that drain into the ocean (Enderlin et al., 2014; Mouginot et al., 2019; King et al., 2020). 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 (Wood et al., 2021; King et al., 2020). 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 (Moon and Joughin, 2008; Catania et al., 2018; Murray et al., 2015a; Carr et al., 2017; Fahrner et al., 2021), 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. Felikson et al., 2017; Aschwanden et al., 2019; Slater et al., 2019).

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. Warren and Glasser, 1992; Warren, 1991; McNabb and Hock, 2014; Moon et al., 2015; Cook et al., 2005; Howat et al., 2008; Howat and Eddy, 2011). 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 (Murray et al., 2015b; Catania et al., 2018; Hill et al., 2018; Bunce et al., 2018; Howat and Eddy, 2011; Wood et al., 2021; King et al., 2020; Fahrner et al., 2021; Black and Joughin, 2022). Outlet glaciers can also change at sub-annual timescales, and the examination of terminus change on shorter timescales (∼ seasonal) aids interpretation of the specific environmental and glaciological processes that influence glaciers (Fried et al., 2018; Moon et al., 2015; Schild and Hamilton, 2013; Cassotto et al., 2015; Ritchie et al., 2008; Howat et al., 2010; Carr et al., 2014; Moon et al., 2014, 2015; Brough et al., 2019; Kehrl et al., 2017; Bevan et al., 2019). 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 (Moon and Joughin, 2008; Carr et al., 2017; Catania et al., 2018; Bunce et al., 2018; Porter et al., 2018). However, determining the variables controlling seasonal variations can be difficult because changes in the climate system occur simultaneously (e.g. Cowton et al., 2018; Fahrner et al., 2021; Wood et al., 2021). 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 (Fried et al., 2018; Chauché et al., 2014). 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 1) 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 Catania et al. (2018) 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. Mohajerani et al., 2019; Baumhoer et al., 2019; Cheng et al., 2021; Zhang et al., 2021). 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 ∼ 94 % of detectable images, under ideal conditions (Cheng et al., 2021). 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.

Andersen et al. (2019)Bevan et al. (2012)Bevan et al. (2019)Bjørk et al. (2012)Black and Joughin (2022)Brough et al. (2019)Bunce et al. (2018)Carr et al. (2013)Carr et al. (2017)Carr et al. (2015)Catania et al. (2018)Cheng et al. (2021)Cowton et al. (2018)Fahrner et al. (2021)Hill et al. (2017)Hill et al. (2018)Korsgaard (2021)Moon and Joughin (2008)Murray et al. (2015a)Raup et al. (2007)Wood et al. (2021)Zhang et al. (2019)

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.

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 1). 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 2), but the bulk (∼ 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. 1) 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.

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 Lea et al., 2014). For consistency in data format, we excluded termini that were identified with only a center point (e.g. King et al., 2020) 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 and Joughin (2008) (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⁻¹, (2) grounding lines below sea level as determined from the BedMachineV3 bed topographic product (Morlighem et al., 2017b), 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 (Mouginot et al., 2019). Using this new ID system, here termed TermPicks ID, we assigned glacier IDs to each glacier in our database (Fig. 1).

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 Mankoff et al. (2019) and glacier naming schemes catalogued by Bjørk et al. (2015) 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. 2c). This was accomplished using the MEaSUREs Greenland Ice Mapping Project (GrIMP) 2000 Image Mosaic (Howat et al., 2014; Howat, 2018) 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 (Howat et al., 2014). 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.

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. 2a). The lack of data at the edges of glacier termini may lead to differences in total retreat using these data compared to other data (Lea et al., 2014). Thus, terminus traces digitized using the box method are flagged in the metadata (Table 3).

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

Table 3 Flags assigned to output terminus trace data, created in conjunction with CALFIN (Cheng et al., 2021). All data in the TermPicks data set have the prefix of X = 0.

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. 1). The TermPicks metadata format was chosen to be consistent with the largest archive of machine-digitized terminus traces from Cheng et al. (2021), 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 3 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 2.

• 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. 3; 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 Cheng et al. (2021) 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 3. 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. 3; Glacier 278).

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

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 (∼ 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. 1). 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. 3; 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 Catania et al. (2018) 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 King et al. (2020), 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 Murray et al. (2015a) where available for the glaciers in our database. Remaining centerlines were manually mapped from the MEaSUREs GrIMP 2000 Image Mosaic (Howat et al., 2014; Howat, 2018) 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 Catania et al. (2018). 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 (Schumm, 1985; Montgomery and Bierman, 2019). 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 (Schumm, 1985); 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 (Chauché et al., 2014; Fried et al., 2018); 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 (Fried et al., 2018; James et al., 2014). 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 (Huttenlocher et al., 1993). 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 Moon and Joughin (2008), 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 ± 190 and 93, respectively. However, the trace count varies depending on an author's interest in a specific glacier or region of glaciers (Fig. 4). Across all glaciers, 32 567 dates have been digitized, of which 4467 have traces from more than one author. This represents duplicated efforts of ∼ 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. 5). 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; Shepherd et al., 2012) which includes basins from Rignot and Mouginot (2012) and Rignot et al. (2011). They are labeled by their geographic location. Region labels are NO = north, NE = northeast, SE = southeast, SW = southwest, CW = center west, and NW = northwest.

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. 6). 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. 7). 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. 8). 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. 9), allowing us to examine changes in seasonality from over ∼ 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. 9).

Figure 5 Heatmap of glacier traces in each regional basin from ESA/NASA ice sheet mass balance inter-comparison exercise 2016 (IMBIE; Shepherd et al., 2012) in this study. Total number of traces per region can be found in Fig. 4. 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 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 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 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 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. 8. 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.

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. 10). 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 ± 0.12 between 1990 and 2020 (Fig. 10). In contrast, Glacier 288 undergoes a two-stage retreat beginning in 1998 with a slower-paced stage of retreat until ∼ 2010 when retreat accelerated through to today. This glacier has a mean sinuosity of 1.35 ± 0.17; however, we observe that the slower period of retreat is tied to a period of increased terminus sinuosity of 1.41 ± 0.18 (Fig. 10), while the period of more rapid retreat experiences a decrease in terminus sinuosity to values of 1.29 ± 0.13 (Fig. 10).

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 Catania et al. (2018). (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 (Morlighem et al., 2017a). 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.

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. 5) 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 (Wulder et al., 2016). 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. 5). 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. 11 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 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.

Regional differences in data availability also exist (Figs. 4 and 5). 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. 11b 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. 4). 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 (Catania et al., 2018) 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 (Cheng et al., 2021). 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. 12). 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 4, 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

Termini traced with different methods or widths of the glacier may have some systemic differences in terminus retreat over time (Lea et al., 2014). For example, Fig. 6 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. Weidick, 1958; Higgins, 1990; Warren and Glasser, 1992; Murray et al., 2015a). 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. 13 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 ∼ 252 km in 225 glaciers (Fig. 7), which is comparable to Murray et al. (2015a) who found ∼ 267 km in 199 glaciers. We find the greatest retreats occur from 1990 to 2010 (Fig. 7), similar to Wood et al. (2021) and Fahrner et al. (2021). Finally, we find a rapid increase in retreat beginning in the 1990s–2000s (Fig. 7), similar to Carr et al. (2017), King et al. (2020), and Fahrner et al. (2021). 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. 9). 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.

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. 10) 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 (Fried et al., 2015, 2018; Jackson et al., 2017), and so we might anticipate increased sinuosity related to local melting associated with these plumes (Chauché et al., 2014; Fried et al., 2018). 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. 10), 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 (Fried et al., 2018).

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 (Cheng et al., 2021) 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; Lea, 2018) within Google Earth Engine’s (GEE) Application Programming Interface (API) (Gorelick et al., 2017). 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 (Fried et al., 2018; Chauché et al., 2014).

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 (Cheng et al., 2021) in our data repository. We find errors in TermPicks on the order of ∼ 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.