Daily variations in Western Greenland slush limits, 2000–2021

Abstract

The marginal areas of the Greenland ice sheet develop streams and lakes each summer, documenting that surface runoff of meltwater is a major component of ice-sheet mass balance. Here we map the slush limit, a proxy for the extent of surface runoff, using daily MODIS data for the years 2000–2021. We develop an automated algorithm capable of detecting daily slush limits, provided sufficient image quality. The algorithm is applied to the ice sheet's western flank (61.7 $^{\circ }$N to 76.5 $^{\circ }$N). We find significant increasing trends in maximum slush limits until the year 2012, but not thereafter. We show that the slush limit typically rises quickly early in the ablation season but stabilizes before melting ceases. The data provide evidence that upward migration of surface runoff in summer 2012 stopped early at the upper margin of the ice slabs. These thick and continuous ice layers are located close to the surface, in the firn, and impede percolation of melt into deeper pore space. Had the ice slabs extended higher, the summer 2012 provided sufficient energy to raise the slush limit by another $\sim$300 m in elevation.

1. Introduction

Surface melt is being observed at increasingly higher elevations on the Greenland ice sheet (e.g. Nghiem and others, 2012; McGrath and others, 2013) and meltwater discharge has increased (Ahlstrøm and others, 2017). The observations underline that melt plays a dominant role in ice-sheet decline (Enderlin and others, 2014; van den Broeke and others, 2016). In total, 50–60% of current mass loss is attributed to increased surface melt (van den Broeke and others, 2016; Team IMBIE, 2020); the rest is mostly due to calving.

Meltwater generated at the ice-sheet surface either runs off or gets retained in snow and firn. Only meltwater that runs off contributes to mass loss. Consequently, it is of scientific interest to quantify the extent of the runoff area, that is the area where at least a fraction of the meltwater runs off and finds its way into the oceans. The runoff area's upper margin is typically labeled runoff limit (e.g. Pfeffer and others, 1991; Reeh, 1991; Braithwaite and others, 1994). The runoff limit is generally located in the accumulation area (Shumskii, 1955, 1964). The same author shows that runoff takes place where either (i) the refreezing capacity of the firn is smaller than the latent heat contained in melt and rain, or (ii) the sum of melt and rain exceeds a certain fraction of firn pore space.

Downslope from the runoff limit follows a zone where runoff is already substantial but its annual sum is clearly smaller than total accumulation. This zone has received limited attention in the literature and lacks established terms. The closest equivalents might be the slush limit (Müller, 1962), the slush line (Greuell and Knap, 2000) or the visible runoff limit (Tedstone and Machguth, 2022). All of these terms refer to water or runoff becoming visible at the surface. However, visible surface discharge develops only in accumulation areas where the refreezing capacity of the firn is not exhausted but there is insufficient pore space to accommodate melt and liquid precipitation (Shumskii, 1955). If pore space suffices but not refreezing capacity, runoff will take place inside the firn (Shumskii, 1955), with little or no visible expression on the surface. Both modes of runoff involve an aquitard that guides lateral motion of water: runoff inside the firn follows high firn densities near the pore close-off depth (tens of meters below the surface; Miller and others, 2018, 2020), while surface runoff is guided by low-permeability icy firn or thick near-surface ice layers (lying on top of otherwise porous firn; Machguth and others, 2016; Leone and others, 2020), termed ice slabs (MacFerrin and others, 2019).

On Greenland, the geographical extents of both modes of meltwater runoff depend on annual amounts of snow fall (MacFerrin and others, 2019) and have limited overlap (Tedstone and Machguth, 2022). Runoff inside the firn prevails where firn aquifers are present (Forster and others, 2014; Koenig and others, 2014). Aquifers have been mapped using Operation IceBridge radar (Miège and others, 2016).

Surface discharge on the ice sheet has been studied already by the US military (ACFEL, 1947; Holmes, 1955) and in the context of hydropower projects (Thomsen and others, 1989). Recent remote-sensing applications have focused on supraglacial lakes (e.g. Box and Ski, 2007; Miles and others, 2017) and mapping of surface streams (e.g. Yang and others, 2015, 2019).

Few studies have explicitly (Holmes, 1955) or implicitly (Poinar and others, 2015; Yang and others, 2019) addressed the runoff limit. Tedstone and Machguth (2022) investigated the visible runoff limit, over the entire ice sheet and most of the satellite era. They define the visible runoff limit as the highest location where supraglacial drainage networks are seen in optical satellite imagery. Using Landsat, they show that Greenland's visible runoff area has increased by 29$^{ + 6}_{-8}$ % since 1985. When compared to Miège and others (2016), the study also shows that surface runoff exceeds aquifer discharge by area.

The Landsat satellites provide one of the oldest remote-sensing archives suited to map visible runoff limits. High spatial resolution allows recognizing individual runoff features. The low temporal resolution, however, limits insight into annual evolution of visible runoff extent. Better suited to obtain such information are the Advanced Very High Resolution Radiometer (AVHRR) or the Moderate Resolution Imaging Spectrometer (MODIS). Both sensors provide daily global coverage, at the cost of relatively low spatial resolution.

Greuell and Knap (2000) demonstrated that AVHRR imagery (1.1 km resolution at nadir) is suitable for detecting the seasonal evolution of the ‘slush line’, which they consider a proxy for abundant surface runoff. The slush line refers to the uppermost occurrence of slush, which is snow with all pore space water-filled (Cogley and others, 2011). Holmes (1955) distinguishes between white slush (slush covered by unsaturated snow) and blue slush, that is a fully water-saturated snowpack which appears blue at the surface. In line with Greuell and Knap (2000) we here refer to slush as ‘blue slush’ which is easily recognizable by eye and from space.

This study builds on the work by Greuell and Knap (2000) by using MODIS to chart the evolution of the daily slush limit over the melt seasons. Instead of slush line we use the term slush limit $\Upsilon _{\rm S}$, coined by Müller (1962) to denote the uppermost visible appearance of ‘material that is water saturated’ (Benson, 1996). We develop an automated algorithm able to reliably detect the slush limit $\Upsilon _{\rm S}$ at daily resolution, for most of Greenland's west coast from 2000 to 2021. We use the output of the algorithm to evaluate whether the rise of $\Upsilon _{\rm S}$ follows identifiable patterns or differs strongly between melt seasons or regions. Finally, we quantify how annual maximum slush limits $\max \Upsilon _{{\rm S}}$ have evolved over time.

2. Study area

We focus on the western flank of the Greenland ice sheet where surface runoff is common (e.g. Holmes, 1955; Thomsen and others, 1989). The study domain comprises latitudes from 61.7 $^{\circ }$N to 76.5 $^{\circ }$N, or $\sim$1700 km in north-south direction. Toward its very south and north, the study area borders with regions where firn aquifers exist (Miège and others, 2016).

3. Data

The MODIS is a multispectral sensor (36 spectral bands) aboard the Aqua and Terra satellites. The satellites were launched in 1999 and 2002. The sensors provide at least daily coverage of each point on Earth at spatial resolutions ranging from 250 to 1000 m, depending on spectral band.

Here we use the data products MOD10A1 (MODIS/Terra Daily Snow Cover at 500 m resolution, version 6.0; Hall and Riggs, 2016) and MOD09GA (MODIS/Terra Surface Reflectance Daily at 500 m, version 6.0; Vermote and Wolfe, 2015). We use daily tiles from day-of-year 120 (29 or 30 April) to 280 (6 or 7 October), covering west Greenland and the entire MODIS era (2000–2021). All scenes are used, regardless of their percentage of cloud cover.

In addition to the MODIS data, we use the Arctic DEM Release7 (100 m resolution mosaic, v.3.0; Porter and others, 2018) and the outlines of the Greenland ice sheet according to Rastner and others (2012). In the analysis of derived $\Upsilon _{{\rm S}}$ we utilize all available (2009 to present) hourly records of air temperature $T_{\rm a}$ measured at the PROMICE automatic weather stations (AWS) KAN_U (1840 m a. s. l.; 67.000 $^{\circ }$N, 47.025 $^{\circ }$W) and KAN_M (1270 m a. s. l.; 67.067 $^{\circ }$N, 48.836 $^{\circ }$W) (van As and others, 2011; Fausto and van As, 2019). Furthermore, we use $T_{\rm a}$ measured at the GC-Net station Dye-2 at $\sim$2120 m a.s.l., 66.480 $^{\circ }$N, 46.279 $^{\circ }$W (Steffen and Box, 2001).

4. Methods

Our approach to detect $\Upsilon _{\rm S}$ is inspired by Greuell and Knap (2000) who used AVHRR imagery to detect slush lines during 1990–1997 for central western Greenland. Greuell and Knap (2000) use the term slush line to emphasize that the rather coarsely resolved AVHRR images (1.1 km at nadir) do not allow to directly map runoff. Their approach focuses on detecting the elevation where spatial variability of surface albedo transitions from low to high. Low spatial variability of surface albedo indicates the monotonous snow covered surface above the slush line. Higher variability of albedo is indicative of a surface where the bright snow cover is intersected with dark meltwater streams, lakes and slush fields.

Our algorithm deviates from Greuell and Knap (2000) by (i) being designed for application to MODIS, (ii) being automated and (iii) incorporating normalized difference water index for ice (NDWI$_{{\rm ice}}$; Yang and Smith, 2013) as an additional detection criterion.

4.1. Preprocessing

4.1.1. Surface albedo from MOD10A1

The snow_albedo_daily_tile data from MOD10A1 (in the following simply referred to as MOD10A1) are stitched together, then reprojected to a Polar Stereographic projection (EPSG:3413) using nearest neighbor interpolation and finally cropped to the study area.

Initially all non-albedo values and clouds in MOD10A1 are masked. Subsequently, the scenes are filtered because the data (i) are subject to residual cloud effects (such as clouds or their shadows retained as albedo values) and (ii) contain frequent artifacts (with a stripy appearance) where albedo $\alpha$ is much lower than the surrounding grid cells (cf. Box and others, 2012).

Albedo values are masked out if at least one of the following three conditions is fulfilled: (i) $\alpha < 12$%, (ii) $\alpha > 90$% and (iii) $\vert \alpha - \alpha _{11days}\vert \ge 30\percnt$. The last condition states that $\alpha$ of each individual pixel is compared to its median albedo $\alpha _{11{\rm days}}$ over an 11-day buffer (5 days before and 5 days after the scene under investigation, excluding the day under investigation). The 11-day buffer is inspired by the MOD10A1 filtering applied by Box and others (2012).

4.1.2. NDWI$_{ice}$ from MOD09GA

We calculate the Normalized difference water index for ice (NDWI$_{{\rm ice}}$; Yang and Smith, 2013) using the red (620–670 nm) and blue (459–479 nm) bands of the MODIS MOD09GA daily surface reflectance product:

(1)$${\rm NDWI}_{{\rm ice}} = {{\rm Blue} - {\rm Red}\over {\rm Blue} + {\rm Red}}.$$

Values of the resulting grids range from 1 to –1. Google Earth Engine (Google Inc) is used to process Daily NDWI$_{{\rm ice}}$ mosaics corresponding to the whole study area.

4.1.3. Digital elevation model

The Arctic DEM is adjusted to represent elevation above the EGM2008 geoid, and is downsampled to the same grid and spatial extent as the preprocessed MOD10A1 and MOD09GA data. Using the ice-sheet mask by Rastner and others (2012), all DEM values outside the ice sheet are masked.

4.2. Calculating the standard deviation of surface albedo

Greuell and Knap (2000) used the standard deviation of albedo ($\sigma _\alpha$) to distinguish areas of high and low spatial variability of albedo. They calculated spatial variability, for each day and pixel, as $\sigma _\alpha$ of a square of 5 by 5 pixels. We base our detection of $\Upsilon _{\rm S}$ on the same metric but calculate $\sigma _\alpha$ differently.

Standard deviation of albedo is calculated based on preprocessed MOD10A1. Each pixel gets attributed a $\sigma _\alpha$ which is calculated along lines of pixels to optimize computational efficiency. Standard deviation is calculated over lines of 11 pixels in length with the pixel under consideration centered. The scanning is done around each pixel in the vertical and in the horizontal, resulting in two arrays of standard deviation. The final $\sigma _\alpha$ of each pixel is calculated as the average of its vertical and horizontal standard deviation. We chose a scan window of 11 pixels because of being similar in width to 5 AVHRR pixels ($\approx$ 5.5 km at nadir). For any pixel, $\sigma _\alpha$ is only calculated when horizontally and vertically a minimum of eight out of 11 pixels contain valid $\alpha$. Otherwise the pixel value is masked. This condition automatically creates a small buffer around masked areas (e.g. ice-free terrain, clouds). The effect is desirable because MOD10A1 close to such regions is often unreliable (e.g. cloud shadows).

4.3. Finding slush limit candidates

The slush limit is the elevation where $\sigma _\alpha$ falls below a certain threshold (Greuell and Knap, 2000). Using a threshold in $\sigma _\alpha$ as the sole criteria for detecting $\Upsilon _{\rm S}$, however, leads to a large percentage of erroneous detections. Hence, further criteria are needed (Greuell and Knap, 2000). Our approach is grouped around the central premise of a threshold in $\sigma _\alpha$. We even found the threshold of 1.25% identified by Greuell and Knap, to perform well on the MOD10A1 derived grids of $\sigma _\alpha$. However, all other conditions were developed independently. In the following, we first explain the basic workflow of the algorithm and then detail our set of conditions applied to find $\Upsilon _{\rm S}$.

The search for daily $\Upsilon _{\rm S}$ is done along 83 longitudinal stripes. These stripes have a north-south extent of 20 km. Grid cells that fall into a stripe are grouped in 20 m elevation bins. In the following, we label the nth elevation bin $z_c( n)$.

For each day and stripe we aim at identifying the elevation bin that corresponds to $\Upsilon _{{\rm S}}$. For a given day and stripe, the algorithm starts searching for $\Upsilon _{\rm S}$ when <40% of the grid cells are cloud covered. The following parameters are calculated for each $z_c( n)$ based on all grid cells that fall into that bin: (1) median $\sigma _\alpha$ ($M_{\sigma \alpha }$), (2) mean $\alpha$ ($\mu _{\alpha }$), (3) the 95th percentile of NDWI$_{{\rm ice}}$ (NDWI$_{95{\rm th}}$) and (4) mean elevation ($\mu _z$). These parameters are subsequently used to determine whether elevation bins can be considered candidates for $\Upsilon _{\rm S}$.

The algorithm scans through all elevation bins. An elevation bin $z_c( i)$ is retained as a candidate for $\Upsilon _{\rm S}$ if it fulfills all of a set of seven conditions. Below we first list the seven conditions, followed by detailed explanations.

1. 1. $z_c( i)$ is at least seven elevation bins (=140 m) below the uppermost elevation bin, and at least 140 m above the lowermost elevation bin.

2. 2. Cloudiness (i.e. percentage of masked grid cells) $\le 25$% for all seven elevation bins above and below $z_c( i)$.

3. 3. All seven elevation bins above $z_c( i)$ have $M_{\sigma \alpha } < 1.25$% AND at least one of the four bins below $z_c( i)$ has $M_{\sigma \alpha } > 1.65$% AND all four bins below $z_c( i)$ have $M_{\sigma \alpha } > 1.25$%

4. 4. The mean of the $\mu _{\alpha }$ of the five elevation bins above $z_c( i)$ is larger than the mean of the $\mu _{\alpha }$ of the seven elevation bins below.

5. 5. The mean of the $\mu _{\alpha }$ of the seven elevation bins below $z_c( i)$ is smaller than 0.72.

6. 6. The mean of the $\mu _{\alpha }$ of $z_c( i) ^{ + 5 \, {\rm bins}}_{-7 \, {\rm bins}}$ is larger than 0.52.

7. 7. The mean of the NDWI$_{95{\rm th}}$ of the seven $z_c( n)$ below $z_c( i)$ is at least 0.0075 larger than the mean NDWI$_{95{\rm th}}$ of the seven $z_c( n)$ above $z_c( i)$.

Condition 1 guarantees that an elevation bin is not too close to the ice margin or the uppermost bin (the algorithm evaluates each elevation bin based on up to seven bins directly above and below). Condition 2 prohibits any detection of $\Upsilon _{\rm S}$ in the vicinity of larger clouded areas (prone to erroneous $\alpha$, $\sigma _\alpha$ and NDWI$_{{\rm ice}}$).

Condition 3 is key as it defines the possible positions of $\Upsilon _{{\rm S}}$; all subsequent conditions serve to either confirm or reject suggested positions. The third condition evaluates where spatial variability of $\alpha$ drops below the threshold value. Figure 1 exemplifies this drop in $M_{\sigma \alpha }$. However, the algorithm has to deal with various realizations of $\Upsilon _{\rm S}$ that are often not as obvious as in Figure 1. Actual $\Upsilon _{\rm S}$ are characterized by low $M_{\sigma \alpha }$ over broad elevation ranges, starting directly above $\Upsilon _{\rm S}$ (condition 3 – all seven elevation bins above $z_c( i)$ have $M_{\sigma \alpha } < 1.25$%). Below $\Upsilon _{\rm S}$, $M_{\sigma \alpha }$ sometimes remains elevated over larger distances (cf. Fig. 1), sometimes drops off over relatively short distances (condition 3 – only four bins below $z_c( i)$ need to have $M_{\sigma \alpha } > 1.25$%). Finally, there needs to be a certain contrast in $M_{\sigma \alpha }$, otherwise there is increased risk of erroneous detection (condition 3 – at least one of the four bins below $z_c( i)$ has $M_{\sigma \alpha } > 1.65$%).

Fig. 1. Example of a slush limit $\Upsilon _{\rm S}$ as mapped on 14 July 2015 close to 67 $^{\circ }$N. The detected $\Upsilon _{\rm S}$ (20 m elevation bin) is shaded in blue in subplots (a) to (c). (a) Filtered MOD10A1 albedo $\alpha$ in %. (b) Standard deviation of albedo ($\sigma _\alpha$; in %) calculated along horizontal and vertical lines of pixels. (c) NDWI$_{{\rm ice}}$ calculated from MOD09GA. (d) Median of $\sigma _\alpha$ (%), mean of $\alpha$ (%) and 95th percentile of NDWI$_{{\rm ice}}$ calculated for all 20 m elevation bins. Gaps in the curves are where cloudiness exceeds a quarter of the gridcells in an elevation bin. The vertical dotted cyan line indicates the chosen $\Upsilon _{{\rm S}}$; gray shading illustrates the maximum width ($\pm$7 elevation bins) of the search window used to detect $\Upsilon _{{\rm S}}$. Horizontal lines illustrate thresholds for $\alpha$ (in blue) and for $\sigma _\alpha$ (in orange). Threshold values are indicated on the figure.

It would be desirable to formulate condition 3 more strictly. The condition expects $M_{\sigma \alpha } < 1.25$% over a range of seven elevation bins above a potential $\Upsilon _{\rm S}$. In reality, spatial variability of surface albedo is consistently low for all elevations above $\Upsilon _{\rm S}$. Quality issues in the MOD10A1 data played a role when formulating conditions. An elevation interval of seven bins was chosen because highest elevations are frequent subject to spurious grid-like or banded artifacts in MOD10A1 that let $M_{\sigma \alpha }$ rise above 1.25%.

Condition 4 searches for a change in $\mu _{\alpha }$, which is the second key-characteristic of $\Upsilon _{\rm S}$. Conditions 5 and 6 constrain the search to a certain range of $\mu _{\alpha }$. Very high values of $\alpha$ typically do not occur directly above $\Upsilon _{\rm S}$. Very low values directly below indicate the presence of a bare ice surface, not typically to be found at $\Upsilon _{\rm S}$. The search window is more narrow above than below a potential $\Upsilon _{\rm S}$ (five vs. seven elevation bins). This choice was made as the area directly above was found to be particularly indicative of potentially erroneous detections. While changes in $\mu _{\alpha }$ are a key indicator of the location of $\Upsilon _{\rm S}$ (Greuell and Knap, 2000), conditions 4–6 are deliberately kept vague: The threshold in condition 4 is dynamic; conditions 5 and 6 define a wide interval of $\mu _{\alpha }$ (0.52–0.72). Formulating exact and narrow thresholds would limit large-scale applicability of the algorithm and conflicts with limitations in data quality.

Surface water is more abundant below $\Upsilon _{\rm S}$. We introduced a change in NDWI$_{{\rm ice}}$ as the third major characteristic of $\Upsilon _{\rm S}$. Several statistics were tested to optimize condition 7 and we found the change of the 95th percentile of NDWI$_{{\rm ice}}$ to be the strongest indicator. Nevertheless, differences are often small and a rather low threshold of 0.0075 was chosen. Condition 7 reduces the number of erroneous detections and removes weakly expressed slush limits. The latter typically occur early in the melt season at low elevations and leave no expression in NDWI$_{{\rm ice}}$ values.

4.4. Filtering of slush limit candidates

The algorithm sometimes detects more than one candidate for $\Upsilon _{\rm S}$ on a given day and stripe. In such a case only the candidate with the highest $\mu _{\alpha }$ (mean over seven elevation bins below and above a candidate) is retained. Remaining candidates for $\Upsilon _{\rm S}$ require filtering to remove false positives (Fig. 2). We apply an automated approach in two stages.

Fig. 2. Two examples of detected outliers among the slush limit candidates. Latitudes refer to the centers of the stripes.

In a first stage we filter the detections of the year with the highest $\Upsilon _{{\rm S}}$. For all latitudinal stripes this is the year 2012 (see Section 5). In each latitudinal stripe we search for conflicts between the $n$ available candidates. We examine all $( n( n-1) /2)$ pairs of candidates within the stripe: we define as conflicting all pairs where $\Upsilon _{{\rm S}}$ decreases substantially ($> 45$ m). Then, the candidate with the most conflicting pairs is excluded from the sample of candidates. The procedure is repeated until no conflicts are left within the latitudinal stripe. All excluded candidates are labeled invalid.

The filtering for conflicts is highly efficient in detecting outliers that are substantially higher or lower than its neighbors. However, the approach can only test whether the last candidate in a certain stripe falls substantially below its precursors. Any large positive deviation does not create a conflict. For this reason, the last valid candidate of each stripe needs separate filtering (thereby valid refers to having passed filtering for conflicts). We initially test whether the last valid candidate is substantially higher ($> 95$ m) than its valid precursor. If the rate of rise between the second last and the last candidate is also high ($> 9.5$ m d$^{-1}$), then the candidate is considered suspicious. Once a suspicious candidate has been found, the algorithm scans for valid candidates that fall into a window four stripes wide both in direction south and north and $\pm$ 8 days around the date of the suspicious candidate. For each of the eight stripes that fall into the search window, the algorithm retains the valid candidate which is closest in time to the suspicious candidate (if there is any within the $\pm$8 days time window). Consequently, there are up to eight neighboring candidates that are consulted whether they support the suspicious candidate. If they fall within 75 m of elevation of the suspicious candidate, they are considered supporting. If there are fewer than two supporting candidates, then the suspicious candidate is labeled invalid.

In a second stage, we apply filtering to all years. Initially, all candidates of a given stripe are compared to the highest valid $\Upsilon _{{\rm S}}$ in the year of highest $\Upsilon _{{\rm S}}$. Candidates that exceed the highest valid $\Upsilon _{{\rm S}}$ by more than 40 m are labeled invalid. The 40 m threshold is introduced as the detections made during summer 2012 might have missed the absolute peak in some cases. Furthermore, comparison to the year of highest $\Upsilon _{{\rm S}}$ is only carried out if for a given stripe the year of highest $\Upsilon _{{\rm S}}$ has at least four valid detections and if there is at least one detection after 15 July. All candidates that pass the comparison to the year of highest $\Upsilon _{{\rm S}}$ are then filtered identical to filtering in stage 1 (described above).

4.5. Annual maximum slush limits

Deriving annual maximum slush limits $\max \Upsilon _{{\rm S}}$ is challenging because the evolution of $\Upsilon _{\rm S}$ varies strongly between years and regions (Section 5.3). We derive $\max \Upsilon _{{\rm S}}$ per stripe and year as follows: (1) $\max \Upsilon _{{\rm S}}$ are only calculated if $n \ge 5$ detections exist for the year/stripe combination under investigation. (2) All $n$ detections are sorted according to elevation in descending order. (3) Then we assign $\sigma _1 = 0$ to the highest detection, followed by calculating the standard deviation $\sigma _2$ for the two highest detections, again followed by calculating $\sigma _3$ over the three highest samples, until the $n$th and lowest sample is reached ($\sigma _n$). (4) We search local minima $\sigma _{i}$ in the array of $\sigma$ values, i.e. $\sigma _{i-1} \ge \sigma _{i} < \sigma _{i + 1}$ (thereby $\sigma _1$ is always treated as a minimum). (5) Only $\sigma _i$ are retained which fulfill the condition $\sigma _{i} \le 25$ m (6) The $\sigma _i$ with the highest standard deviation is then considered $\max \Upsilon _{{\rm S}}$.

The highest $\sigma _i$ corresponds to the largest group of $\Upsilon _{\rm S}$ whose standard deviation is still acceptable. The approach prefers larger groups of similar $\Upsilon _{\rm S}$ over isolated peaks (Fig. 3). In years where $\Upsilon _{\rm S}$ does climb more or less continuously, $\sigma _1$ typically is the only $\sigma _{i}$ and selected as $\max \Upsilon _{{\rm S}}$. Groups of similar $\Upsilon _{\rm S}$ that lie far below the highest individual value are excluded by the condition $\sigma _{i} \le 25$ m.

Fig. 3. Two examples of calculated annual maxima of the slush limit $\max \Upsilon _{{\rm S}}$. Latitudes refer to the centers of the stripes.

Calculated $\max \Upsilon _{{\rm S}}$ can be too early in the melt season to be reliably considered annual maxima. We exclude all derived annual maxima which take place before 10 July. In case the maxima is based on a group of detections, we test whether the latest date of the group is after 10 July. If this is the case, the $\max \Upsilon _{{\rm S}}$ is retained.

4.6. Comparison to Landsat-derived visible runoff limits

We compare $\Upsilon _{{\rm S}}$ to observations of the visible runoff limit $\Upsilon _{{\rm R}}$ identified from Landsat 30 m near-infrared imagery from 1985 to 2020 (Tedstone and Machguth, 2022). Briefly, the Landsat-derived observations are based on the extraction of supraglacial drainage networks to identify the locations at which visible meltwater runoff occurs. Landsat observations consist of up to 8000 retrievals of the visible runoff limit in any one Landsat scene ($\sim 8000 \times 8000$ pixels).

We compare our individual $\Upsilon _{{\rm S}}$ with $\Upsilon _{{\rm R}}$. As Landsat retrievals use GIMPDEM (Howat and others, 2014) for elevation we first convert MODIS retrievals to use the same source. Then, for each $\Upsilon _{{\rm S}}$, we identify all Landsat retrievals made within the MODIS stripe on the day of observation. If there are Landsat detections available, this usually yields at least several hundred visible runoff limit locations. We aggregate the Landsat locations by taking the median values of x, y and elevation. This yields a dataset of coincident $\Upsilon _{{\rm R}}$ and $\Upsilon _{{\rm S}}$ at the native spatial resolution of the MODIS retrievals.

Next, we filter the dataset to remove $\Upsilon _{{\rm R}}$ which are noisy or not comparable. (i) In each stripe, days on which the $\Upsilon _{{\rm R}}$ has an elevation median absolute deviation >100 m are removed (affecting 5.5% of observations). (ii) It is possible that at any given point in time, runoff features are detected on a Landsat scene situated below the actual $\Upsilon _{{\rm R}}$, while no scene exists that covers the true $\Upsilon _{{\rm R}}$. The highest Landsat detected runoff features could then be mistakenly considered the true $\Upsilon _{{\rm R}}$. In such cases, the identified $\Upsilon _{{\rm R}}$ is located close to an edge of the scene. We therefore remove all paired observations whose Landsat retrieval was <500 pixels ($\sim$15 km) from the edge of the imaged scene. (iii) We filter out inter-comparisons where the MODIS coordinate falls outside of the WRS path/row bounding box of the Landsat scene being examined.

4.7. Meteorological conditions near the K-Transect slush limit

We investigate drivers of annual evolution of $\Upsilon _{\rm S}$ on the example of the K-Transect (67 $^{\circ }$N), a region with good availability of meteorological observations (Section 3). We use cumulative positive degree hours (PDH) for a straightforward characterization of the energy input at the elevation of $\Upsilon _{\rm S}$. While typically positive degree days are used (e.g. Reeh, 1991; Braithwaite, 1994), we here rely on PDH as hourly meteorological observations are available and the use of positive degree days at higher elevation on the ice sheet is subject to issues (van den Broeke and others, 2010).

Cumulative PDH over melt seasons 2009–2021 are calculated based on hourly records of air temperature $T_{\rm a}$ measured along the K-Transect at the PROMICE automatic weather stations (AWS) KAN_U (1840 m a. s. l.) and KAN_M (1270 m a. s. l.). For each year, the PDH at both stations are linearly interpolated to the elevation of the annual $\max \Upsilon _{{\rm S}}$.

5. Results

5.1. Overview

The algorithm determines 20 900 slush limit candidates, out of which the filtering keeps 19 800 ($\sim$95%) as valid $\Upsilon _{\rm S}$.

On average, there are 900 $\Upsilon _{\rm S}$ per year or 10.8 per year and stripe. Detection is most frequent between 65.5 $^{\circ }$N and 70.5 $^{\circ }$N at 15–21 $\Upsilon _{\rm S}$ per year and stripe. For most of the remaining latitudinal stripes there are around 10 $\Upsilon _{\rm S}$ per year. Two intervals (71.5 – 72.25 $^{\circ }$N and 75.25 – 76.5 $^{\circ }$N) have only few $\Upsilon _{\rm S}$.

Maximum slush limits $\max \Upsilon _{\rm S}$ were derived in 72% of all year/stripe combinations (1323 successful retrievals). Success rates vary between years and latitudes and are generally higher for latitudes 65.5 – 69.75 $^{\circ }$N (Figs 4a and c) as well as in warm summers with relatively low cloud cover (e.g. 2012, 2019, cf. Fig. 5).

Fig. 4. Overview of maximum annual slush limits ($\max \Upsilon _{\rm S}$) retrieved for the years 2000–2021. (a) Greenland's west coast and location of median $\max \Upsilon _{\rm S}$. Colors indicate the number of successfully retrieved annual $\max \Upsilon _{\rm S}$. (b) Box plots of elevation (m a. s. l.) of all annual $\max \Upsilon _{\rm S}$. (c) Years with successful retrievals of $\max \Upsilon _{\rm S}$.

Fig. 5. Annual maximum slush limits ($\max \Upsilon _{{\rm S}}$) for the years 2000–2021. (a) Box plot of all $\max \Upsilon _{{\rm S}}$ per year. (b) Number of latitudinal stripes with $\max \Upsilon _{{\rm S}}$. The dotted red line indicates the number of latitudinal stripes (83). (c) Latitudes of retrieved annual $\max \Upsilon _{{\rm S}}$.

5.2. Spatial and temporal trends in annual maximum slush limits

Plotting $\max \Upsilon _{{\rm S}}$ against latitude shows a clear decrease toward the north (Figs 4a and b). The median of $\max \Upsilon _{{\rm S}}$ peaks at $\sim$2080 m a.s.l. at around 63.25 $^{\circ }$N and descends below 1300 m a. s. l. at the northern margin of our study area. Two regions divert from this general trend: (i) south of 63 $^{\circ }$N $\max \Upsilon _{{\rm S}}$ vary strongly with latitude and are on average lower than the adjacent stripes north, and (ii) 69.25–73 $^{\circ }$N where $\max \Upsilon _{{\rm S}}$ show no latitudinal trend.

Highest and lowest $\max \Upsilon _{{\rm S}}$ in each stripe differ by 277 m on average. Inter-quantile distance and spread vary with latitude and peak in the south (Fig. 4b). Plotting $\max \Upsilon _{{\rm S}}$ along the time axis (Fig. 5a) shows the highest median as well as absolute highest values in the year 2012; the lowest median is observed in the year 2001.

The irregular spatio-temporal distribution of retrieved $\max \Upsilon _{{\rm S}}$ (Section 5.1) challenges the assessment of temporal trends in $\max \Upsilon _{{\rm S}}$. We calculate linear trends of annual median $\max \Upsilon _{{\rm S}}$ over all latitudes (most affected by aforementioned issues) and for three regions where data coverage is more constant over time. These are ‘central’ (65.5–69.75 $^{\circ }$N, 24 latitudinal stripes) where data coverage is optimal, ‘south’ (61.75–64.25 $^{\circ }$N; 15 latitudinal stripes) and ‘north’ (72.5–75.0 $^{\circ }$N; 13 latitudinal stripes). Furthermore, we assess trends over the full time period of the simulation as well as for the years 2000–2012 and 2013–2021. These periods were defined based on Ryan and others (2019) who found a significant trend in bare-ice exposure of the ice sheet for the years 2000–2012 but not afterwards.

Table 1 summarizes the results of the regression analysis. All regions show significant positive trends at 95% confidence for 2000–2012. There are no significant trends for 2013–2021. Over the full time period we observe significant trends (at 90% confidence) for ‘south’ and ‘north’.

Table 1. Linear regression of median annual $\max \Upsilon _{\rm S}$ for the time periods 2000–2012, 2013–2021 and all years (2000–2021)

Slope is in m a$^{-1}$. Regressions significant at the 95% confidence level are highlighted in bold, at 90% in normal font. The remaining regressions are in italic.

5.3. Annual evolution of slush limits

The evolution of $\Upsilon _{\rm S}$ during individual summers shows substantial variability. We first assess differences by latitude, followed by the identification of characteristic temporal patterns.

Regionally averaged evolution of $\Upsilon _{\rm S}$ (Fig. 6) north of $\sim$73 $^{\circ }$N indicates a relatively uniform and gentle rise of $\Upsilon _{\rm S}$. The observed melt season is relatively short. Further south ($\sim$64 to $\sim$73 $^{\circ }$N) the observed melt season is substantially longer. In May, the initial progression of $\Upsilon _{\rm S}$ is slow, followed by a steeper rise throughout June. Subsequently the rate of rise decreases and plateaus in early August. South of $\sim$64 $^{\circ }$N the rise of $\Upsilon _{\rm S}$ appears more monotonous. Again, in early August a plateau is reached.

Fig. 6. Annual progression of the slush limit $\Upsilon _{{\rm S}}$ in six regions of the Greenland west coast. Each region comprises 12 latitudinal stripes. Navy blue lines and circles show average behavior calculated from all $\Upsilon _{{\rm S}}$ that fall into a region; pale blue circles denote averages based on <15 individual $\Upsilon _{{\rm S}}$. Cyan lines and dots illustrate progression of $\Upsilon _{{\rm S}}$ in individual stripe/year combinations; for clarity only every 5th stripe/year combinations has been selected randomly. Slush limit progression in red to orange colors is selected manually to illustrate frequent behavior: FL, flat progression; PT, plateauing; SR, steep rise; ON, oscillation.

Figure 6 illustrates how progression of $\Upsilon _{\rm S}$ in individual years differs from average behavior. There are examples where the progression of $\Upsilon _{\rm S}$ appears jagged or very long phases of absolutely linear ‘growth’ occur. There are also cases where $\Upsilon _{\rm S}$ is at $\max \Upsilon _{{\rm S}}$ throughout the entire observed melt season. Such ‘behavior’ is generally associated with years that have extended phases without detection or very few detections. We therefore focus on years that have good data coverage to qualitatively identify reoccurring pattern.

With the exception of the southern-most region in Figure 6, the spread around the mean decreases toward the end of the melt season. Across most latitudes and in many years $\Upsilon _{\rm S}$ reaches a plateau where it remains until detection ceases in autumn (Fig. 6). Less frequently, $\Upsilon _{\rm S}$ continuously climbs until there is no more detection. A peculiar behavior of $\Upsilon _{\rm S}$ is observed at certain latitudes in the south. There, detected $\Upsilon _{\rm S}$ can switch suddenly from a lower to a substantially higher level. The reasons behind this behavior are not yet understood. It could be related to regions where surface discharge networks develop less clearly and runoff also occurs through aquifers. It might be that surface discharge develops quite abruptly when subsurface discharge is overwhelmed by high meltwater input.

5.4. Comparison to Landsat-derived visible runoff limits

We compared $\Upsilon _{{\rm S}}$ to observations of the visible runoff limit $\Upsilon _{{\rm R}}$ identified from Landsat 30 m near-infrared imagery from 1985 to 2020 (Tedstone and Machguth, 2022). Aggregating and filtering (Section 4.6) yields an inter-comparison dataset with 1335 paired observations, or $\sim 7\percnt$ of all 19 800 valid $\Upsilon _{{\rm S}}$. Figure 7 shows the paired observations in a scatter plot, together with a linear regression. The latter is highly significant ($p < 0.0001$) and yields R $^2$ of 0.87.

Fig. 7. Linear regression of Landsat-derived visible slush limits $\Upsilon _{{\rm R}}$ against MODIS-derived slush limits $\Upsilon _{\rm S}$. Number of samples = 1335, slope of the linear regression is 0.857, $R^2 = 0.87$, $p < 0.0001$.

5.5. Slush limit ‘plateauing’ and the summer of 2012

We observe frequent ‘plateauing’ of $\Upsilon _{\rm S}$. To investigate the causes of this phenomenon, we utilize PDH data 2009–2021 from the K-Transect (a region where the phenomenon is frequent). We find that $518\pm 152$ PDH (average$\pm 1\sigma$) have occurred at the location of $\max \Upsilon _{S}$ by the time it is reached (Table 2). The minimum is 314 PDH (summer 2020), the maximum is 872 PDH (summer 2010). In the following, PDH $_\Uparrow$ designates the sum of PDHs consumed to bring $\Upsilon _{\rm S}$ to $\max \Upsilon _{{\rm S}}$.

Table 2. Cumulative positive degree hours at the K-Transect before (PDH $_{\Uparrow }$) and after (PDH $_\Rightarrow$) the maximum slush limit ($\max \Upsilon _{{\rm S}}$) has been reached

All PDH values refer to the elevation of $\max \Upsilon _{{\rm S}}$ in the respective years. The number of days between the first and last detection of $\max \Upsilon _{{\rm S}}$ is labeled plateau. For $\max \Upsilon _{{\rm S}}$, the day of year (DoY) of its first detection and elevation (m a.s.l.) are given.

We investigate whether the plateauing is caused by a lack of additional PDHs once $\max \Upsilon _{{\rm S}}$ has been reached. We calculate PDH $_\Rightarrow$ as the amount of PDH that occur after when $\max \Upsilon _{{\rm S}}$ is first reached and until $T_{\rm a}$ falls permanently below 0 $^{\circ }$C. Excluding 2012, we find that PDH $_\Rightarrow$ is on average$\pm 1\sigma$ 142$\pm$61 PDH. This corresponds to $\sim$27% of PDH $_\Uparrow$.

In contrast, during the summer of 2012, PDH $_\Rightarrow$ amounted to 737 PDH. This value differs strongly from all other years (cf. Table 2) and is even substantially higher than the 530 PDH that occurred in 2012 before $\max \Upsilon _{{\rm S}}$ was reached. Furthermore, $\max \Upsilon _{{\rm S}}$ was reached very early in the season, earlier than in any other year (Table 2).

6. Discussion

6.1. Performance of the algorithm

The data present a clear picture of the spatial and temporal evolution of $\Upsilon _{{\rm S}}$ along the western flank of the ice sheet and the 22 years from 2000 to 2021. The data also provide a detailed image of the annual evolution of $\Upsilon _{\rm S}$ in the majority of the latitudinal stripes and years. The degree of detail, however, varies across latitudes and years.

The total number of retrieved $\Upsilon _{\rm S}$ is limited by frequent cloudiness but also data quality. The latter is most obvious in the years 2000 and 2021 where MOD10A1 data are affected by severe striping and an erroneous cloud mask, respectively. It is unknown by how much the number of $\Upsilon _{\rm S}$ detection could be increased with an optimized algorithm. Dropping conditions increases detection but also the number of false positives. Consequently, stronger filtering would be needed. For example, before the NDWI$_{{\rm ice}}$ threshold was added to the algorithm, the number of $\Upsilon _{\rm S}$ candidates was higher by $\sim$40%. Compared to the final algorithm, the filtering (similar approach as in the final version) removed twice as many candidates, but a substantial number of false positives went undetected. Preference was given to a stricter algorithm that yields a somewhat smaller set of high-quality $\Upsilon _{{\rm S}}$.

The filtering of the $\Upsilon _{{\rm S}}$ candidates works efficiently and there are virtually no obvious outliers left after filtering. Consequently, the data are here presented and used without manual removal of outliers. The filtering is relatively strict and the automated labeling of almost all outliers comes at the costs of (i) having spent a substantial amount of time on optimizing the filter and (ii) removal of some valid candidates.

Our algorithm records the rise of $\Upsilon _{\rm S}$, it does not document lowering of $\Upsilon _{\rm S}$. Decreasing $\Upsilon _{\rm S}$ could occur due to summer cold spells and at the end of the melt season. Indeed, candidates of $\Upsilon _{\rm S}$ often show a decrease in elevation at the end of summer. Such candidates are removed by the filtering for two reasons. Primarily, we consider upward migration of $\Upsilon _{\rm S}$ a clear indicator of a rise in the extent of surface runoff; decreasing $\Upsilon _{\rm S}$ are not necessarily indicative of a lowering limit of surface runoff. Lowering of $\Upsilon _{\rm S}$, as recorded from MODIS data, is often related to summer or late ablation season snow fall events.

While able to mask the hallmarks of $\Upsilon _{\rm S}$, snowfall does not immediately end the physical processes that take place around $\Upsilon _{\rm S}$ (i.e. water saturation in snow, lateral flow in the snow matrix and runoff in streams, cf. Fig. 8). The hydrologic system has substantial inertia as described by Holmes (1955). At $\sim$66.3 $^{\circ }$N/$\sim$47.8 $^{\circ }$W they observed the surface hydrology at the runoff limit over the entire melt season of 1953. They report that melt ceased at 3 August 1953, but main rivers continued to flow at reduced discharge until 21 August when they had frozen over. We therefore believe that optical remote sensing is not optimally suited to observing decreases in slush limits. Secondly, and on a purely pragmatic level, filtering candidates of $\Upsilon _{\rm S}$ becomes easier when decreasing $\Upsilon _{\rm S}$ are removed.

Fig. 8. The hole shown in the photo was dug on 29 July 2020 at a location (47.2391 $^{\circ }$W, 66.9913 $^{\circ }$N; $\sim$1760 m a.s.l.) where a few days earlier (around July 22) a slush field had started to form (all pore space in the surface snow was water filled). Intense snowfall between 24 and 28 July had covered the newly formed slush field. While now hidden below the layer of fresh snow, water continued to flow slowly downhill. Depth of the hole 62 cm (equal to total snow depth), water depth 42 cm. The bottom of the hole and of the snow pack is formed by the ice slab which acts as aquitard and is at the location more than 10 m thick.

6.2. Comparison to Landsat visible runoff limit and other data sets

The comparison to Landsat mapped visible runoff limits $\Upsilon _{{\rm R}}$ (Tedstone and Machguth, 2022) yields a very good agreement with 95% of differences $< 135$ m. This encouraging result is somewhat disturbed by outliers where $\Upsilon _{{\rm S}}$ and $\Upsilon _{{\rm R}}$ deviate by up to 1.2 km (Fig. 7). The character of these outliers was examined manually on a selection of 30 outliers. Absolute deviations of the selection range from 110 to 1170 m, random selection has been done separately for 15 strong outliers (absolute deviations $> 200$ m) and 15 moderate outliers (absolute deviations between 200 and 110 m). We found that 14 of the strong outliers are caused by issues in the comparison (see Section 5.4), not by failure of the Landsat or MODIS detection algorithms. One strong outlier was caused by failure of filtering of $\Upsilon _{{\rm S}}$ candidates. The causes of moderate outliers are more difficult to quantify. About half of them seem related to flawed comparisons, the remainder concerns vaguely expressed slush limits where Landsat $\Upsilon _{{\rm R}}$ is likely more accurate and a couple of cases where the reason is unclear. The analysis shows that it would be justified to remove most of the extreme outliers from the comparison of $\Upsilon _{{\rm S}}$ and $\Upsilon _{{\rm R}}$. Consequently, the statistics presented in Section 5.4 and Figure 7 are conservative. Removal of the 32 outliers $> 200$ m that are not known to be caused by failure of the MODIS algorithm yields R $^2 = 0.96$ and a regression slope of 0.93 ($n = 1303$).

We find good qualitative agreement between $\max \Upsilon _{{\rm S}}$ and AVHRR mapped slush lines of 1990 and 1995 by Greuell and Knap (2000) (Fig. 9). Greuell and Knap (2000) highlight the two years because of strong melting. Indeed, the comparison shows that the slush lines of the two strong melt summers of 1990 and 1995 coincide mostly with the median of $\max \Upsilon _{{\rm S}}$. We restrict this comparison to a qualitative level because the slush lines according to Greuell and Knap (2000) needed digitizing from relatively low resolution figures and accuracy of their digital elevation model is poorly known.

Fig. 9. Comparison of annual maximum slush limits to the runoff limit as calculated by Reeh (1991) as well as AVHRR derived slush lines (Greuell and Knap, 2000) for the years 1990 and 1995.

Reeh (1991) used positive degree day modeling to calculate the snow line, the equilibrium line and the runoff limit for most of Greenland. Greuell and Knap (2000) argue that their slush lines are most representative for the runoff limit and compare their data to Reeh's (1991) runoff limits. We do the same (Fig. 9) and find that the two data sets agree well. The runoff limit according to Reeh (1991) shows fewer details because of the coarse resolution of the model input (mainly accumulation distribution and air temperatures) available at the time.

6.3. Slush limit, visible and actual runoff limit

The Landsat visible runoff limits agree well to MODIS $\Upsilon _{{\rm S}}$, but the question remains to what degree they represent the actual runoff limit, that is the highest elevation where a minimal fraction of melt finds a way to the ocean. Holmes (1955) wrote of ‘white slush’ that is hidden under a layer of unsaturated snow and thus typically invisible to optical remote sensing. While not measured directly by Holmes (1955), their study indicated that lateral runoff in white slush reached at least 4.2 km (2.6 miles) above the highest visible slush fields. Clerx and others (2022) directly measured lateral flow velocity inside the snow, in a shallow layer of meltwater on top of the ice slab. Multiplying measured velocities with a rough estimate of days per melt season where water can flow yields a value similar to Holmes (1955). Tedstone and Machguth (2022) note that their visible runoff limits are 10–20 m lower in elevation than the upper limit of the runoff area at the EGIG line ($\sim$69.7 $^{\circ }$N, $\sim$48.5 $^{\circ }$W), as quantified for the years 2007 and 2008 by Humphrey and others (2012). Based on these references, we cautiously conclude that MODIS $\Upsilon _{{\rm S}}$ provide a reasonably good estimate of the actual runoff limit. Compared to the latter, $\Upsilon _{{\rm S}}$ are likely biased low by a few kilometers. This estimate refers to regions with well-developed ice slabs.

6.4. Temporal trends in maximum slush limit elevation

Our data show weak or insignificant trends in $\max \Upsilon _{{\rm S}}$ (Section 5.2) when focusing on the full time period. This appears to contradict Tedstone and Machguth (2022) who find a substantial increase in the visible runoff limit. However, they investigate 1985 to present, 15 years more than this study. They also show that the increase in the visible runoff limit levels off after the year 2012. The absence of a trend in $\max \Upsilon _{{\rm S}}$ might also be perceived to be at odds with studies that charted the inland progression of supraglacial lakes (e.g. Howat and others, 2013; Gledhill and Williamson, 2017). Similar to Tedstone and Machguth (2022), these studies consider longer time periods. In addition, they exclude the most recent years (e.g. Howat and others, 2013, end with the year 2012).

Ryan and others (2019) assessed the extent of bare ice exposure on the entire Greenland ice sheet from 2001 to 2017. They found no significant linear trend in end-of-summer snowline elevation and bare ice extent over their full study period. In contrast, they observed a significant trend when focusing on 2000–2012. We observe a similar behavior of $\max \Upsilon _{\rm S}$ with a significant rise until the year 2012. In the north and south of our study area $\max \Upsilon _{\rm S}$ increased by $\sim 20$ m a$^{-1}$.

The emerging pattern of a relatively steep increase of $\max \Upsilon _{\rm S}$ until 2012, followed by stagnation, is in broad agreement with the observed evolution of ice content in firn at the K-Transect (Rennermalm and others, 2021). A peak in ice content was measured directly after the series of intense melt years had peaked in 2012. Since then ice content has decreased. However, the analysis of Rennermalm and others (2021) ends with the year 2018, excluding the strong melt summer of 2019.

6.5. ‘Plateauing’ and the slush limit in the summer of 2012

We calculated PDH at K-Transect $\max \Upsilon _{\rm S}$ to investigate reasons behind frequent ‘plateauing’ of $\Upsilon _{\rm S}$. We find that $518\pm 152$ PDH are consumed to bring $\Upsilon _{\rm S}$ to its annual maxima. There are no extreme outliers. Varying amounts of snow could be a reason for the differences between years: the more snow, the more melt is needed to fill the pore space until water becomes visible at the surface. Inaccuracies in the data (see below) likely also contribute to the differences. We calculate that $\sim 27\percnt$ of total annual PDH (excluding the year 2012) occur after $\max \Upsilon _{\rm S}$ has been reached. Apparently, this amount of PDH is insufficient to lift $\Upsilon _{\rm S}$ above $\max \Upsilon _{\rm S}$ but suffices to stabilize $\Upsilon _{\rm S}$, sometimes for more than a month.

It is beyond the scope of this study to investigate how the stabilization of $\Upsilon _{\rm S}$ exactly works. It is possible that, once the hydrologic network is established, the slush limit remains visible on MODIS imagery even with reduced water supply. As observed by Holmes (1955), meltwater mobilized earlier in the season contributes to maintaining late-season $\Upsilon _{\rm S}$. PDH that occur late in the season likely assist in such processes. Alternatively, one might suggest that ‘plateauing’ of $\Upsilon _{\rm S}$ marks the upper extent of a near-surface aquitard, such as solid glacier ice or an ice slab. Indeed, Greuell and Knap (2000) observed that in three warm summers melting continued while their AVHRR-derived slush lines had already stopped at a certain elevation (roughly at 1720 m a.s.l. at the K-Transect). They suggested that this behavior is related to sudden changes in the density of the subsurface. In the following we discuss this hypothesis in the context of the 2012 melt season.

The year 2012 reached the highest $\max \Upsilon _{\rm S}$ ($\sim 1840$ m a.s.l. at the K-Transect; Table 2 and Fig. 5), but stands out even clearer in terms of PDHs at $\max \Upsilon _{\rm S}$. There was an extraordinary sum of PDH that occurred after $\max \Upsilon _{{\rm S}}$ had been reached (Table 2). Why did $\Upsilon _{\rm S}$ not rise any further, given abundant melt energy? There are either data issues or a physical reason for $\Upsilon _{\rm S}$ stagnating.

Along the K-Transect, all years have a high number of $\Upsilon _{\rm S}$ retrievals without substantial data gaps around the date when $\max \Upsilon _{\rm S}$ is reached. There are years where detection of $\Upsilon _{\rm S}$ stops relatively early (among them the year 2012). Nevertheless, undetected rise in $\Upsilon _{\rm S}$ is unlikely because neighboring latitudes do not show marked increase of $\Upsilon _{\rm S}$ after detections along the K-Transect ceased. Estimated PDH rely on only two AWS and interpolation to $\max \Upsilon _{\rm S}$ is done linearly. Still, errors in the interpolation are probably small because the elevation difference between $\max \Upsilon _{\rm S}$ and the KAN_U AWS is $< 140$ m for all but one year. In the year 2012, $\max \Upsilon _{\rm S}$ was directly at KAN_U. Consequently, we assume that uncertainties in PDH are of limited influence.

Spatial variations of firn properties, namely a decrease in ice content, are very likely to have stopped the 2012 rise in $\Upsilon _{\rm S}$. At elevations below $\max \Upsilon _{\rm S}$ the ice slab acted as near-surface aquitard, provoking visible runoff. At elevations $> \max \Upsilon _{{\rm S}}$ the ice slab is discontinuous or absent. In these areas meltwater percolated vertically into the porous firn rather than forming slush fields and streams. Our assumption is based on $\sim$40 firn cores, drilled at and above KAN_U in 2012, 2013 and later years (Machguth and others, 2016; MacFerrin and others, 2019; Rennermalm and others, 2021). Radar profiles from spring 2013 show the ice slab extending approximately as high as KAN_U (Machguth and others, 2016; MacFerrin and others, 2019).

6.6. Long-term trends in slush limit elevation

The year 2012 is the only year with clear evidence of firn properties controlling $\max \Upsilon _{{\rm S}}$. Firn data from around $\max \Upsilon _{{\rm S}}$ and before 2012 are very scarce. This makes it difficult to prove whether in other years before 2012 $\max \Upsilon _{{\rm S}}$ was collocated with ‘a sudden change in the density profile’ (Greuell and Knap, 2000). For the time period 1990–1995, Greuell and Knap estimated the ‘abrupt change in the density profile’ at $\sim$1720 m a.s.l. Since then, the ‘sudden change in the density profile’ has migrated upslope to 1840 m a.s.l. Greuell and Knap (2000) emphasize that positions of their slush lines are subject to an uncertainty of $\pm 2$ km (furthermore, the DEM was of 2 km spatial resolution). Nevertheless, the increase of $\sim$120 m at the K-Transect appears substantial and might answer a question formulated by Greuell and Knap (2000), namely to what degree the location of the ‘abrupt change in the density profile’ varies as a function of climate.

How did the ‘abrupt change in the density profile’ look like in the early 1990s? Ice slabs of the western flank of the ice sheet are typically transient phenomena (Leone and others, 2020) in a warming climate. The climatic conditions preceding the measurements of Greuell and Knap are characterized by an extended period (1930 s to early 1990 s) of slight cooling (Van As and others, 2016; Cappelen, 2021). The slush limit of the early 1990s might have been associated with the upper margin of continuous glacial ice, not with ice slabs. Indeed, a radar profile from spring 2013 (Machguth and others, 2016; MacFerrin and others, 2019) shows the uppermost appearance of continuous ice (down to at least 20 m depth) at $\sim$1700 m a.s.l. At higher elevations, the ice slab overlies porous firn.

It is tempting to estimate the potential $\max \Upsilon _{{\rm S}}$ of summer 2012 at the K-Transect. At the elevations of $\max \Upsilon _{\rm S}$, the sum of PDH is $657\pm 159$ at end of melt season (average$\pm 1 \sigma$; all years excluding 2012). Using linear interpolation we find that in 2012 the same amount of cumulative PDH occurred at $2154\pm 82$ m a.s.l. Assuming the ice slab would have extended to higher elevations, $\max \Upsilon _{{\rm S}}$ could have risen $\sim 300$ m above the observed $\max \Upsilon _{\rm S}$ (1841 m a.s.l., cf. Table 2). This estimate is supportedby a total of 850 PDH measured in 2012 at the GC-Net station Dye-2 at $\sim 2120$ m a.s.l., $\sim$60 km south of the K-Transect. Our analysis underlines the exceptional character of the summer 2012 and shows that surface runoff from deep within the current accumulation area can take place as soon as firn density profiles have become favorable.

7. Conclusions

We designed an automated algorithm to detect the slush limit along most of the western margin of the Greenland ice sheet. The approach was applied to daily MODIS data covering the study area and all years of the MODIS record (2000–2021). We find very good agreement with visible runoff limits mapped from higher resolution Landsat data (Tedstone and Machguth, 2022).

Our results document spatial variations of the slush limit with high accuracy and yield good temporal detail. We observe frequent ‘plateauing’ of the slush limit, that is the slush limit stops rising clearly before the end of the melt season. While we find weak or absent trends in the slush limit over the entire time period, there is a clear and significant linear trend that peaks and ends in the year 2012.

Our algorithm is ‘classical’ in the sense that we do not rely on computer learning or artificial intelligence (AI). Advantages are replicability and potentially process understanding, disadvantages could be fewer detected slush limits and a higher number of false positives. There might be potential in exploring our concept in an AI framework.

Our data shed new light on the extraordinary characteristic of the summer of 2012. On the example of the K-Transect (67 $^\circ$N) we showed that the rise of the slush limit, and thus also the extent of surface runoff, stopped early in the season at the upper margin of the ice slabs. If firn conditions would have allowed, the slush limit would have risen further, potentially up to 2154$\pm$82 m a.s.l. While this number is still hypothetical, it underpins concerns about the future of the Greenland ice sheet.