Estimating Supraglacial Lake Depth in Western Greenland 1 Using Landsat 8 and Comparison with Other Multispectral 2 Methods 3

Abstract


Conclusions References
Tables Figures

Introduction and rationale
Supraglacial lakes in Greenland play a crucial role in the glacial hydrological system.Together with supraglacial streams (Smith et al., 2015), lakes temporarily store large quantities of meltwater which can promote the opening of conduits to the bed through hydrofracture and thus influence ice dynamics (Das et al., 2008;Phillips et al., 2013;Selmes et al., 2011;Tedesco et al., 2013).Surface runoff has been increasing in Greenland over the past years (Seo et al., 2015), and supraglacial lakes have been shown to Introduction

Conclusions References
Tables Figures

Back Close
Full establish hydrological connections to the glacier bed and speed up flow, both through observational tools (Joughin et al., 2013) and models (Parizek and Alley, 2004).In Antarctica, supraglacial lakes can contribute to ice shelf disintegration, possibly through positive feedback mechanisms (Banwell et al., 2013;Glasser and Scambos, 2008).Supraglacial lakes are also important for surface heat fluxes by storing latent heat near the surface of the ice sheet.In addition, lakes reduce surface albedo, hence further promoting melt through a positive feedback mechanism (Leeson et al., 2015).
Several multispectral remote sensing tools and methods exist for both classifying (Johansson and Brown, 2013; Leeson et al., 2013;Sundal et al., 2011) and estimating the depth of supraglacial lakes (Sneed and Hamilton, 2011).MODIS (the MODerate Resolution Imaging Spectroradiometer) is able to provide large spatial coverage (2330 km swath width), moderate resolution (∼ 250 m) images twice per day (e.g., Box and Ski, 2007;Fitzpatrick et al., 2013).ASTER (the Advanced Spaceborne Thermal Emission and Reflection Radiometer, e.g.Sneed and Hamilton, 2007) and Landsat (e.g.Banwell et al., 2014;Morriss et al., 2013) have higher spatial resolution (10-30 m) but lower spatial coverage and fewer acquisitions (16 day repeat).Commercial sensors, such as WorldView-2, provide high resolution multispectral measurements (∼ 2 m) that can be used to image relatively small water features such as streams over smaller areas (17 km wide swath), at both high temporal and spatial resolution during the daylight period (Chu, 2014;Legleiter et al., 2014;Smith et al., 2015).However, commercial imagery is collected largely "on demand" and cloud cover can still be a confounding factor.For supraglacial lake depth estimation with all of the above multispectral sensors, a large scale validation has not been undertaken to date.Lake depth retrieval is based upon the intuitive understanding that deeper water will absorb more energy and therefore will have lower reflectance values.Some methods use one band for a reflectance-depth relationship, while others use a ratio of reflectances from two different spectral bands (see Sect. 2).However, satellite retrieval of supraglacial lake depth faces several limitations, including difficulty measuring the true reflectance of dark/deep lakes, assumptions inherent in the method about minimal Introduction

Conclusions References
Tables Figures

Back Close
Full quantities of suspended and dissolved matter in lake water, a smooth (i.e., not windroughened) lake surface, and homogeneous and low-slope lake bottoms (Sneed and Hamilton, 2011).We further assume that we can extrapolate locally calibrated coefficients to new areas (e.g., Legleiter et al., 2014), and can ignore minor variations in effects of atmospheric path radiance.
With the successful launch of Landsat 8 in 2013 (Irons et al., 2012;Roy et al., 2014), a new and improved multispectral sensor is available for lake depth estimation.The Landsat 8 Operational Land Imager (OLI) has enhanced radiometric resolution (12-bit vs. 8-bit), a higher signal to noise ratio, and an expanded dynamic range as compared to Landsat 7's Enhanced Thematic Mapper Plus (ETM+).While published studies have largely used red and green wavelengths, OLI's two additional bands (coastal, 0.433-0.453 µm; and cirrus, 1.360-1.390 µm) and narrower multispectral and panchromatic bands relative to ETM+ will provide more spectral information and more unique (i.e., less auto-correlated) reflectance values, respectively.These properties lead to improvements for methods based on band ratios.Finally, an increased scene collection rate by Landsat 8 will lead to more opportunities to observe ice sheets and their supraglacial lakes.
In this paper we investigate refinements on retrieval methods for supraglacial lake depth from Landsat 8 imagery.We use in situ spectral measurements from a supraglacial lake in Greenland to emulate satellite reflectance, and then we compare them with depth data from a supraglacial lake to test several techniques to extract lake depth.Based on this analysis, we then apply the best methods to Landsat 8 OLI imagery for case study areas in Northwest Greenland and the Sermeq Kujalleq (Jakobshavn Isbrae) area.We validate depth estimates using digital elevation models (DEMs) derived from stereo sub-meter imagery.We discuss best practices for deriving lake depths using Landsat 8 and the implications of these conclusions for other multispectral sensors.Analysis of 2014 imagery yields information about supraglacial lake size, distribution, and seasonal behavior.Introduction

Conclusions References
Tables Figures

Back Close
Full

Physically based lake depth
The depth of a supraglacial lake can be approximated as (after Philpot, 1989): where z is lake depth (units: m), R lake is the reflectance of a lake pixel, R ∞ is the reflectance of optically deep water, A d is the lake bottom albedo, and g is related to the losses in upward and downward travel through the water column (units: m −1 ).This method has been used successfully in several studies in both Greenland and Antarctica (e.g., Banwell et al., 2014;Sneed and Hamilton, 2007).Since it is based upon a description of the processes that take place as light enters, passes through, and exits a lake, it is physically based and therefore easy to adjust if measurements of lake water and lake bed properties are available.However, it requires optically deep water (i.e., no detectable reflected component from the lake bed) to be present in the images.It further assumes that lake water has little to no dissolved or suspended matter and would be severely impacted by surface waves (wind-driven ripples, choppy waves, etc.).Additionally, it requires that the lake bottoms have low slopes and a homogeneous albedo (Sneed and Hamilton, 2011).While most of these assumptions hold for supraglacial lakes in Greenland (Sneed and Hamilton, 2011), lake bed surfaces are known to be too inhomogeneous to support the approach generally.In addition, optically deep water is not always available in inland Landsat scenes.The effects of these shortcomings on supraglacial lake depth retrievals have not been quantified.In this study, for application to Landsat 8 imagery, R ∞ was obtained from dark ocean or lake water in the scene, following Sneed andHamilton (2007, 2011).If no coast was available in the scene containing the lake, R ∞ was obtained from another scene further along the path (with an implicit assumption of similar atmospheric conditions).The g parameter was calculated following earlier studies (Smith and Baker, 1981;Sneed and Introduction

Conclusions References
Tables Figures

Back Close
Full  , 2007), but with an updated absorption coefficient from Pope and Fry (1997, Table 3); for more details, see the Supplement.

Hamilton
As the spatially closest available estimate for lake bottom albedo, A d was obtained from the reflectance immediately outside identified lake areas.However, in order to test this assumption, we also solve for lake bottom albedo rather than assuming it to be the same as the surrounding ice.In an application of a form of spectral mixture analysis (Lillesand et al., 2007), we define a fractional coverage of ice (r i ) and cryoconite (r c = 1 -r i ) in each lake bottom pixel.(Cryoconite is accumulated small surface debris; Wientjes et al., 2011.)We then use these coefficients and end-members of ice and cryoconite reflectances to calculate the lake bottom albedo (A d ).To create a determinable equation after introducing this new unknown (r i ), we use reflectances from two OLI spectral bands (indicated with subscripts 1 and 2, below), and derive end-member reflectances for ice (R i1 or R i2 ) and cryoconite (R c1 or R c2 ) using glacier reflectance spectra from Pope and Rees (2014b) in conjunction with OLI spectral response functions in both bands (Barsi et al., 2014).We input these parameters into Eq.( 1) and then combine the expressions by equating lake depth, thus obtaining: After Eq. ( 2) is solved for r i , the bottom albedo for one OLI spectral band can be calculated and subsequently used to compute lake depth: where R lake1 is water leaving reflectance (as in Eq. 1) for the first band in the pair used and z is lake depth.Introduction

Conclusions References
Tables Figures

Back Close
Full

Empirically derived lake depth
The second method we consider uses spectral band ratios and has been used to derive water depth in shallow marine settings (e.g., Dierssen et al., 2003) and alluvial rivers (e.g., Legleiter and Overstreet, 2012), and has been also adapted for use on the Greenland ice sheet (Legleiter et al., 2014).While the physically based method above is highly dependent on A d and g, earlier studies show that the band ratio method is expected to be more robust to variations in these parameters (Legleiter et al., 2009;Stumpf et al., 2003).This is because the method relies of relative behavior in two different wavelengths, as opposed to absolute optical behavior.This band ratio method employs an empirically derived quadratic formula to relate lake depths to the ratio of the reflectance of two spectral bands (R 1 and R 2 ): This empirical/band ratio method requires the derivation of calibrated coefficients (i.e.a, b, and c), and coefficients vary depending based on which sensors and bands are used (Legleiter et al., 2014).We calculate these coefficients using a known set of reflectances and depths (from in situ measurements, see Sects.3.1 and 4.1).

Data
Three datasets are used in this study: in situ reflectance spectra and lake depth, Landsat 8 multispectral imagery, and DEMs derived from stereo WorldView imagery.We use in situ data to test different lake retrieval methods for a range of spectral bands.Then, we calculate lake depth with a range of the most promising methods using OLI imagery.
We then use WorldView DEMs to validate the OLI-derived lake depths.The detailed workflow of software (including MATLAB and GDAL code) used for data analysis and Introduction

Conclusions References
Tables Figures

Back Close
Full presentation in this study will be fully described and documented in a subsequent paper.

In situ data
In situ concurrent and bathymetric measurements were collected using a small remotecontrolled boat with a payload including a compact spectroradiometer and a small sonar to collect in situ, coincident lake-bottom reflectance and depth over one lake in west Greenland in the summer of 2010 (Tedesco and Steiner, 2011).We use 2226 unique sample points from that study to evaluate the performance of the remote sensing methods described above.Field spectra are convolved to account for the spectral response of the spaceborne sensors as follows: where r nb is the narrowband reflectance, r(λ) is the spectral reflectance, R(λ) is the relative spectral response (Barsi et al., 2014), and λ is the wavelength.In order to emulate sensor dynamic range and radiometric resolution, we impose minimum and maximum reflectances and round reflectance values to the appropriate precision (i.e., 8-bit or 12-bit; see Pope and Rees, 2014a).We then regress the convolved reflectances and in situ depth measurements to test the goodness of fit of the physically based relationship presented in Eq. ( 1) and the empirical method described in Eqs. ( 5) and (6).

Landsat 8 imagery
Landsat 8 was launched 11 February 2013 and became operational on 30 May 2013 (Roy et al., 2014).OLI collects spectral data gridded at 30 m spatial resolution (15 m for panchromatic data).Top-of-atmosphere (TOA) reflectance is calculated using calibration coefficients provided in the image metadata and a solar elevation cosine cor-Introduction

Conclusions References
Tables Figures

Back Close
Full  , 2013).Based on a sensitivity analysis of path radiance to water vapor and ozone using an atmospheric radiative transfer model (see Sect. 5), images were not atmospherically corrected.Two study areas are chosen for applying OLI imagery (see Fig. 1).One site located in northwest Greenland (including Sverdrup Gletsjer, Dietrichson Gletsjer, Sermersuaq,     and Kjer Gletsjer, on Melville Bay; 56.2966-58.7186• W, 74.9685-75.7808• N) is examined in 2013 as it as an area with a high concentration of lakes and was imaged 4 times by Landsat 8 throughout the summer.A larger region further to the south is examined in 2014 using all available Landsat 8 scenes collected over the Sermeq Kujalleq (Jakobshavn Isbrae) region in West Greenland.For a list of all OLI scenes used in this study, see Table S2.
Using the calculated TOA reflectances, we define supraglacial lake extent using the ratio between the blue and red bands (Banwell et al., 2014;Box and Ski, 2007).However, since OLI bands are slightly different from those of past sensors, we could not use published thresholds for extent.We set the threshold for this ratio at 1.5 based (vs.1.05-1.25 for ETM+ in Banwell et al., 2014) upon visual comparison with the imagery.We then visually inspected and manually adjusted the threshold mask to remove coastal water areas (i.e., not on the ice sheet) and clouds.Although Leeson et al. (2013) describe such thresholding as too coarse for low resolution imagery (i.e.MODIS), they do acknowledge its utility for higher resolution imagery (i.e.ASTER, Landsat, etc.).Regions 4 pixels or smaller (i.e.small lakes likely comprised solely of mixed pixels) or less than 2 pixels wide (i.e.linear features likely to be channels, not lakes) were removed from the lake mask.
We interpolate the lake mask using a nearest neighbor algorithm, in order to apply the physically based method to the higher resolution panchromatic band.Where both panchromatic and spectral bands were used together, we bilinearly interpolate the panchromatic image to 30 m resolution.Introduction

Conclusions References
Tables Figures

Back Close
Full

Worldview DEMs
We use submeter (∼ 0.5 m pixel −1 ) stereo imagery from DigitalGlobe's WorldView-1 and WorldView-2 to create DEMs of lake areas both before filling and after drainage.These high resolution DEMs are generated using the open source NASA Ames Stereo Pipeline tool (Moratto et al., 2010).For both the Sermeq Kujalleq (Jakobshavn) and northwest sites, DEMs from 6 different days were used, for a total of 12 DEMs (see Table S2).
WorldView-1 image data have a geolocation accuracy of better than 4.0 m horizontal 90 % circular error of probability and WorldView-2 better than 3.5 m (DigitalGlobe, 2014).Thus, the imagery and DEMs are more precisely positioned than the 15-30 m Landsat 8 pixels.
The vertical accuracy of the derived DEM products is less than 5.0 m 90 % vertical error of probability with submeter relative vertical precision (Mitchell, 2010).Differencing a WorldView DEM with an airborne topographic mapper lidar profile over a pronounced basin in northeast Greenland provided a standard deviation over the spread of elevations of 0.25 m.Considered conservatively, differencing one WorldView DEM with a second DEM collected one year later provided a standard deviation of 0.58 m for the elevation differences (Willis et al., 2015).A stack of 13 overlapping WorldView-1 and WorldView-2 DEMs (∼ 17 km × 17 km) over Summit Station, Greenland provides absolute vertical accuracy estimates of ∼ 2.0 m (±1.2 m) relative to airborne LiDAR measurements, according to preliminary studies.After removing absolute offsets, the per-pixel standard deviation of elevations for the stack are ∼ 15-30 cm.
We resample the DEMs to the same grid as Landsat imagery using cubic interpolation.The Landsat and WorldView acquisitions are from different dates, and although lake basins do ablate during the summer, this should not have significant impact on the results presented here, as most supraglacial lakes in Greenland remain fixed over bedrock-controlled surface depressions (Lampkin and VanderBerg, 2011).Using the lake mask, we identify a shoreline for a given date (see Sect. to derive lake depth.We remove outliers of impossibly shallow (i.e.negative depth) or deep (> 65 m) values as errors in the DEM.In addition, we remove lakes having a standard deviation in lake elevation along the shoreline of larger than 1.5 m.These steps also mitigated any potential bias caused by temporal offset between DEM and spectral depth measurements.
After all filtering, over 250 000 pixels (30 m) from 6 days in 2013 and 6 days in 2014 remained for spectral lake depth validation.

In situ results
The results of depth-reflectance regressions for all methods are displayed in Table 1 and Fig. 2. The bands tested here using in situ data were based upon those identified in the literature (e.g., Box and Ski, 2007;Sneed and Hamilton, 2007;Tedesco and Steiner, 2011), as well as the OLI's new coastal band and the significantly narrowed panchromatic band (0.500-0.680 µm, at 15 m spatial resolution).ETM+ high and low gain results are virtually indistinguishable, and so only low gain results are shown here.
For each regression, we use the correlation coefficient (r) and the root mean square error (RMSE, relative to sonar depths) to assess the success of each method.The results of the physically based method show that the OLI blue and coastal bands do not perform well relative to other bands (RMSE of 3.10 and 11.03 m, respectively; r of 0.29 and 0.05, respectively).The OLI Band 3 (green, 0.525-0.600µm; 0.78 m, r = 0.78) performs as well as legacy ETM+'s Band 2 (green, 0.525-0.605µm; 0.77 m, r = 0.79).Finally, both OLI Band 4 (red, 0.640-0.670µm) and Band 8 (panchromatic, 0.500-0.680µm) bands outperform their analogous Landsat 7 bands (RMSE of 0.28 and 0.63 m, respectively; r of 0.96 and 0.84, respectively).
Red light attenuates more strongly in water than green or blue light.So, for the same lake depth, there will be a larger (and easier to measure) change in net reflectance for Introduction

Conclusions References
Tables Figures

Back Close
Full red wavelengths than shorter wavelengths.However, the rapid attenuation of red light means that only shallower lakes may be measured in this band.The maximum in situ lake depth measurement is ∼ 5 m, well within the red light limit, but deeper lakes may exist in the overall study area.We address this issue below by using many Landsat scenes and WorldView DEMs.
We investigate the two-band physically based method (where A d was calculated) with a range of emulated OLI bands (see Table 1).We find similarly high correlation coefficients (r = 0.94) to the regression method.Nevertheless, only the combination of blue and green bands had an RMSE below 1 m.This method appears to slightly overestimate lake depths.We investigate the reasons for this with the Landsat and WorldView data below.
Applying the empirical method using field data (see Table 1, Fig. 2) indicates that the more continuous bands of the ETM+ outperform the narrower (less spectrally autocorrelated) bands of the OLI when estimating lake depths.However, the addition of the coastal band should allow the OLI to still perform quite well (r > 0.92, RMSE < 0.38), in particular when paired with the green or panchromatic bands.
Our analysis shows the potential of using Landsat 8 to retrieve supraglacial lake depth as well as or better than with Landsat 7. We identify the best methods for OLI (identified with asterisks in Table 1) based on the highest correlation coefficients and lowest RMSEs.We then apply these methods to Landsat 8 data and validate them with WorldView stereo DEMs.

2013 Northwest Greenland results
In the northwest Greenland study area, we identified 694 lakes on 2 July 2013 (day of year 183) with a total area of 27. most constant as lake growth areas moved higher in elevation over the following three weeks, and then decreased again toward the end of August as cooler conditions prevailed (see Fig. 3).While all methods show the same pattern of surface water storage, the total water volumes returned by the different methods differ by over a factor of 2. We investigate which method agrees best with the WorldView DEMs, providing some confidence in the retrievals of that method.

Comparison with DEMs
We difference overlapping areas of DEM-and Landsat-derived lake depths.The statistics of this comparison are shown in Fig. 4. As seen in the northwest Greenland case study, the methods are divided into two groups.Landsat-derived depths using band 3, bands 2 and 3, a ratio of bands 1 and 3, and a ratio of bands 1 and 8 all considerably overestimate lake depth relative to the DEMs.However, the physically based single band method for the red band (OLI Band 4) only slightly underestimates lake depth (−0.1 ± 1.7 m), while the panchromatic band (OLI Band 8) slightly overestimates lake depth (0.1 ± 1.4 m).
Combining these two best-performing bands, the resulting spectral and DEM-derived lake depths are in very good agreement, showing a difference of 0.0 ± 1.6 m.We infer that the optimal method for estimating supraglacial lake depth with Landsat 8 is to take an average of the physically based (see Eq. 1) depths as derived from the red and panchromatic channels (bold in Table 1).It is likely that the spread in depths is the result of a combination of factors including temporal offset between DEM and spectral data collection, ice flow, image coregistration, and atmospheric effects, as well as uncertainties inherent in the lake depth retrievals.Therefore, while having a fairly large (1.6 m) uncertainty for individual pixels, the mean lake depth derived from these methods agrees well.

Conclusions References
Tables Figures

Back Close
Full The lake depth algorithm (i.e., average of single band depths from OLI red and panchromatic bands) was applied to 34 Landsat 8 scenes from the summer of 2014 over the Sermeq Kujalleq (Jakobshavn) area (see Figs. 1 and 5).The total meltwater storage in supraglacial lakes peaked near 3 km 3 in mid-July 2014.A histogram of all lake depths shows that there are many shallow lakes (0.3 to 1.5 m depth), the distribution has many depths at 2.5 to 4 m, with few lakes exceeding 5.5 m depth (see Fig. 6a).
The preponderance of shallow lake pixels reflects the fact that the observed lakes have low surface slopes at their edges.
The Sermeq Kujalleq (Jakobshavn) dataset provides a timeseries that shows lake growth and drainage/freezing (see Fig. 5a).There are many factors that contribute to lake growth and drainage, including temperature, insolation, albedo, topography, and ice dynamics.These complex drivers are related to the more easily quantified mean elevation and latitude of each scene.For example, isolating the coastal scenes shows the delayed onset of melt and earlier shutdown in the north compared to the south (see Fig. 5b).
To further refine our investigation of geographic factors associated with lake depth over the summer season, we examine single swaths of Landsat imagery through time.Path 008 (in the WRS-2 reference scheme, Irons et al., 2012) which transects the lower Sermeq Kujalleq (Jakobshavn Isbrae) shows a strong influence of both elevation and latitude in rates of lake growth and water storage (Fig. 5c).Isolating Path 006, on the other hand, conflates the effects of elevation and altitude on surface meltwater storage, but because we have more temporal coverage (see Fig. 5d) we see the decline of total lake volume as summer progresses toward autumn.Again, higher latitude and elevation delay melt onset (i.e.Path 006, Row 012).For 006/013 and 006/014, it is likely that the reduced ice sheet area within 006/014 is the explanation for the reduced meltwater volume.Rates of increase and decay of lake volume are similar for this pair.Introduction

Conclusions References
Tables Figures

Back Close
Full The distribution of lake depths (by pixel) with elevation is shown in Fig. 6b.Lakes are distributed from ∼ 300 to ∼ 2100 m elevation.Maximum lake depths occur at about 1200 m a.s.l.At lower elevations, lake depths recorded by our method vary significantly, likely due to rapid lake growth and drainage across a range of dates at lower elevations, vs. the higher elevation maximum depths mostly derived from a Landsat scene on 30 July 2014.From 1200 to 2100 m, measured lake depths decline steadily with less variation.This likely reflects a combination of factors, including the variations in induced surface topography of the ice sheet as it flows over undulating bedrock (Lampkin and VanderBerg, 2011).At higher elevations, slow flow leads to low-amplitude ice surface topography as well as less available meltwater.In addition, while lakes are less likely to variably fill and drain at higher elevations, there was also reduced imagery available from ∼ 30 July 2014 onwards.Therefore, the more consistent maximum depths at higher elevations are a combination of incomplete temporal coverage and.Further down, more melt and higher amplitude topography from faster ice flow facilitate lake formation.However, below 1200 m, increased ablation begins to reduce this topography.In addition, the volume of melt available will determine whether depressions are large enough to hold lakes or instead drain via connecting supraglacial channels.The melt volume and therefore the relationship between lakes and channels will thus vary both seasonally and with elevation as well (Lampkin and VanderBerg, 2014).

Discussion
In the above results section, we presented data showing promising lake depth retrieval techniques, used Landsat imagery and WorldView DEMs to test these methods, and applied imagery towards an understanding of the supraglacial hydrology of a section of the Greenland ice sheet.In discussion below, we consider factors which control the success and/or failure of both lake depth retrieval methods (i.e.physically based and empirically derived) as well as the merits of particular wavelengths of light for lake depth estimation.We then apply this understanding to assess earlier lake depth retrievals in Introduction

Conclusions References
Tables Figures

Back Close
Full the literature that used ETM+ and ASTER data.We next include a sensitivity analysis of atmospheric corrections of OLI imagery and finish with a discussion of limitations to OLI lake depth retrievals.The depths returned by the empirical (band ratio) method considerably overestimate lake depths relative to the WorldView DEMs.The method is entirely dependent upon the calibration of the input parameters (i.e., a, b, and c).The parameters used in this study are in turn based solely upon extrapolation from in situ measurements at a single lake.Therefore, it is possible that the lake used for calibration is not representative of lakes in our study region.Legleiter et al. (2014) note that the coefficients for the empirical method may be scale-dependent, and values calculated from field data may not be appropriate for the 30 m pixels of Landsat 8. Indeed, preliminary work (a separate study led by co-author Moussavi) both calibrates and validates spectrally derived depths with WorldView DEMs to show that the band-ratio/empirical method and singleband/physically based method perform similarly well.The use of a ratio of coastal and green reflectances performed well for lake depth retrieval using WorldView-2 imagery (Legleiter et al., 2014).Therefore the band ratio method may, with better parameters, perform similarly to the physically-based single-band approaches.
Nevertheless, the physically based depth retrievals show a large spread in total water volume returns.Physically based depth retrievals rely on accurate bottom albedos (A d ) and water absorption coefficient (g).While A d is derived from the imagery, g is always calculated for each spectral band based on laboratory measurements and is therefore consistent across all OLI scenes.Comparison of laboratory-measured g values with g values derived from in situ data (see Table 1) shows that when the laboratorymeasured g is higher than the one obtained from regressing in situ data, lake depths are overestimated and vice versa.For example, OLI Band 3 (green) shows a 70 % difference in directly measured and regressed g, and it overestimates lake depths by a mean of 2.4 ± 2.1 m relative to WorldView DEMs.By contrast, Band 4 (red) and Band 8 (panchromatic) have very small differences between measured and regressed Introduction

Conclusions References
Tables Figures
Water absorption properties also vary with wavelength.For example, poor performance in blue and coastal bands is related to very low absorption.Red wavelengths attenuate relatively quickly in water, and this is described by a relatively high g (0.7507 m −1 ) compared to green (0.1413 m −1 ).This high g value for red light makes it less sensitive to errors in g than green wavelengths.Lake depth estimates using a red channel are also less sensitive to A d than with a green channel (Tedesco and Steiner, 2011), again due to the high absorption for longer wavelengths.Ultimately, as long as the sensor radiometry is able to measure the return from deep-water pixels, longer wavelengths (i.e., red) can return generally more accurate lake depths because they are less sensitive to the input parameters.
To evaluate other studies in the literature and compare them to our results, we applied the same methods we use (i.e., lab-measured absorption/scattering parameters and appropriate spectral response functions) to calculate g values for Landsat 7 ETM+ bands (see Table 1).Tedesco and Steiner (2011) studied the accuracy of ETM+'s green band for lake depth estimation.They tested different multipliers of the diffuse attenuation coefficient for downwelling light to get the water absorption coefficient g.They showed that for Landsat 7 ETM+'s green band, sonar and spectral depths correlated better when a larger multiplier was used.This is broadly consistent with the 70 % offset between observed and theoretical values that we observe (Table 1).They also find that this offset "cannot be easily explained, aside from a possible chlorophyll concentration in the water, currently considered to be unlikely."Morriss et al. ( 2013) used ETM+'s red band and extracted a higher value of g (0.86 m −1 ); this is very close to the regressed value we observe of 0.83 m −1 (see Table 1), and so we expect their depth estimates to be slightly overestimated.Banwell et al. (2014) and Arnold et al. (2014) also used Landsat 7's green band with a g of 0.1954 m −1 , ∼ 40 % percent higher than our regressed value of 0.14 m −1 , probably leading to significantly overestimated lake depths, estimated to be ∼ 30 %.Because

Conclusions References
Tables Figures

Back Close
Full the comparisons of Greenland and Antarctic lakes (Banwell et al., 2014) are based on relative depths, their conclusions are likely still valid.Arnold et al. (2014) believed that their model under-predicted water depths, which could in reality mean that their model is behaving correctly but their validation data (i.e.Landsat lake depths) were biased.Shallower lake depths would also likely reduce the identified threshold for initiation of hydrofracture -if indeed lake volume initiates drainage as opposed to other basal hydrological factors (i.e., reduction of local compressional stresses, Stevens et al., 2015).
Using the same process as for Landsat sensors, we calculated g values for ASTER, MODIS, and WorldView-2 bands (see Table S1).Sneed andHamilton (2007, 2011) used ASTER's green band for lake depth estimation (g = 0.1180 m −1 ).This is ∼ 20 % smaller than the regressed value of 0.15 m −1 (see Table S1).They will therefore have likely underestimated lake depth.In all three studies discussed in the previous paragraphs, the regressed g values are much closer to the updated lab-based g values than those used previously used in the literature (see Sect. 2.1).Adoption of the new g values presented here in Tables 1 and     S1 would therefore likely lead to improved lake depth estimates.
In addition, historical Landsat imagery has been used to study the spatial evolution of melt ponds over the past decade (Howat et al., 2013); no revision of these conclusions are relevant, but improved confidence in lake depth estimates now warrants investigation of possible changes in surface meltwater storage in supraglacial lakes using the historical record.
Both supraglacial lakes and channels can contribute significantly to regional water storage an transport (Smith et al., 2015).If the water stored in supraglacial lakes in row 12 of path 008 in mid-July were spread across the whole 25 246 km 2 of ice in the scene, it would have an average depth of almost 3 cm.In other scenes, calculations provide average depths of 0.5 to 1.5 cm.Our maximum observed value is almost as high as the volume in supraglacial streams measured by Smith et al. (2015), reinforcing the potentially daily turnover of a well-connected surface system they observed.west Greenland of ∼ 2.5-3 cm per day, similar to those observed by van den Broeke et al. (2011).This implies that lakes are storing on the order of one day's worth of melt (or less), indicating daily or subdaily residence times, depending on connectivity.Lakes, therefore, provide only a transient role in water storage.So, what is their importance to the ice sheet?Lakes are projected to advance inland under a warming climate, with a minimal effect on overall ice sheet albedo but the potential to increase water transfer to the bed in areas without efficient drainage networks, therefore speeding up ice flow (Leeson et al., 2015).Indeed, episodic transport of water to the bed is expected to have a larger effect than continuously increased fluxes which set up an efficient drainage system (Flowers, 2015).In addition to episodically increased velocities, the enhanced hydrofracture propagated by the water stored in supraglacial lakes may play a key role in heat transport to the glacier bed, contributing another mechanism to increased flow.While the limit for fracture propagation on the Greenland ice sheet is at ∼ 1600 m elevation (Poinar et al., 2015), with current maximum lake depths ∼ 1200 m, this still leaves considerable room for lakes to exert their influence on the ice sheet.
Thus, although only transiently storing water, supraglacial lakes still have an important role to play in the evolving Greenland supraglacial hydrological system.
For all sensors, wavelengths, and input parameters, an important consideration for reflectance-derived lake depth is the atmospheric correction used to prepare the multispectral imagery.All imagery is processed to TOA reflectance, which means that there is some extraneous path radiance remnant in the data.Therefore, TOA values will slightly overestimate the true reflectance.This offset will not be the same between bands, and will influence the retrieved lake depths as discussed below.
The single band physically based model requires that the reflectance of optically deep water be derived for each scene separately.Effectively, this shifts the exponential decay curve of light in lake water but does not change its shape.Therefore, as long as path radiance is assumed to be homogeneous across the 185 km wide Landsat scene, TOA reflectance is sufficient for lake depth estimation.To test this assumption, the MODTRAN radiative transfer model (Berk et al., 2005)  Full radiance on a day for which Landsat data were used in northwest Greenland (18 July 2013) to investigate variations associated with variable water vapor and ozone across a Landsat scene (see Fig. S3).According to MODIS retrievals (accurate to 30 DU; Borbas et al., 2011), ozone varies within a Landsat scene on the order of approximately ±50 DU, which translates to a path radiance of ±1.6 % in the red channel.For lake depth, this can propagate on to a ∼ 20 % error in lake depth.However, because much of this error appears largely random (for a given point in time and space, although coherent patterns on a given day), although it decreases confidence in individual lake depth retrievals (and likely contributes to the ±1.6 m spread seen with WorldView validation of lake depths), averaged water volume retrieval should not be biased.For water vapor there was a 0.3 % change in path radiance between the minimum and maximum Landsat scene values, making it a small contributor to overall error.Between days, however, path radiance effects due to water vapor may vary by an order of magnitude more.
For the multiple band methods, the differential change in path radiance has larger effects.Sensitivity tests showed that a 3 % change in path radiance for one or both bands changed water volumes on the order of 10-30 %.Therefore, a more rigorous atmospheric correction is necessary in order to apply multi-band lake depth algorithms.
Still, for the study here, because validation is conducted across 12 non-consecutive days in both spring and autumn, we do not expect atmospheric conditions to bias our conclusions.
There are additional limitations to our method.As discussed above, OLI lake depth estimates (average single-band estimates from red and panchromatic bands) are robust for regional averages but not single pixels.In addition, the threshold used to identify lake extent may need to be adjusted for different regions and scenes (e.g.Banwell et al., 2014;Box and Ski, 2007).Lake depth retrievals are also sensitive to variations in ice albedo, as well as to the presence of ice lids on the surface of supraglacial lakes, which can be common in both early and late summer.Cloud cover and Landsat's 16 day revisit time also limit the conclusions that can be drawn from OLI lake depths.Many Introduction

Conclusions References
Tables Figures

Back Close
Full studies have used daily MODIS data to identify and track supraglacial lakes (e.g.Liang et al., 2012;Selmes et al., 2011;Sundal et al., 2011).Fusing the higher temporal resolution of MODIS (or additional sensors such as ESA's upcoming Sentinel-2) and higher spatial resolution of Landsat, along with more in situ calibration and validation data, should lead to unique insights in to supraglacial water storage.

Conclusion
Examination of the evolution of water storage on the surface of ice sheets and glaciers is important for understanding mass balance, dynamics, and heat transport throughout the ice mass.In this study, in situ data were used to test the capability of Landsat 8's OLI to estimate supraglacial lake depth.Promising methods were applied to two sets of Landsat observations.Patterns of water storage were similar from the two methods, but a factor of two difference was calculated for the total water volume.World-View DEMs were used to assess which of the methods was most accurate.The best method identified for Landsat 8 OLI was an average of the depth derived from singleband physically-based retrievals of Band 4 (red) and Band 8 (panchromatic); the mean difference between spectrally-derived and DEM-derived lake depths is only 0.0±1.6 m, showing no bias but some spread.Therefore, this method is recommended for future lake depth retrievals with OLI, especially for regional studies.This is the first time supraglacial lake depths have been validated across multiple dates and regions.Discrepancies between spectrally-and DEM-derived depths appear to be explained by differences between lab-measured and in situ-derived water absorption coefficients (g).The success of other sensors and bands in deriving supraglacial lake depth can thus be inferred from these g values.With this insight, multispectral lake depth estimates in the literature were revisited in the discussion.Lake extent studies can now be expanded to include lake volume with higher confidence.Updated g values are provided (see Tables 1 and S1), but further in situ data collection and satellite-based studies are needed to build more robust methods.Full The recommended depth retrieval method was applied to all available Landsat 8 imagery for summer 2014 for the Sermeq Kujalleq (Jakobshavn) region of west Greenland.Seasonal and regional trends in lake depth, evolution, and distribution were observed.Further work moving forward will need to contextualize Landsat data with other remote sensing imagery, fieldwork, and model outputs.
The Supplement related to this article is available online at doi:10.5194/tcd-9-3257-2015-supplement.Introduction

Conclusions References
Tables Figures

Back Close
Full    Full  S2).The background is elevation from the Greenland Ice Mapping Project (GIMP) DEM, courtesy BPRC Glacier Dynamics Research Group, Ohio State University (Howat et al., 2014).Introduction

Conclusions References
Tables Figures

Back Close
Full Statistics for all regressions are reported in Table 1.Introduction

Conclusions References
Tables Figures

Back Close
Full .Discrepancies in lake depth estimation for physically based retrievals can be traced to differences between lab-measured and in situ-regressed water absorption coefficients (see Table 1).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | rection (USGS Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3.2), which is then used Introduction Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 2013 (Day 215), and 274 lakes totaling 8.6 km 2 on 19 August 2013 (Day 231).Lake depths were calculated with all previously discussed methods, as well as an average between the two best single-band depth estimates.Total lake volume in the study area increased in early July, stayed al-Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Indeed, Tedesco et al. (2012) observe bare ice melt rates next to supraglacial lakes in Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .
Figure 1.Regional map showing the two study regions for lake depth estimation using Landsat 8 OLI imagery.The northwest Greenland study region is identified with a single box indicating a subscene area.The Sermeq Kujalleq (Jakobshavn) study region shows WRS-2 path/row outlines for Landsat scenes color-coded and dashed to indicate the mean latitude and average elevation of ice within the scenes (see Sect. 4.4 and TableS2).The background is elevation from the Greenland Ice Mapping Project (GIMP) DEM, courtesy BPRC Glacier Dynamics Research Group, Ohio State University(Howat et al., 2014).

Figure 2 .
Figure 2. Regression plots for in situ measured reflectance spectra used to emulate Landsat OLI and ETM+ reflectance and sonar-measured depths, including OLI single band (a), ETM+ low gain single band (b), OLI coastal and panchromatic (c), and OLI coastal and green (d).Statistics for all regressions are reported in Table1.

Figure 4 .
Figure 4. Statistics for the difference in supraglacial lake depth from physically-based and empirical methods derived from OLI imagery and WorldView DEMs, including mean/standard deviation (solid lines) and median/quartiles (dotted lines).An average of the Band 4 and Band 8 methods is used for our mapping (Figs.5 and 6).The method showing the least bias and lowest errors is an average of Band 4 (red) and Band 8 (panchromatic) single band physically based retrievals, with a mean offset of 0.0 ± 1.6 m (as indicated by the bar at the bottom of the diagram).Discrepancies in lake depth estimation for physically based retrievals can be traced to differences between lab-measured and in situ-regressed water absorption coefficients (see Table1).

Figure 5 .
Figure 5.Total water stored in supraglacial lakes over the 2014 summer using single Landsat 8 scenes covering the Sermeq Kujalleq (Jakobshavn) region (see Fig. 1, Table S1).All scenes are shown together in (a).(b) shows only the low elevation, coastal scenes, demonstrating delayed lake formation at higher latitudes.(c) shows both elevation and latitude effects in driving supraglacial water storage for scenes in WRS-2 path 8.(d) shows latitude and elevation effects for scenes in WRS-2 path 6.All sub-figures are on the same grid as part (a).
was used to simulate path

Table 1 .
Laboratory-based and in situ-derived water absorption coefficients for lake depth estimation using the physically based method (g, see Eq. 1) and empirical method (a, b, and c, see Eqs. 5 and 6).Regression statistics (correlation coefficient and root mean squared error) for lake depth estimates using field spectra convolved to emulate multispectral bands are also included.Asterisks indicate the methods applied to OLI data in this paper.Bold text indicates recommended bands for lake depth estimation with OLI.See TableS1for results from other multispectral sensors.