Using ICESat-2 and Operation IceBridge altimetry for supraglacial lake depth retrievals

Supraglacial lakes and melt ponds occur in the ablation zones of Antarctica and Greenland during the summer months. Detection of lake extent, depth, and temporal evolution is important for understanding glacier dynamics. Previous remote sensing observations of lake depth are limited to estimates from passive satellite imagery, which has inherent uncertainties, and there is little ground truth available. In this study, we use laser altimetry data from the Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) over the Antarctic and Greenland ablation zones and the Airborne Topographic Mapper (ATM) for Hiawatha Glacier (Greenland) to demonstrate retrievals of supraglacial lake depth. Using an algorithm to separate lake surfaces and beds, we present case studies for 12 supraglacial lakes with the ATM lidar and 12 lakes with ICESat-2. Both lidars reliably detect bottom returns for lake beds as deep as 7 m. Lake bed uncertainties for these retrievals are 0.05–0.20 m for ATM and 0.12–0.80 m for ICESat-2, with the highest uncertainties observed for lakes deeper than 4 m. The bimodal nature of lake returns means that high-confidence photons are often insufficient to fully profile lakes, so lower confidence and buffer photons are required to view the lake bed. Despite challenges in automation, the altimeter results are promising, and we expect them to serve as a benchmark for future studies of surface meltwater depths.

1 Introduction

The ice sheets of Antarctica and Greenland modulate rates of sea level rise, contributing 14.0 ± 2.0 mm (Antarctica) and 13.7 ± 1.1 mm (Greenland) since 1979 (Mouginot et al., 2019; Rignot et al., 2019). Current trends indicate greater melt in the coming decades, leading to the contributions from both ice sheets to overtake the contribution of thermal expansion to sea level rise (Vaughan et al., 2013). Meltwater plays vital roles in ice sheet evolution (e.g., van den Broeke et al., 2016), including aggregation on ice sheets as supraglacial lakes, many of which are several meters deep (Echelmeyer et al., 1991). When unfrozen, these lakes exhibit a lower albedo than that of the surrounding ice, allowing them to absorb more incoming solar radiation and melt ice more efficiently, thus generating a positive feedback (Curry et al., 1996). Supraglacial lakes are significant reservoirs of latent heat (Humphrey et al., 2012), and their spectral emissivity in the infrared (IR) spectrum also differs from bare ice (Chen et al., 2014; Huang et al., 2018), which can lead to potentially significant impacts on the surface energy balance of ice sheets.

A substantial portion of meltwater eventually drains into supraglacial streams or moulins (drainage channels), where it can flow to the ice bed (Banwell et al., 2012; Catania et al., 2008; Selmes et al., 2011). During catastrophic lake drainage events, meltwater penetration into the ice can also lead to hydrofracture, a mechanism through which meltwater facilitates full ice fracture as a result of the stresses induced by the density contrast between liquid water and ice (Das et al., 2008). Meltwater injection to the bed can also modify basal water pressures which in turn modify the resistance to ice flow and thus can impact sliding velocity and ice discharge. (Parizek and Alley, 2004; Zwally et al., 2002). Hydrofracture can lead to significant ice loss for outlet glaciers and ice shelves (Banwell et al., 2013). Current observations and modeling efforts indicate a propagation of supraglacial lakes farther inland as the climate warms (Howat et al., 2013; Leeson et al., 2015; Lüthje et al., 2006), raising further concerns for accelerated mass loss. For these reasons, knowledge of supraglacial lakes is important for our understanding of ice sheet evolution.

Previous studies developed techniques for detecting supraglacial lakes and retrieving depth, areal coverage, and volume. In situ observations employed sonar and radiometers to approximate lake depth and albedo (Box and Ski, 2007; Tedesco and Steiner, 2011). However, the harsh conditions of Antarctica and Greenland, the transience of meltwater, and the sheer size of the ice sheet ablation zones restrict the potential for extensive in situ measurements, encouraging lake depth and areal coverage estimates from passive remote sensing data such as Landsat-8, MODIS, and Sentinel-2 A/B. Supraglacial water is darker than surrounding ice in the visible and IR bands, allowing the use of band ratios between blue and red reflectance (Stumpf et al., 2003). The normalized water difference index (NWDI) and dynamic thresholding techniques have also been considered for lake detection (Fitzpatrick et al., 2014; Liang et al., 2012; Moussavi et al., 2016; Pope, 2016; Williamson et al., 2017; Moussavi et al., 2020). Other methods implemented radiative transfer models (Georgiou et al., 2009) or positive-degree-day models (McMillan et al., 2007) to estimate lake albedo and meltwater volume. By comparing surface reflectance data of supraglacial water to that of ice and optically deep water, empirical relationships have been derived to approximate lake depth (Philpot, 1989; Sneed and Hamilton, 2007).

Image-based empirical techniques rely on approximations of lake bed albedo and an attenuation parameter, both of which are subject to uncertainties from lake heterogeneity and cloud cover (Morassutti and Ledrew, 1996). Furthermore, Pope et al. (2016) found that band ratios were insensitive to lakes deeper than 5 m, leading to errors that may exceed 1 m. Parameter fitting in the empirical equations requires supplementary depth retrievals, often from in situ sources. More accurate methods for supraglacial lake detection are needed to improve image-based estimates.

In September 2018, the Ice, Clouds, and land Elevation Satellite-2 (ICESat-2) was launched with the primary objective of obtaining laser altimetry measurements of the polar regions (Abdalati et al., 2010; Markus et al., 2017; Neumann et al., 2019b). Observations using the Airborne Topographic Mapper (ATM) and Multiple Altimeter Beam Experimental Lidar (MABEL) indicated the potential for shallow water profiling with laser altimetry (Brock et al., 2002; Brunt et al., 2016; Jasinski et al., 2016), and ICESat-2 applications were recently demonstrated by Ma et al. (2019) and Parrish et al. (2019). In this study, we identify test cases from ICESat-2 and ATM altimetry data and use these pilot cases to develop an algorithm for detecting supraglacial lakes and retrieving lake depth. The algorithm is designed as a semiautomatic method to find supraglacial lakes within select altimetry granules.

2 Data description

2.1 ICESat-2

ICESat-2 is a polar-orbiting satellite with an inclination of 92^∘ that carries the Advanced Topographic Laser Altimeter System (ATLAS), a 532 nm micro-pulse laser that is split into six distinct beams with names based on the ground track: GT1L/R, GT2L/R, and GT3L/R. The beams are configured in pairs with a 90 m separation between beams within a beam pair and 3.3 km between pairs. With an operational altitude of ∼ 500 km and a 10 kHz pulse repetition rate, ICESat-2 records a unique laser pulse approximately every 0.7 m along-track over a 91 d repeat cycle.

The ATLAS product used here is the ATL03 Global Geolocated Photon Data V002 (Neumann et al., 2019a), which consists of retrieved photons tagged with latitude, longitude, received time, and elevation. Each photon pulse also carries a classification as either signal or background (noise). The differentiation between signal and background is performed using a statistical algorithm outlined by Neumann et al. (2019b). Signal photons are further classified by confidence level, such that photons labeled as “high confidence” are most likely to originate from the surface. Generally, cloudy or variable profiles exhibit “medium/low confidence” or noise photons, whereas low-slope surfaces, such as water and ice sheets, result in more high-confidence photons (Neumann et al., 2019b). In thin layers of water, high-confidence photons are observed from both the water surface and the underlying ice.

Our study focused on the central strong beam (GT2L), as the number of lakes was deemed sufficient for our purposes. While we recognize that the other strong beams could be useful for depth retrievals, we did not consider them here. We speculate that the weak beams may avoid issues with multiple scattering and specular reflection, but their power is too low to reliably detect lakes deeper than 4 m. Ground-based validation by Brunt et al. (2019b) indicates an accuracy of <5 cm in ATL03 photons over ice sheet interiors. The use of medium, low, and “buffer” photons slightly decreases measurement precision, but a less truncated transmit pulse gives better agreement with ATL06 and ground-based data (Brunt et al., 2019b).

2.2 Airborne Topographic Mapper

The Operation IceBridge (OIB) campaign was designed to fill the gap in polar altimetry between ICESat and ICESat-2. Its scientific payload included the Airborne Topographic Mapper, a 532 nm lidar that has been used for ice sheet and shallow water measurements since 1993. The ATM lidar conically scans at 20 Hz, providing a 400 m swath width along-track (Brock et al., 2002; Krabill et al., 2002). The ATM Level-1B Elevation and Return Strength (ILATM1B) product converts analog waveforms into a geolocated elevation product to emulate ATLAS data (Studinger, 2013). Although it lacks a statistical confidence definition, ATM applies a centroid model to digitized waveforms to retrieve high-confidence photons. Brunt et al. (2019a) found that ATM errors were −9.5 to 3.6 cm relative to ground-based measurements. Here, the ATM results presented serve as a proof of concept for the lake detection algorithm.

3 Methods

3.1 Lake detection

Supraglacial lake surfaces are much flatter than surrounding terrain. We thus performed topography checks with the expectations that (i) lake surfaces would be identifiable in photon histograms and (ii) lake beds may be found via statistical inference in the region of the lake surface. To simplify the identification of lake features, we separated them into two arrays: one for the surface and one for the bed, which we refer to as “lake surface–bed separation” (LSBS). For both lidars, the procedure for separation was identical and is as follows (see Fig. 1 for a schematic view).

Figure 1 Schematic for the workflow of the lake surface–bed separation algorithm.

• i.

  We divided each data granule into discrete along-track windows to reduce the data volume to ∼ 10⁴–10⁵ photons per window. This photon count is equivalent to ∼ 1–10 km in the along-track distance for ICESat-2 and ∼ 0.15–1.5 km for ATM. If a supraglacial lake appeared on the edge of the window, the window size was adjusted to include the full observed water feature.

• ii.

  Each data window was binned into elevation-based histograms. We assumed that the lake surface dominates the total bin count within each window of photons. We check the flatness of the window by computing the standard deviation (σ) of high-confidence signal photons within the upper 85th percentile of bin count. We define a “flat” surface for regions where σ≤0.05 m for ATL03 data and ≤ 0.02 m for ILATM1B data. We selected these values by comparing the flatness of lake surfaces to that of surrounding ice topography. If data were within the appropriate flatness threshold, they were verified as a lake surface using Landsat-8 OLI imagery. This step was included to filter non-glacial features, such as ocean or fjords.

• iii.

  If the satellite image(s) confirmed the presence of a lake, the data were assigned to a new array for the height of the lake surface (h_sfc). The horizontal extent of the lake surface served as a constraint for where the lake bottom data could be defined. Within these horizontal bounds, photons were defined as a lake bottom if they satisfied the condition $h_{\mathrm{sfc}}-a{\mathit{\sigma }}_{\mathrm{sfc}}\le h\le h_{\mathrm{sfc}}-b{\mathit{\sigma }}_{\mathrm{sfc}}$, where σ_sfc is the standard deviation of lake surface photons. The constraints a and b were derived through trial and error, such that a=1.0(1.8) and b=0.5(0.75) for ICESat-2 (ATM). We set these constraints to reduce the impacts of multiple scattering and specular reflection on depth estimates. If these conditions were met, then the data were placed in an array for the height of the lake bottom, h_btm.

• iv.

  A series of filters were applied to improve surface–bed estimates. For ICESat-2, lakes shallower than 1.3 m or smaller than 200 m in horizontal extent were considered too noisy or ill-defined for further analysis (see Sect. 5.2 for more details). To remove water bodies with deep bed returns (e.g., oceans or fjords) or with no bed returns, the algorithm counted the number of bed photons present for both lidars. If the number of bed photons was very small (100 or less), then the scene was marked as a probable false positive.

• v.

  If the data were obtained from ICESat-2, then we followed a photon refinement routine that is described in more detail in Sect. 3.2. Calculations for lake depth were then performed for both ATM and ICESat-2 retrievals and corrected for refraction (Sect. 3.3).

3.2 ATL03 refinement

The above steps were sufficient to obtain lake profiles within the ATM data, but melt lake bottoms observed by ICESat-2 were significantly noisier as a consequence of higher background (noise) photon rates. After the initial LSBS procedure, we manually assessed bed estimates for each lake. For lakes that did not pass qualitative assessment, we adopted photon refinement procedures initially used for the ATL06 surface-finding algorithm (Smith et al., 2019). In short, ATL03 photon aggregates within overlapping 40 m segments were used to estimate lake surfaces and beds with greater precision via least-squares linear fitting applied to the aggregates. These linear fits were used to approximate a window of acceptable surface or bed photons for every 20 m along-track. A more detailed description of the ATL06 algorithm is given in Smith et al. (2019).

The linear regression in ATL06 accounts for all ATL03 photons (background or signal), and the technique performs a background-corrected spread estimate to narrow the range for acceptable photons. This is an iterative scheme; the refinement process repeats its acceptable photon filter until no photons are removed. As a consequence, the ATL06 algorithm assumes a single returning surface, so over a melt lake it will compute a height for either the lake bottom or the lake surface, depending on their return strengths.

The condition for acceptable surface photons in ATL06 is given by

$$\begin{array}{}\text{(1)}& \left|r-r_{\mathrm{med}}\right|<\mathrm{0.5}{H}_{\mathrm{w}}.\end{array}$$

Within a 40 m photon segment, r is the residual of a photon relative to the linear regression, r_med is the median residual, and H_w is window height. The height of the window is taken as the maximum of the observed photon spread, the window height (if any), and 3 m, and photons within the window range are defined as the surface. The LSBS algorithm follows a similar procedure, but the flatness of the lake surface and relatively low photon density of the corresponding beds rendered iteration unnecessary. The lake bed is then defined as photons not within the window and below the surface. In other terms, lake bed photons satisfy the conditions

$$\begin{array}{}\text{(2)}& \left|r-r_{\mathrm{med}}\right|>\mathrm{0.5}{H}_{\mathrm{w}},h<h_{\mathrm{sfc}}.\end{array}$$

As with the initial guess, the lake bottom was only defined within the horizontal bounds of the lake surface, and the improved guesses were assigned to h_sfc and h_btm.

As a final adjustment to lake photons, we applied a refraction correction algorithm to account for slowing down of light as it enters water. The correction follows the methods utilized by Parrish et al. (2019) by approximating refractive biases as a function of depth and beam elevation angle. The center strong beam for ICESat-2 is near the nadir, so the horizontal offset was determined to be small relative to the size of lakes (∼ 3 cm, far below the horizontal geolocation uncertainty for ICESat-2). However, vertical offsets of 1 m or more were found for lakes ≥4 m in depth, necessitating the use of refraction correction.

3.3 Lake depth and extent estimations

Once we obtained h_sfc and h_btm, lake depth from the altimeter signal (z_s) was estimated using

$$\begin{array}{}\text{(3)}& z_{\mathrm{s}}={\stackrel{\mathrm{‾}}{h}}_{\mathrm{sfc}}-{\stackrel{\mathrm{‾}}{h}}_{\mathrm{btm}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{sfc}}$ and ${\stackrel{\mathrm{‾}}{h}}_{\mathrm{btm}}$ represent the moving mean of the surface elevation and the bottom elevation, respectively. The moving mean was used to account for signal attenuation and scattering at the lake bottom, a problem most evident for ICESat-2 retrievals.

For deep or inhomogeneous lakes, attenuation of photon energy in water resulted in fewer signal photons observed at lake bottoms (Fig. 4). In these situations, we fitted polynomial or spline fits to all lake profiles with bounds at the lake edges. Lakes observed by ATM typically featured “bowl” shapes and attenuation at the deepest parts, so third-order polynomials were sufficient. In ICESat-2 data, the retrieved lake beds showed greater complexity, so we tested polynomial fits and splines on a case-by-case basis. Lake depths approximated with curve fitting were denoted as z_p. We compare z_s and z_p over lakes with well-defined bottoms, and we show in Sect. 4 that the two generally agree to within 0.88 m.

To test the limits of the algorithm relative to lake size, we utilized the great-circle formula (ATM) or predefined along-track distance (ICESat-2) to approximate along-track extent L. We acknowledge the desire to retrieve lake volume from laser altimetry, but we leave the development of such an algorithm for a future study. For example, depth retrievals from ICESat-2 could potentially be combined with lake radius and shape estimations determined from visible satellite imagery to derive water volume.

3.4 Case study locations

We present cases over the Amery Ice Shelf on 2 January, 2019 (ICESat-2 Track 0081; 68.271–73.798^∘ S, 63.057–78.620^∘ E), the western Greenland ablation zone for 17 June 2019 (ICESat-2 Track 1222; 66.575–69.582^∘ N, 48.284–49.239^∘ W), and Hiawatha Glacier on 19 July 2017 (ATM; 77.780–79.3119^∘ N, 65.279–67.484^∘ W) (Fig. 2). Comparisons between Landsat-8 imagery and ICESat-2/OIB flight tracks confirmed supraglacial lake overpasses for study. In spring 2019, an early onset of the Arctic melt season resulted in both ICESat-2 and Operation IceBridge surveying supraglacial lakes near Jakobshavn Isbræ in May. However, there were no lakes sampled at the time by both ICESat-2 and OIB.

Figure 2 True-color Landsat-8 composites of Hiawatha Glacier on 18 July 2017 (a), the Amery Ice Shelf on 1 January 2019 (b), and the western Greenland ablation zone on 17 June 2019 (c). Flight tracks for Operation IceBridge (a) and ICESat-2 (b, c) are shown in dotted orange.

4 Results

We detected 12 melt lakes with sufficient bed returns from the ATM data and 16 potential melt lake surfaces overall. The melt lake profiles are shown in Fig. 3, with maximum depths of 0.98–7.38 m and extents of 180–730 m. The algorithm reliably distinguishes between lake surfaces and the surrounding ice terrain. The mean spread among lake surface photons is 0.0087 m, or well within the flatness threshold of 0.02 m. Lake bottoms are well-defined when d_s<8 m. Lake bottoms deeper than 8 m exhibit fewer signal returns, for the associated return signal is below the threshold required to be digitized (Martin et al., 2012). The average depth estimate for the lakes in Fig. 3 was 1.95 m (Table 1), and lakes at this depth typically featured adequate bed returns. In deeper lakes, the polynomial estimate produced reasonable guesses for the lake bed location, with the most effective fitting seen in lakes 3e, 3g, and 3h. With the polynomial-based depths, mean lake depth increased to 2.15 m, and the maximum modeled depth was 8.83 m.

Figure 3 ATM lake profiles from 17 July 2017 fitted using lake surface–bed separation, including the raw ILATM1B product, the lake surface signal, the lake bottom signal, the polynomial- and spline-fitted bottom, and the point of maximum depth. Along-track distance is relative to the start of a data granule.

Table 1 Cumulative statistics for ATM supraglacial lakes explored in this study, including mean and maximum signal-based depth (d_s) and polynomial-based depth (d_p), along-track extent L, mean lake depth uncertainty (${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{d}}$), and mean polynomial estimation error (${\stackrel{\mathrm{‾}}{\mathit{ϵ}}}_p$). Units are in meters.

The spread in ATM lake bed photons is low (Table 1, column 7), with a maximum of 0.2 m for lake 3g. The highest uncertainties are observed for lake depths greater than 3 m, perhaps influenced by low signal-to-noise ratios or the conical scanning of the OIB lidar instrument. Polynomial estimation errors are 0.41 m on average. Several depth errors are below this mean, but a strong standard error (1.03 m) in lake 3g, due to difficulties in capturing its steep bed slope, slightly skews the mean error. Excluding this value, the mean error among ATM polynomial estimates reduces to 0.35 m.

We examined an additional 12 supraglacial lakes with ICESat-2, eight in Greenland and four on the Amery Ice Shelf in Antarctica. Three of the Antarctic melt lakes (4a, 4b, 4d) are highlighted in Magruder et al. (2019) and Fricker et al. (2020). The refined algorithm captures lake surfaces and beds reasonably well (Fig. 4), with a mean uncertainty of 0.015 m for surface photons and 0.38 m for bed photons. The lake edges partially account for the bed photon uncertainty, for the limited number of acceptable photons produces a slight bias in bed estimates. Antarctic melt lakes were generally shallower than those seen in Greenland (Table 2) – only lake 4a exceeded 3 m in depth, whereas the mean maximum depth over Greenland was 4.08 m. Melt lakes on the Amery Ice Shelf were 3–8 km in extent, thus facilitating detection in histograms. Greenland lakes exhibited a wider range of sizes, but the algorithm successfully performed retrievals for lakes as small as 200 m in extent.

Figure 4 Supraglacial lakes and melt ponds detected by ICESat-2 over the Amery Ice Shelf (a–d) and western Greenland (e–l), using Tracks 0081 and 1222, respectively.

Table 2 As with Table 1, but for ICESat-2 tracks.

On average, the noisier data from ICESat-2 produce uncertainties greater than 0.2 m for the Antarctic lakes and 0.3 m for the Greenland lakes, as seen in Table 2, column 8. The inclusion of lower-confidence photons increases uncertainty despite the restricted bed photon criteria, for the larger photon cloud increases the spread of the entire lake profile. The curve fits improved depth estimates for lakes 4b, 4f, and 4i. Of these lakes, only 4i used a polynomial estimate due to poor spline fitting. The inclusion of interpolants increased the mean depth estimates of 4b, 4f, and 4i by 0.08 m, 0.04 m, and 0.03 m, respectively. The spline fitting significantly increased the maximum observed depth in lake 4b from 2.67 m to 3.27 m. The remaining lakes featured more complete bed profiles, meaning that the fitting estimates were less important.

5 Discussion

5.1 Algorithm performance

The conical scanning of the ATM lidar produced oscillations in 1D elevation profiles that dampened over lake surfaces, so lakes generally were easier to identify with the airborne retrievals. Flights conducted during the OIB campaign actively avoided cloudy conditions, reducing attenuation sources and further simplifying the lake-finding process over common melt regions. The data volume per granule was lower than ATL03, resulting in less time needed to run the algorithm. However, the number of retrievals possible with ATM is limited, so observations with the lidar best serve as a validation and correction tool for ICESat-2 and other retrieval methods.

The laser power and detector sensitivity of the ATLAS instrument on board ICESat-2 are sufficient to reliably detect lake beds, and a high along-track resolution will correspond to improved estimates of lake bed topography, water depth, and water volume. Despite strong advantages, significant difficulties must be considered before automatic lake detection is feasible. At its operational altitude, the ATLAS laser is subject to first-photon bias, solar background radiation, and scattering and absorption by blowing snow and clouds. Clouds are common over the fringes of Antarctica and Greenland (Bennartz et al., 2013; Lachlan-Cope, 2010; Van Tricht et al., 2016), and often their optical depth is sufficient to render the surface undetectable. Handling the large data volumes in ATL03 granules also presents a significant challenge. A single granule provides coverage over hundreds of kilometers, so the running time of the algorithm increases relative to ATM granules. Lakes smaller than 1 km are difficult to automatically detect with the algorithm, but LSBS may still be performed for lakes as small as 200 m if the location of a lake is known through other means (e.g., Landsat-8 imagery or ATM retrievals).

We observed differences in lake topography for ICESat-2 lakes, and we attribute them to the underlying ice surfaces. Supraglacial lakes in Greenland typically form into smooth basins within depressions formed by the underlying bedrock, and their location is independent of ice motion (Echelmeyer et al., 1991). In contrast, meltwater on the Amery Ice Shelf originates from the blue ice zone, propagating along the ice surface in streams. The location of lakes and ice topography is thus tied to the flow lines of the ice shelf surface. These features are flooded in the Antarctic melt season, producing melt lakes and streams up to 80 km in length (Mellor and McKinnon, 1960; Phillips, 1998; Kingslake et al., 2017).

A potential issue for lake depth retrievals concerns specular reflection. When photons interact with a flat water surface, they may reflect directly back to the detector with minimal energy loss. The excessive return energy produces a “dead time” in the ATLAS detector, and the return signal is represented by multiple subsurface returns below the actual surface (Neumann et al., 2020). An example of this phenomenon may be seen in Fig. 4f, where a prominent subsurface return 1 m below the true surface is featured along the lake extent. However, because the subsurface echo is smaller than the true surface when viewed through histograms, the LSBS algorithm is able to avoid biases caused by specular reflection.

The success of this method for lake depth retrievals is governed by spatial and temporal sampling of the instruments across the lakes when they are full. The methods presented here are most effective when the altimeter passes directly over the deep part of a lake rather than at its edge. This provides a lake depth profile that is more representative of the complete lake, allowing for improved estimates of lake depth and extent. A complete lake profile also provides sufficient information to the LSBS algorithm, reducing the risk of false negatives that occur with small lakes or incomplete profiles. The temporal sampling of ICESat-2 and ATM is infrequent (every 91 d for ICESat-2 and random for ATM), and so the same lakes will not always be present every time these data are required. Therefore, coincident satellite imagery is desirable to simplify the lake-finding process.

5.2 Automation challenges

The identification of lake beds in the LSBS algorithm is based on a window of acceptable photons. The photon window is constrained by the coefficients a and b (for ICESat-2, a=1.0, b=0.5). Lake beds detected in this manner had a height uncertainty of 0.38 m (Table 2). The coefficients for ATM (a=1.8, b=0.75) resulted in more accurate retrievals on an individual basis. However, implementing varying a and b values proved difficult to automate, as other values may produce more accurate depths. The challenges in full automation are related to three key issues. First, the observed extent of lakes varied considerably, especially over Greenland. The diversity in lake sizes complicated attempts to derive a universal flatness check. Smaller lakes present fewer lake surface photons, so a smaller data window (∼ 10⁴ photons) is required to prevent false positives. However, larger lakes may not be fully represented in smaller windows. A larger data window (∼ 10⁵ photons) will fully capture the largest lakes, but smaller lakes may then be overlooked.

Second, multiple scattering at the lake bed increases the photon spread and thus also increases the uncertainty of depth retrievals. Most supraglacial lakes observed by ATM featured smooth beds, so photons experienced one or few scattering events before returning to the detector. The instrument digitizer automatically filters return signals with low photon counts, reducing the spread of bed photons, at the cost of deep lake bottom detection. In contrast, the lakes observed with ICESat-2 exhibited more heterogeneous beds, leading to increased scattering events by photons and delays in return pulses. In these cases, the given values for a and b may not produce the most accurate bed solution. Furthermore, if the return is significant for a given photon window, then it may lead to a false negative for a portion of the lake (Fig. 4i). To reduce uncertainty in lake depth retrievals, future improvements in working with ICESat-2 data should focus on identifying and filtering multiple scattering.

Finally, the ATL03 signal-finding algorithm is conservative in that it accepts false positives (background photons classified as signal photons) to ensure that all signal photons are passed to higher-level products. Thus, uncertainties in the ATL03 photon classification contribute to noise in the water column and the lake bed. The classification algorithm uses predefined surface masks to allocate statistical confidence to ATL03 photons for multiple surface types (e.g. inland water, land ice, land), with overlap possible between masks (Neumann et al., 2020). Melt lakes are categorized as land ice (lake surface) and land (lake surface and bed). Because the land classification also includes the bed, it includes more potential signal photons than land ice, so our recommendation is to only use land photons for supraglacial lake depth retrievals. It must be noted, however, that a lake bed profile is fully resolved only with the inclusion of low- and medium-confidence and buffer photons. The buffer photons ensure that all photons identified as surface signal are provided to the appropriate upper-level data product algorithms. However, they can introduce greater noise to the profile, so more sophisticated filtering techniques are needed to distinguish between signal photons and the solar background.

6 Conclusions

We present a method to detect supraglacial lakes and estimate lake depth from 532 nm laser altimetry data. We establish test cases for lake detection over two regions of Greenland (Hiawatha Glacier, 19 July 2017 and Jakobshavn Isbræ, 17 June 2019) and East Antarctica (Amery Ice Shelf, 2 January 2019), and our results demonstrate that depth retrievals are possible using laser altimetry. Verification of lake detection is given with lake surface flatness tests, where we observe low topographical variance over lake surfaces relative to surrounding ice. Lake bottoms are easy to identify once lake surfaces are established, given that the lakes are not deeper than 7 m.

We introduce a lake surface–bed separation scheme for ATM and ICESat-2 geolocated photon data to determine the maximum depth of lakes. Our results indicate that altimetry signals reliably detect bottoms as deep as 7 m, after which absorption of the photons in water reduces the number of reflected photons. Heterogeneity at the lake bed also produces attenuation, complicating retrieval attempts for lakes with rough bed topography or with high impurity concentration. Additional work is required to assess the impacts of lake impurities and geometry on altimetry signals and to improve estimates for such cases. Despite these shortcomings, we anticipate retrieval capability to improve as observations from the 2019 and 2020 Arctic melt seasons are released.

We establish the feasibility for estimates of supraglacial lake depth over Antarctica and Greenland. The high accuracy of 532 nm laser altimeters allows these results to serve as a benchmark for future retrieval studies. Future studies need to examine the accuracy of ICESat-2 lake retrievals relative to ATM where applicable, with additional comparisons to depth estimates from passive imaging sensors.