Automated Water Level Monitoring at the Continental Scale from ICESat-2 Photons

Abstract

Of the approximately 6700 lakes and reservoirs larger than 1 km² in the Contiguous United States (CONUS), only ~430 (~6%) are actively gaged by the United States Geological Survey (USGS) or their partners and are available for download through the National Water Information System database. Remote sensing analysis provides a means to fill in these data gaps in order to glean a better understanding of the spatiotemporal water level changes across the CONUS. This study takes advantage of two-plus years of NASA’s ICESat-2 (IS-2) ATLAS photon data (ATL03 products) in order to derive water level changes for ~6200 overlapping lakes and reservoirs (>1 km²) in the CONUS. Interactive visualizations of large spatial datasets are becoming more commonplace as data volumes for new Earth observing sensors have markedly increased in recent years. We present such a visualization created from an automated cluster computing workflow that utilizes tens of billions of ATLAS photons which derives water level changes for all of the overlapping lakes and reservoirs in the CONUS. Furthermore, users of this interactive website can download segmented and clustered IS-2 ATL03 photons for each individual waterbody so that they may run their own analysis. We examine ~19,000 IS-2 derived water level changes that are spatially and temporally coincident with water level changes from USGS gages and find high agreement with our results as compared to the in situ gage data. The mean squared error (MSE) and the mean absolute error (MAE) between these two products are 1 cm and 6 cm, respectively.

1. Introduction

Lakes and reservoirs provide an important water resource for human use. These resources range from recreational activities, power generation, drinking water, agricultural irrigation, or other commercial and industrial uses. They are also notable indicators of climatic changes [1,2] as well as local and regional drought/flooding events [3,4]. Furthermore, monitoring changes in lake and reservoir water levels are of a benefit to both local and regional water managers so that they may make more informed decisions about water management policies. This is especially true at the continental scale where the spatiotemporal changes in water levels are notably diverse across the landscape [5,6]. To this end, continental scale monitoring of lake and reservoir water level changes is particularly difficult to accomplish solely by the utilization of in situ water level gages. This statement holds true for countries that have a meaningful in situ gaging and monitoring system in place. The United States Geological Survey’s (USGS) National Water Information System (NWIS) is a notable example. To that end, Figure 1 highlights the spatial distribution of the ~6700 lakes and reservoirs >1 km² in the contiguous United States (CONUS). The USGS monitored water levels for only 430 (~6%) of those waterbodies within the past three years. In comparison, the water level monitoring technique described herein provides lake and reservoir water elevation measurements for 6163 (~92%) waterbodies. This highlights a stark need to utilize other methodologies and datasets for monitoring water level changes in the CONUS.

Remote sensing techniques have been used for many years to augment the use of in situ surface water level measurements [7,8]. For example, Ref. [9] used early spaceborne radar altimetry data from the United States Navy’s Geosat platform to monitor temporal changes in water levels within large lakes and inland seas. Furthermore, work from the early to mid-2010s utilized a multi-platform approach for longer temporal scale lake level changes [10,11,12,13] at notably larger spatial scales. More recent water level studies using radar altimeters have improved upon past research [14,15,16,17]. However, most spaceborne radar altimeters are limited in their ability to meaningfully resolve water level changes for a large quantity of lakes. This is mostly due to spatial and temporal gaps in data coverage as well as the ground footprint size of the altimeter’s energy source. The latter is especially true for radar altimeter platforms (e.g., Topex/Poseidon’s ~1 km footprint). Spaceborne laser altimeters (like NASA’s first ICESat mission) have notably smaller footprint diameters (~70 m) in comparison to their radar altimeter counterparts. The reduced footprint size of the altimeter allows for water level changes to be derived for waterbodies with smaller areal extents. This of course increases the number of lakes and reservoirs where meaningful measurements can be acquired and allows for a more complete picture of surface water changes. Several researchers have successfully utilized highly accurate ICESat (IS-1) laser altimetry products to monitor water level changes over the period of the sensor’s lifetime [18,19,20,21].

In 2018, NASA launched IS-1’s successor, ICESat-2 (IS-2) into polar orbit. Like its predecessor, IS-2’s main scientific objectives revolve around cryospheric measurements in the polar regions [22]. However, secondary mission objectives do involve the monitoring of inland surface water height changes. The Advanced Topographic Laser Altimeter System (ATLAS) is the primary sensor onboard the IS-2 platform and is a notable improvement over IS-1’s Geoscience Laser Altimeter System (GLAS). The ATLAS sensor utilizes three different pairs of beam tracks (six in total), which enables an increased spatial coverage as compared to the single track of the GLAS sensor onboard the IS-1 mission. Furthermore, the ATLAS sensor utilizes a novel photon counting approach that allows for an increase in both the precision and the accuracy of the vertical time-of-flight elevation measurements [23]. The ground footprint and the laser posting for the IS-2 platform have both been vastly improved from the GLAS sensor as well. The ATLAS sensor’s footprint is ~17.5 m with an along-track posting of about 70 cm [22]. The specifications of the ATLAS sensor allow it to resolve more and smaller waterbodies than its predecessor. This is especially true when comparing it to spaceborne radar altimeter platforms with their notably larger footprints.

Recently, more water level analyses have begun to utilize IS-2 ATLAS datasets into their studies. For example [24,25], employ ATLAS products to accurately monitor water level changes for several lakes on the Tibetan Plateau (TP) and for ~220,000 global waterbodies, respectively. A recent study also notes that utilizing IS-2 data increases the quantity of measurable lakes on the TP by a factor of two as compared to IS-1 datasets [24]. An accuracy comparison of in situ water level gage readings with levels derived from IS-2 and a modern spaceborne altimeter (Satellite with ARgos and ALtika [SARAL]) for around 30 reservoirs in China shows that the relative altimetric accuracy from IS-2 data is nearly two decimeters better than SARAL’s [26]. Some studies have compared IS-2 water level measurements with in situ gage data and have found high relative accuracies [26,27,28,29]. These studies further highlight the quality of lake and reservoir monitoring products from IS-2 data as compared to previous spaceborne altimeters. However, there are some limitations to consider when utilizing IS-2 data to monitor water levels. In particular, there is a known issue that occurs over smooth open water where false multiple surface returns are found in the raw ATL03 data. These double echoes are likely caused by after-pulses or electronic noise after the primary surface return [30]. The workflow proposed herein helps to alleviate this issue.

This study seeks to provide a novel automated workflow that utilizes the latest spaceborne altimetric products in order to monitor lake and reservoir water level changes for all waterbodies >1 km² in the CONUS. The vast number of waterbodies in the CONUS requires the need for the automated workflow aspect of this work. Furthermore, users of the readily available surface water height products derived from IS-2’s ATLAS sensor (i.e., ATL13) are reliant on their static water body extents that are built into their processing pipeline. This necessitates the need for the novel automated workflow proposed herein. An added objective of this work is to provide accuracy assessments of these remotely sensed water level products as compared to thousands of temporally and spatially overlapping water level changes from USGS gage readings. Furthermore, another goal of this study is to disseminate these results via an interactive website where an interested reader may better understand the spatiotemporal differences in lake and reservoir level changes at this scale. Lastly, this work seeks to provide these IS-2-based water level products to other users in the community in hopes that they will be of a benefit for future studies. This research was undertaken in order to provide baseline water level monitoring on a much larger scale than is traditionally possible with in situ monitoring. Furthermore, knowledge of the accuracy of these remotely sensed spaceborne laser altimeter water level products validates the quality of the IS-2 ATLAS platform for use in future water level monitoring studies. This work presents the largest (to the best of our knowledge) comparison of IS-2 water level measurements to in situ gage measurements. We utilize careful processing techniques and cluster computing workflows with a validated waterbody extent product and all spatially overlapping ATL03 photons for lakes and reservoirs >1 km² in the CONUS along with Landsat Dynamic Surface Water Extent (DSWE) products in order to satisfy our goals as laid out above.

2. Data and Methods

2.1. ICESat-2 ATLAS ATL03

This study uses photons from the IS-2 ATLAS platform in order to derive water level changes for all overlapping lakes >1 km² in the CONUS. The spaceborne laser altimeter platform was launched in 2018 and the first available products were acquired on 12 October of the same year. The ATLAS sensor uses a beam splitting approach to separate photons into three different pairs of beams (six in total). These individual beams have been named GT1L/GT1R, GT2L/GT2R, and GT3L/GT3R. The beams in each pair (e.g., GT1L/GT1R) are separated by ~90 m in the across-track direction and by ~2.5 km in the along-track direction and each of the three beam pairs are separated by ~3.3 km. Each pair consists of a “strong” beam and a “weak” beam where the “strong” beams have ~4 times the energy as the “weak” beams [22]. The ground footprint size and along-track spacing of the photons is ~17.5 m and ~70 cm, respectively.

Photon data from the ATLAS sensor are provided to the scientific community in different data products at varying levels of post-processing. There is an inland water surface height Level-3A product (ATL13) derived by the IS-2 science team that provides along-track water surface heights from overlapping lakes, reservoirs, bays, estuaries, and rivers [31]. These ATL13 data are derived from a Level-2 global geolocated photon data product (ATL03). The ATL13 dataset consists of estimated mean water surface heights within different segment lengths (~100 m and ~1–3 km). These segmented heights were derived from individual ATL03 photons that are within an inland body mask used in their processing workflow [32]. Users of this ATL13 product are reliant on the inland waterbody mask used in that workflow for the actual bodies of water that can be investigated using the ATL13 surface water heights. Previous versions of the ATL13 product did not allow for a complete analysis of all lakes >1 km² in the CONUS due to their selection of that particular waterbody mask.

This study utilizes individual photons from the more robust ATL03 Level-2 product so that we have complete control over which photons will be processed in our water level workflows. The ATL03 products for this study were acquired from the National Snow and Ice Data Center Distributed Active Archive Center (NSIDC DAAC) [33]. ATL03 photons are separated into 14 different data granules that encompass around 1/14 of a full IS-2 orbit. We acquired all ATL03 (Version 3) data granules overlapping the CONUS (Regions 01 and 02: ascending; Regions 06 and 07: descending) from the NSIDC servers. These consisted of ~7900 different ATL03 granules with a temporal range of October 2018 to November 2020. The ATL03 products consist of processed geolocated photons with the main fields of interest being the height above the ellipsoid (WGS84), time of photon event, various confidence flags, and geodetic latitude and longitude for the individual photons [34]. The IS-2 processing methodology employed for this research is further described in Section 2.3.

2.2. Waterbody Extents

This study uses spatial extents of all lakes and reservoirs larger than 1 km² in the CONUS as the initial waterbody mask for the ATL03 photons. The mask was derived from techniques described in [6] and underwent significant human-aided quality assurance and quality control in order to increase the quality of the waterbody extent product. The extents from this initial water mask were buffered inward by 30 m in order to reduce the quantity of edge photons (i.e., photons that capture elevation from the water/land periphery). This in-buffered extent file consists of 6690 different waterbodies and the centroids from those extents are displayed in Figure 1. Using a custom lake mask allows complete control over the size and location of waterbodies used as an initial filter for the ATL03 photons. This is in comparison to the segmented photons of the ATL13 product that rely on their own proprietary waterbody masks. Our custom waterbody mask is used in the water level workflow as described in Section 2.3 to initially filter ATL03 photons within full IS-2 granules into their spatially overlapping water boundaries.

The dynamics of surface water hydrologic systems (e.g., lakes and reservoirs) necessitate a need to further spatially filter ATL03 photons by waterbody extents that are as temporally coincident to each IS-2 data acquisition date and time as is possible. The waterbody mask described in the previous paragraph can be thought of as an initial static extent filter for lakes >1 km² in the CONUS. Landsat Level-3 DSWE datasets are used in this study to further spatially filter ATL03 photons into their actual active bodies of water. DSWE data are a gridded Landsat Level 3 science product derived from Landsat Analysis Ready Data (ARD) that provides cell-by-cell information pertaining to the existence and condition of surface water extents [35]. Landsat ARD products consist of the most geometrically accurate Landsat 4–5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI)/ Thermal Infrared Sensor (TIRS) data that are processed to the highest scientific standards and level of processing required for direct use in assessing and monitoring landscape changes [36]. These DSWE products are used in the water level workflow as described in Section 2.3 to better spatially filter ATL03 photons into the actual water boundaries of the temporally closest DSWE grid.

2.3. Altimeter Water Levels

Each of the 7898 ATL03 granules acquired from the NSIDC that overlap the CONUS were used as the basis of the IS-2 ATL03 processing workflow described in this subsection (Figure 2). A parallel computing and spatial indexing approach was used to spatially filter the IS-2 granules due to the vast quantities of ATL03 photons as well as the number of vertices in the waterbody mask. To this end, each granule was sent to its own core within UCLA’s Hoffman2 Cluster in an embarrassingly parallel computing workflow in order to split the granule photons into their respective spatially overlapping waterbody extents using the lake/reservoir mask as described in the first portion of Section 2.2. To speed up processing times, the initial static water mask is clipped based on the buffered extent of the granule that is being processed prior to the assignment of the photons into their corresponding overlapping waterbody. This step uses spatial indexing techniques to rapidly determine which of the waterbodies in the mask are likely to contain the photons in the granule. This tremendously speeds up the precise spatial determination of the intersecting photons into their respective waterbodies. This initial photon intersection analysis is done for each of the six tracks (three “weak” and three “strong”) for the granule and only photons with land, inland ice, and inland water “signal_conf_ph” values flagged as high confidence (i.e., a value of 4) are saved for further processing. The outputs of this initial processing step are arrays of high confidence photons from all of the 7898 granules for every track for each waterbody in the inward-buffered mask.

The extents of several waterbodies in the lake/reservoir database exhibit notable inter- and intra-annual variations in areal extent. However, the initial waterbody mask provides only a snapshot of these dynamic areal extents. That said, care must be taken in order to filter out ATL03 photons that made it through the initial static spatial filtering process but do not actually fall within regions in which surface water was present during the time of the individual IS-2 granule acquisitions. These photons should be further filtered using a dynamic waterbody extent mask in order to remove the false surface water photons from the database so that the final water level products will have an increased accuracy. Each waterbody that has photons within its extent after the initial static spatial filter was then passed to its own core for the dynamic spatial filtering process. This parallel filtering approach dramatically speeds up processing times. The individual ATL03 track photons for every date for a given waterbody were looped through and the temporally closest DSWE scene where less than 20% of that track’s photons fell within the “cloud, cloud shadow, and snow” flagged regions were used to filter out photons that did not fall within the DSWE “Water–high confidence” and “Water–moderate confidence” flagged regions. The DSWE scene with the lowest percentage of cloudy photons was utilized for tracks where the 20% cloud threshold was not met by any DSWE scene acquired within ±4 years of that particular track’s acquisition date. This dynamic spatial extent filtering process was utilized in order to remove as many of the false surface water photons over waterbodies that see larger spatial extent changes.

The individual tracks were then processed into “weak” beam and “strong” beam water levels for each date for every lake or reservoir in the database using an outlier filtering, segmenting, and clustering technique. To that end, photons in a given waterbody’s track were converted from WGS84 ellipsoidal height into orthometric height using the EGM2008 geoid model. Next, outlier photons were filtered out using the “dem_h” data within the ATL03 product where all photons outside of the “dem_h” mean value minus 200 m and the “dem_h” mean value plus 100 m elevation range were excluded from further processing. These photons were then histogrammed into 1 m bins and the water level for the most-frequent bin was attained. This max bin level was used to further filter out photons where only those that were within +3 m and −2 m of the max bin level were kept for further processing.

The photons were then quantitatively and spatially segmented based on a maximum number of photons in a segment and a ~100 m distance threshold value. Segments consist of 50 photons for “strong” beams and 25 photons for “weak” beams where the photons in each segment are within ~100 m from each other. Segments containing less than 50 (“strong” beam) or 25 (“weak” beam) photons were removed from further processing. A histogram peak filter was run on each set of segmented (50 or 25) photons in order to remove false subsurface water signals (displayed in Figure 3). This peak filter entails histogramming the photons for each segment into 5 cm bins and then determining the three highest frequency bins. The bins containing less than 33% of the maximum frequency bin were removed. The water level of the maximum frequency bin of the remaining three bins was selected in order to utilize the bin that is typically associated with the actual surface reflectance (and not the false sub-surface photons). However, the water level of the second most frequent bin was selected if it was >55 cm than the most frequent bin. This was done in case the segment’s photons exhibited a higher number of pulses within a lower return bin as compared to the upper return bin. Only photons in this segment that were within ±50 cm of the water level from the selected maximum frequency bin were kept for additional processing. A subsequent filtering step was applied on the segments where the photons with an absolute deviation outside of the median absolute deviation (MAD) were removed from further processing. The final segment water levels were assigned by taking the mean of the remaining photons within each segment.

These segments were then clustered by along-track distance and water level using a density-based spatial clustering of applications with noise (DBSCAN) method with an epsilon value of 50 and a minimum segment sample setting of 1 [37]. The epsilon value is analogous to the maximum distance between two segments for one segment to be considered within the neighborhood of the other segment. Prior to clustering, the segments’ latitudes were normalized such that an along-track distance of ~500 m was akin to a ~50 in the data. This allowed the epsilon value of 50 to appropriately cluster the segments by their spatial distance values of 50 cm (water level) and ~500 m (along-track distance). Clusters consisting of only one segment were removed from the processing workflow (minimum sample = 1). The mean cluster water levels were determined from the individual segments within each cluster. Clusters whose means that were outside of two standard deviations from the all-segment mean and outlier clusters that bring the all-cluster standard deviation above 20 cm were removed from further processing. Clusters with a mean absolute deviation greater than 2.5 cm were further filtered such that segments in a cluster whose water levels were outside of ±5 cm from the water level at the peak of the Gaussian-smoothed non-parametric kernel density estimation derived from the intra-cluster segment water levels were removed, and a new cluster mean was derived. Otherwise, the previous cluster mean was used. The final lake level for the track was determined by the median of all remaining cluster water level values for each beam type. Figure 3 highlights the output of this workflow and displays photons, segments, and clusters for two example “strong” beam tracks.

2.4. In Situ Gage Level and Altimeter Water Level Comparison

USGS surface water gage locations and their IDs were acquired from the USGS NWIS for twenty parameter codes related to surface water levels (i.e., 61055, 62615, 00062, 72292, 99064, 72277, 62600, 72293, 62614, 99065, 72264, 99020, 62617, 00065, 72214, 72020, 30211, 62616, 30207, 72275). These 765 different gages were manually matched to their corresponding lake or reservoir in the initial static waterbody mask that was described at the beginning of Section 2.2. An automated spatial matching approach was considered, but a manual approach was eventually decided upon. This was mostly due to the occasional poor accuracy of the reported gage coordinates used as an input into the minimum distance based spatial matching approach. The 351 different centroids for each lake in our custom static waterbody mask that have both IS-2 and gage water levels are plotted in light blue within Figure 1. All available temporally coincident in situ water level readings for each waterbody with an active water level gage were acquired from the NWIS via the hydrofunctions python package from Martin Roberge and contributors (https://github.com/mroberge/hydrofunctions, accessed on 21 March 2021). Initial analysis of the gage levels sought to convert the many different water height datums to a uniform vertical datum to directly compare to the spatially and temporally overlapping water levels derived from the IS-2 processing workflow. However, this method was abandoned due to the poor quality of many of the gages’ metadata. Most of the gages’ vertical datum offsets were not precise or accurate enough to directly compare to the IS-2 water levels. Instead, we compared temporally coincident relative water level changes from the gages and the IS-2 water level products.

Relative water level changes were derived for each of the gages’ temporally overlapping IS-2 date pairs. The IS-2 acquisition times were converted to local gage times in order to acquire the proper coincident gage readings. Furthermore, only gage readings acquired within ±24 h from their corresponding IS-2 granule’s acquisition time were utilized in the comparison. USGS gage data are served via two temporal scales: (1) daily values (DV) and instantaneous values (IV). Typically, DVs for a given gage are derived using some statistical analysis (e.g., daily mean or median) of IV values, or they are simply an IV at a given time of day (e.g., noon). All gage/IS-2 comparisons were done with the temporally closest IV data as compared to the IS-2 acquisition time. This allowed for the utilization of the temporally closest gage readings to be compared to their IS-2 water level counterparts so as to more accurately compare the water level readings between the two methods. These relative water level changes from temporally overlapping USGS gage data were ultimately used to compare to the IS-2 relative water level changes in order to assess the accuracy of our ATL03-derived water level changes as described in the workflow in the previous subsection.

3. Results

3.1. Altimeter Water Levels

Water levels were successfully derived using the workflow described in Section 2.3 for 6163 different lakes within the CONUS. Approximately 84% (5177) of those lakes had more than a single day’s IS-2 water level measurement allowing for a time series of water level changes. It is not practical to display to the reader the entirety of our water level results from this study in static graphs. However, the results from all ~6000 lakes and reservoirs examined in this study are viewable via an online interactive website where an interested reader may better understand the spatiotemporal differences in water level changes at this scale (https://icesat2waterlevelmapping.s3-us-west-1.amazonaws.com/overViewMap.html, accessed on 10 September 2021). We invite all interested readers to explore the website. Furthermore, all water level segments and clusters derived by this research can be downloaded from the above website. The link provides an overview map showing the location of all 6163 lakes in this study. Selecting an individual lake allows the user to view the interactive water level change time series separated by track and beam type as well as the interactive daily data plots for each date in the time series. These daily plots display the ATL03 photons, the segments, and their clusters as described in Section 2.3 and are separated by track and beam type for a single date for the selected lake. Furthermore, each individual daily date pages links to an interactive map where the beam segments and clusters for this date are overlaid on a true-color basemap.

Two static examples of the daily data plots are provided in Figure 3 for a “strong” beam of a single track for two separate waterbodies. This figure shows two example tracks broken down by their filtered photons, segmented photons, and clustered segments. We also provide sample water level time series from four reservoirs and four lakes examined in this study along with their in situ gage water levels in Figure 4. The location of the waterbodies plotted in Figure 3 and Figure 4 are denoted within Figure 1 (white circles). The complete IS-2 derived water level results from our processing workflow are available at the previously mentioned link. An examination of lake IDs 1742 and 3370 within Figure 4 provides evidence that our processing methodology is capable of resolving water level changes at least at the decimeter scale.

3.2. In Situ Gage Level and Altimeter Water Level Comparison

There are 18,690 temporally overlapping in situ and IS-2 derived (IV gage readings and “strong” beam only) relative water level change measurements from 708 unique dates and 330 different gaged lakes in our database. A scatter plot showing the relative water level measurement comparisons for all ~19,000 of these points is plotted in Figure 5 (R² = 0.99) with a notable positive correlation. Two zoomed insets (−1 to +1 m and −0.25 to +0.25 m) are provided in order to show a more focused comparison centered on the origin. The IS-2 water level changes were less than their in situ gages’ counterparts in roughly 49.4% of the gage and “strong” beam IS-2 comparisons while around 50.6% were more than the gaged water level measurements. The mean squared error (MSE) and the mean absolute error (MAE) between these two products are 1 cm and 6 cm, respectively. The mean (median) absolute residual difference between the two products is ~6 cm (~4 cm) with a standard deviation of ~8 cm. Furthermore, these residuals are similar to or better than other IS-2 water level studies (e.g., [28]). This provides further evidence about the quality of our ATLAS derived water level products described herein.

We plot a frequency histogram of those differences in Figure 6 in order to better examine the spread of the relative water residuals. This plot shows that ~97% (18,214) of the 18,690 temporally overlapping “strong” beam IS-2 derived water level change measurements are within ±25 cm of their gaged (IV) counterparts. Furthermore, ~81% (15,085) and ~55% (10,372) of the IS-2 water level changes are within ±10 and ±5 cm of their corresponding in situ gaged measurements, respectively. It is evident that the vast majority of our IS-2 derived water level changes are quite accurate as compared to their temporally overlapping in situ gaged counterparts. These results are proof positive that IS-2 ATLAS photons can meaningfully monitor water level changes at the continental scale for a variety of surface waterbody types.

4. Discussion

4.1. Altimeter Water Levels

Water levels were successfully determined for 6163 lakes within the CONUS using the methodology described in Section 2.3. We note that this methodology uses the global geolocated ATL03 photons as the main IS-2 input and that these photons are the basis for higher-level IS-2 products (e.g., ATL08 and ATL13). The ATL03 photons were utilized in this study as it allows for complete control over the associated waterbody mask used as the initial photon spatial filter as well as for control over segmenting and clustering of the photons. However, this does come with drawbacks as the sheer volume of ATL03 data is much larger than the post-processed and segmented higher-level products (e.g., ATL08). These higher-level IS-2 products would be less computationally expensive to download and process as compared to the ATL03 data. However, the user would then be reliant on their post-processing segmenting algorithms (e.g., ATL08 or ATL13) as well as their waterbody mask (e.g., ATL13). We note that some researchers have had success using the ATL08 products for water level studies [28]. Using the ATL08 product still allows for the utilization of a custom waterbody mask. However, the user is dependent on the segmented photons provided in that product. To this end, we have compared segments derived from our workflow with segments from the ATL13 and ATL08 science products for an individual track from an example lake and reservoir. These comparisons are plotted in Figure S1 and they highlight a few of these points.

It is evident in the top plot that the ATL13 segments include photon returns that are not over active surface water. This is attributed to the fact the ATL13 workflow uses a static (as opposed to a static and a dynamic mask as is used within this study) surface water mask and the water level on this particular IS-2 overpass was lower than when the static water mask was generated (or the ATL13 mask is inaccurate in this region). This further highlights how the methodology within this research can improve upon ATLAS-derived water level monitoring products. The second and third plots within Figure S1 show the similarity between the ATL03 segments as derived from our workflow as compared to the actual surface water ATL13 segments and the ATL08 segments. These similarities are proof positive that our workflow is capable of producing meaningful water level segments. A closer look at the bottom plot in Figure S1 shows that there is a region of missing ATL13 segments. We note that our ATL03 derived segments (as well as the ATL08 segments) pick up this surface water area, but the ATL13 processing workflow does not. We attribute this ATL13 data issue to errors in their processing workflow or water mask.

We note that water level errors can occur during the second spatial filtering step in our workflow where the temporally closest DSWE product is used to further filter out ATL03 photons that do not fall within areas actively classified as surface water. These errors can appear for waterbodies that see large fluxes in areal extent between IS-2 data acquisitions and their temporally closest DSWE acquisition (e.g., some reservoirs and ephemeral lakes). In particular, this can happen when ATL03 photons fall on regions of inactive surface water within the initial waterbody mask and those regions subsequently fill in with water during the time between the IS-2 photon acquisition and the temporally closest DSWE acquisition. In those cases, the DSWE will not filter out the photons because the workflow believes that they are over active water regions. Future work involves increasing the number of remotely sensed datasets (e.g., Planet and Sentinel) used in the secondary spatial filtering step in order to reduce the temporal lag between the IS-2 data acquisition and the acquisition of the satellite derived active surface water data. This will help to alleviate the aforementioned issue.

The benefits of our photon filtering, segmenting, and clustering technique are highlighted in Figure 3. These two different examples show how this methodology is able to determine the individual track water level by automatically removing outlier photons and appropriately grouping the remaining photons into meaningful segments and clusters. For example, the bottom plot shows the along-track water levels for a single “strong” beam acquired over a reservoir in Georgia. The data gaps present in the along-track direction of the IS-2 path for this track denote regions where the IS-2 overpass goes through several “fingers” of the reservoir where the pulses received at the sensor alternate between water photons and land photons. These types of overpasses are more complicated than a single continuous track through water (as in the top plot within the same figure). This alternating land–water overpass example highlights the robustness of our segmenting and clustering methodology. It is also evident in both plots within Figure 3 that there are multiple surface returns apparent around −2 m and −4 m from the actual water level. These double echoes are a known issue that often occurs over smooth open water surfaces and is likely caused by after-pulses or electronic noise after the primary surface return [30]. These two plots within Figure 3 further highlight the robustness of our IS-2 water level workflow.

Our workflow separates water level change products into both “strong” and “weak” beams. We compared the 43,216 different IS-2 water level measurements from both the overlapping “strong” and “weak” beams and found that the mean (median) absolute residual difference between the “strong” and “weak” beam water levels is ~4 cm (~3 cm) with a standard deviation of ~10 cm. Figure 7 shows a frequency histogram of the spread of “strong” and “weak” beam water level residuals. This plot shows that ~93% (40,372) of the 43,216 overlapping “strong” and “weak” beam IS-2 derived water level measurements are within ±10 cm of each other. Furthermore, ~48% (20,820) and ~20% (8796) of these “strong” and “weak” beam IS-2 water levels are within ±2.5 and ±1 cm of each other, respectively. The “strong” beam water levels are less than their “weak” beam water level counterparts in around 69% (29,708) of these beam comparisons. The majority of the “strong” beam water levels are likely lower than their “weak” beam counterparts due to the increased energy (~4x) of the “strong” photons as compared to the “weak” photons. This increased energy could allow for more penetration into the water’s subsurface and slightly decreases (or pulls down) the water levels for the “strong” beams. Regardless, this comparison provides robust evidence for the utilization of both beam types for water level change analyses, assuming that the potential several centimeter difference between the two beams is acceptable for the particular use case.

4.2. In Situ Gage Level and Altimeter Water Level Comparison

Figure 5 and Figure 6 show a scatter plot of the relative water level changes from in situ gages (IV) and IS-2 derived water levels from “strong” beams only as well as a frequency histogram of those relative water level changes, respectively. We compared the relative water level changes instead of the individual water levels due to the fact that the absolute water level readings from many of the gages were not precise enough to directly compare to there IS-2 derived water level counterparts. For example, many gage readings are given in a local vertical datum, and in the best cases, the datum offsets are provided in the gages’ metadata. These offsets would allow those measurements to be directly comparable to the IS-2 derived water levels that are in a known vertical datum. However, we find that even in the gages that do provide vertical datum offsets the accuracy of those offsets is not known or is not precise enough for a trustworthy IS-2 comparison. Therefore, we have chosen to compare the relative water level differences between the temporally overlapping gage measurements and our IS-2 water level change results. These gage datum differences are evident in Figure 4 where the IS-2 water level and gage water level time series are plotted for eight bodies of water. The second y-axis (right) shows the gage water level, and it is evident that the gage datums are offset from the orthometric heights of the IS-2 time series. We expect these two y-axes to be offset as the gage datums have not been properly converted due to the issues mentioned at the beginning of this paragraph. Using the relative water level differences between two coincident dates allows us to disregard this datum-offset issue and still provides for a meaningful and direct accuracy assessment between the gage derived and the IS-2 derived water levels.

Furthermore, in most cases, we do not expect the water level changes from the IS-2 data and the in situ gage data to be the same. The gage readings are merely a water level snapshot in time at a single location on the waterbody’s surface. Lake and reservoir water levels are highly dynamic across the water’s surface, and this dynamism leads to uneven water elevations for the waterbody. To that end, along-track IS-2 photons are collecting surface water return photons through a cross-section of the waterbody, and these surface water photons will have different water elevations as caused by the multitude of factors that allow for the uneven surface water levels. These factors include wind speed and direction, gravity effects, seiches, pressure oscillations, as well as lake connectivity and are the cause for uneven water levels and slopes across a given body of water. A waterbody is not expected to have a uniformly flat surface and a water level measurement at one location within a lake or a reservoir is not expected to be the same measurement at a different location within the same body of water. Therefore, we expect there to be some differences in the IS-2 water levels and the gage water levels presented in Figure 4, Figure 5 and Figure 6.

The gage and IS-2 comparisons shown in Figure 5 and Figure 6 and discussed in the preceding paragraphs as well as in Section 3.2 utilize IS-2 water level changes from ATLAS “strong” beams only. In comparison to those results, there are 16,525 temporally overlapping in situ gage and IS-2 derived (IV gage readings and “weak” beam only) relative water level change measurements from 700 unique dates and 328 different gaged lakes in our database. A scatter plot showing the relative water level measurement comparisons for all ~17,000 of these “weak” points is plotted in Figure S2 (R² = 0.99). Similar to the “strong” gage comparisons, the IS-2 water level changes were less than their gaged counterparts for roughly 49.9% of the gage and “weak” beam comparisons, while around 50.1% were more than the gaged water level measurements. The MSE and the MAE between these two products are 1 cm and 7 cm, respectively. The mean (median) absolute residual difference between the two products is ~7 cm (~5 cm) with a standard deviation of ~8 cm. These results are around 1 cm greater than their “strong” beam gage comparison counterparts.

Figure S3 shows a frequency histogram of the spread of relative water level change residuals. This plot shows that ~97% (16,083) of the 16,525 temporally overlapping “weak” beam IS-2 derived water level change measurements are within ±25 cm of their gaged counterparts. Furthermore, ~78% (12,891) and ~52% (8640) of the IS-2 water level changes are within ±10 and ±5 cm of their corresponding gaged measurements, respectively. This comparison between the “weak” beam and gage derived relative water level changes as well as the IS-2 derived “strong” and “weak” beam water level residual comparisons in Figure 7 show that the methodology presented herein is capable of providing accurate water level changes from both types of ATLAS beams. This is important as it means that both “strong” and “weak” beams can be utilized to derive water level changes, which implies that more waterbodies can be investigated and that more water level data points can be derived as compared to a single “strong” beam analysis only.

5. Conclusions

This research has highlighted an automated cluster computing workflow that utilizes IS-2 ATLAS ATL03 photons to accurately determine water level changes for all waterbodies >1 km² in the CONUS. To this end, water levels were successfully derived for 6163 different lakes and reservoirs using the methodology described herein. A comparison of relative water level changes between in situ gage (IV) measurements and our IS-2 water level changes showed an absolute mean (median) residual of ~6 cm (~4 cm) with a standard deviation of ~8 cm. The MSE and the MAE between these two products were 1 cm and 6 cm, respectively. Around 97% of the temporally overlapping “strong” beam IS-2 derived water level change measurements were within ±25 cm of their gaged counterparts. Furthermore, ~81% and ~55% of these IS-2 water level changes were within ±10 and ±5 cm of their corresponding in situ gaged measurements. A comparison of the “strong” and “weak” beam water level measurements showed that the mean (median) absolute residual differences between the two beams was ~4 cm (~2 cm) and had a standard deviation of ~4 cm. The MSE and the MAE between these two products were 1 cm and 7 cm, respectively. Around 93% of the “strong” and “weak” beam IS-2 derived water level measurements were within ±10 cm of each other. An interactive website (https://icesat2waterlevelmapping.s3-us-west-1.amazonaws.com/overViewMap.html, accessed on 10 September 2021) was created so that interested readers may better understand the spatiotemporal differences in lake and reservoir level changes at this scale. Future work entails increasing the number of spaceborne platforms (e.g., Planet and Sentinel) used during the dynamic secondary spatial filtering step, and expanding the spatial extent of the study region from the CONUS to the global scale. This research was undertaken in order to provide baseline water level monitoring on a much larger scale than is traditionally possible with in situ monitoring. Furthermore, knowledge of the accuracy of these remotely sensed spaceborne laser altimeter water level products validates the quality of the IS-2 ATLAS platform for use in future water level monitoring studies.