Validation of ICESat-2 ATLAS Bathymetry and Analysis of ATLAS’s Bathymetric Mapping Performance

NASA’s Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2) was launched in September, 2018. The satellite carries a single instrument, ATLAS (Advanced Topographic Laser Altimeter System), a green wavelength, photon-counting lidar, enabling global measurement and monitoring of elevation with a primary focus on the cryosphere. Although bathymetric mapping was not one of the design goals for ATLAS, pre-launch work by our research team showed the potential to map bathymetry with ICESat-2, using data from MABEL (Multiple Altimeter Beam Experimental Lidar), NASA’s high-altitude airborne ATLAS emulator, and adapting the laser-radar equation for ATLAS specific parameters. However, many of the sensor variables were only approximations, which limited a full assessment of the bathymetric mapping capabilities of ICESat-2 during pre-launch studies. Following the successful launch, preliminary analyses of the geolocated photon returns have been conducted for a number of coastal sites, revealing several salient examples of seafloor detection in water depths of up to ~40 m. The geolocated seafloor photon returns cannot be taken as bathymetric measurements, however, since the algorithm used to generate them is not designed to account for the refraction that occurs at the air–water interface or the corresponding change in the speed of light in the water column. This paper presents the first early on-orbit validation of ICESat-2 bathymetry and quantification of the bathymetric mapping performance of ATLAS using data acquired over St. Thomas, U.S. Virgin Islands. A refraction correction, developed and tested in this work, is applied, after which the ICESat-2 bathymetry is compared against high-accuracy airborne topo-bathymetric lidar reference data collected by the U.S. Geological Survey (USGS) and the National Oceanic and Atmospheric Administration (NOAA). The results show agreement to within 0.43—0.60 m root mean square error (RMSE) over 1 m grid resolution for these early on-orbit data. Refraction-corrected bottom return photons are then inspected for four coastal locations around the globe in relation to Visible Infrared Imaging Radiometer Suite (VIIRS) Kd(490) data to empirically determine the maximum depth mapping capability of ATLAS as a function of water clarity. It is demonstrated that ATLAS has a maximum depth mapping capability of nearly 1 Secchi in depth for water depths up to 38 m and Kd(490) in the range of 0.05–0.12 m−1. Collectively, these results indicate the great potential for bathymetric mapping with ICESat-2, offering a promising new tool to assist in filling the global void in nearshore bathymetry.


Introduction
Nearshore bathymetry, which we define loosely as submerged topography spanning the 0-10 m depth range, is notoriously difficult to map. Conventional, terrestrial surveying technologies, such as total stations and global navigation satellite system (GNSS) receivers, do not work underwater. Hence, these technologies are only viable in the shallowest waters, in which a person can stand with a survey rod, and, even then, such methods can be challenging, dangerous (especially in high-energy nearshore environments), and slow. Meanwhile, boat-based acoustic hydrographic surveying technologies, such as multibeam echosounders (MBES), enable efficient and accurate data acquisition in deeper waters (>4-5 m). However, ships and even small boats can be dangerous to operate in the shallowest areas, especially in the presence of breaking waves and submerged hazards, such as corals and rocks. These hazards prompted the National Oceanic and Atmospheric Administration (NOAA) hydrographic

Introduction
Nearshore bathymetry, which we define loosely as submerged topography spanning the 0-10 m depth range, is notoriously difficult to map. Conventional, terrestrial surveying technologies, such as total stations and global navigation satellite system (GNSS) receivers, do not work underwater. Hence, these technologies are only viable in the shallowest waters, in which a person can stand with a survey rod, and, even then, such methods can be challenging, dangerous (especially in high-energy nearshore environments), and slow. Meanwhile, boat-based acoustic hydrographic surveying technologies, such as multibeam echosounders (MBES), enable efficient and accurate data acquisition in deeper waters (>4-5 m). However, ships and even small boats can be dangerous to operate in the shallowest areas, especially in the presence of breaking waves and submerged hazards, such as corals and rocks. These hazards prompted the National Oceanic and Atmospheric Administration (NOAA) hydrographic surveying fleet to implement a policy of not surveying shoreward of a pre-defined boundary, termed the navigable area limit line (NALL), which is generally the offshore-most of the following: a) an offset of 0.8 mm from the mean high water (MHW) shoreline at the scale of the largest-scale nautical chart of the area (e.g., 64 m for a 1:80,000 sale chart), b) the 3.5 m depth contour, or c) the inshore limit of safe navigation, due to kelp, rocks, breaking waves, etc. [1].
Unfortunately, the lack of nearshore bathymetric data hinders a number of coastal science, management, and engineering applications. Important uses of nearshore bathymetry include analysis of nearshore hydrodynamics, morphology, and sediment transport [2], tsunami inundation modeling [3], and benthic habitat mapping [4], among others [5]. Due to the importance of shallow nearshore bathymetric data and the challenges in acquiring the data with conventional surveying technologies, remote sensing methods of bathymetric mapping have been of interest since at least the 1970s [6,7], and the interest appears to be growing, based on the rapidly increasing number of publications on this topic (Figure 1). The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2), which was launched on September 15, 2018, may afford a powerful new tool for addressing the challenges in nearshore bathymetric mapping. ICESat-2 is a follow-on to the original ICESat mission, which was operational for seven years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). While the original ICESat was used in at least one notable study to estimate lake water storage capacity and reconstruct bathymetry [8], the surface mapping laser in its Geoscience The Ice, Cloud, and Land Elevation Satellite-2 (ICESat-2), which was launched on September 15, 2018, may afford a powerful new tool for addressing the challenges in nearshore bathymetric mapping. ICESat-2 is a follow-on to the original ICESat mission, which was operational for seven years (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). While the original ICESat was used in at least one notable study to estimate lake water storage capacity and reconstruct bathymetry [8], the surface mapping laser in its Geoscience Laser Altimeter System (GLAS) instrument used a near infrared wavelength, which does not penetrate water to any appreciable depth. ICESat-2 carries a single instrument, the Advanced Topographic Laser Altimeter System (ATLAS), a green (532 nm), photon-counting lidar with a 10kHz pulse repetition rate, and nominal 17 m diameter footprint [9][10][11]. ATLAS employs six beams arranged in three pairs, where each pair consists of a "weak" and "strong" beam (separated by~90 m in the across track direction) and provides a nominal 0.7 m along track point spacing from a 496 km mean orbital altitude [10]. The relative energy levels within the pair allows for improved response to variations in surface reflectivity. While nearly all methods of mapping bathymetry from satellite remote sensing are based on passive, multispectral imagery [12][13][14][15][16], ATLAS provides a more direct form of bathymetric mapping from an active sensor as the first Earth-orbiting laser capable of penetrating water at a high resolution in the along-track direction. Due to the more direct bathymetric measurement process, ATLAS-derived bathymetry is less susceptible to false readings than bathymetry retrieval from multispectral imagery (e.g., erroneous readings due to change in substrate type or other confounding variables introduced by the environment) and also does not require reference or "seed" depths. However, bathymetry from multispectral imagery does offer the advantage of providing data over large spatial extents, not limited to the discrete profiles of ATLAS. Thus, the synergistic fusion of active (specifically, ICESat-2) and passive (multispectral satellite imagery) data may ultimately provide the optimal solution for shallow nearshore bathymetric mapping by leveraging the strengths of each [17].
Because water clarity, and, hence, depth mapping performance for any bathymetric lidar tends to vary temporally in many locations, ICESat-2's repeat cycle of 91 days and the pointing strategy for distributed measurements in the mid-latitudes provide opportunities to explore seasonal dynamics and temporal characteristics for sites of interest. Previous work by the research team [17][18][19] investigated the bathymetric mapping potential of ICESat-2, using data from the Multiple Altimeter Beam Experimental Lidar (MABEL), an engineering testbed for ATLAS, and laser-radar equation-based modelling. However, prior to launch, a number of sensor parameters were only roughly known, such that precise quantification of ATLAS's bathymetric mapping capability was limited.
With ICESat-2 now on orbit and ATLAS data available to the scientific community, the unanswered questions from the project team's previous research can be addressed. Specifically, the goals of this study are to empirically address the following questions:
If so, how accurate is the ATLAS-derived bathymetry? 3.
What is ATLAS's maximum depth mapping capability, as a function of water clarity?
The answers to these questions will enable coastal researchers, coastal zone managers, and nautical charting offices worldwide to assess the utility of ATLAS for filling in prevalent shallow-water data gaps and for using the bathymetry for applications ranging from storm inundation modeling to benthic habitat mapping to nautical charting and coastal change analysis.

Refraction Correction
The ATLAS ATL03 data product provides the geolocated photon detections [20]. These data are the input to the higher-level geosphysical data products that use customized algorithms for specific surfaces (land ice, sea ice, vegetation, etc.). However, none of the higher-level products contain the required algorithmic steps to produce valid bathymetric measurements. More precisely, the ICESat-2 products do not currently account for the refraction and the corresponding change in the speed of light that occurs at the air-water interface. This produces both horizontal and vertical errors in the geolocation recorded in ATL03, resulting in locations that are deeper and further off nadir than the true measurement. The impact of refraction on the ICESat-2 collection is depicted in Figure 2 4. Angle of incidence, θ1, for each photon, obtained from the elevation of the unity pointing vector for the photon. More specifically, θ1 = π/2 -(ref_elev), where ref_elev is one of the parameters available in the ATL03 output. (Note, this is actually a simplification, which does not account for the curvature of the Earth across ATLAS's swath, as will be discussed later.) 5. Azimuth of the unit pointing vector for the photon, κ, obtained from the parameter, ref_azimuth, which is one of the output parameters from the ATL03 algorithm. The approach used in the refraction correction is able to compute horizontal and vertical offsets that can be applied to the uncorrected bottom return photon's coordinates, as illustrated in Figure 3. The advantages of this method are that it is both computationally efficient and geometrically intuitive. While the NASA ATL03 algorithm works in an inertial frame, our refraction correction is derived in a spacecraft-centered coordinate system defined as follows: Z is the vertical direction (parallel to and opposite the local direction of gravity), and Y is perpendicular to Z (i.e., in the horizontal plane) and oriented along the azimuth of the pointing vector. Note that while the Z axes are the same in Figures 2 and 3, the Y axes may differ. They are the same in the case that the azimuth of the photon pointing vector is exactly aligned with the across-track direction. If desired, it is possible to consider this as a 3D coordinate system by defining X to be mutually perpendicular to Y and Z and oriented so as to complete a right-handed coordinate system (i.e., out of the page in Figure 3). Once obtained, the horizontal offset ΔY can be projected onto the local (E, N) coordinate axes using The input to the refraction correction algorithm developed and used in this work consists of the following: 1. Geolocated seafloor photon returns, subset from the ATL03 geolocated photon returns through either manual or automated analysis to segment seafloor points; 2. A water surface model, W ij , derived from the water surface photon returns; 3. Estimates of the refractive indices of air, n 1 , and of water, n 2 ; 4. Angle of incidence, θ 1 , for each photon, obtained from the elevation of the unity pointing vector for the photon. More specifically, θ 1 = π/2 -(ref_elev), where ref_elev is one of the parameters available in the ATL03 output. (Note, this is actually a simplification, which does not account for the curvature of the Earth across ATLAS's swath, as will be discussed later.) 5. Azimuth of the unit pointing vector for the photon, κ, obtained from the parameter, ref_azimuth, which is one of the output parameters from the ATL03 algorithm.
The approach used in the refraction correction is able to compute horizontal and vertical offsets that can be applied to the uncorrected bottom return photon's coordinates, as illustrated in Figure 3. The advantages of this method are that it is both computationally efficient and geometrically intuitive. While the NASA ATL03 algorithm works in an inertial frame, our refraction correction is derived in a spacecraft-centered coordinate system defined as follows: Z is the vertical direction (parallel to and opposite the local direction of gravity), and Y is perpendicular to Z (i.e., in the horizontal plane) and oriented along the azimuth of the pointing vector. Note that while the Z axes are the same in Figures 2  and 3, the Y axes may differ. They are the same in the case that the azimuth of the photon pointing vector is exactly aligned with the across-track direction. If desired, it is possible to consider this as a 3D coordinate system by defining X to be mutually perpendicular to Y and Z and oriented so as to complete a right-handed coordinate system (i.e., out of the page in Figure 3). Once obtained, the horizontal offset ∆Y can be projected onto the local (E, N) coordinate axes using the azimuth of the unit pointing vector, κ, which is one of the output parameters (ref_azimuth) from the ATL03 algorithm [20].
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 20 the azimuth of the unit pointing vector, κ, which is one of the output parameters (ref_azimuth) from the ATL03 algorithm [20]. The refraction correction can be considered a simple rotation and scaling of the in-water range vector, but it is most easily derived and intuited geometrically with reference to Figure 3. In the figure, is the angle of incidence; is the angle of refraction; D is the uncorrected depth, which is obtained by differencing the Z coordinates of the water surface model and uncorrected bottom return photons; S is the slant range to the uncorrected bottom return photon location; R is the corrected slant range; P is the distance between the uncorrected and corrected photon return points in the Y-Z plane, and other angles and distances are as labeled. It can be readily recognized that = − , and = − . From Snell's Law, the angle of refraction, , is Due to the change in the speed of light that occurs at the air-water interface (which is not accounted for in the current ATL03 geolocation procedure), the relationship between R and S is The remaining relationships needed to compute (ΔY, ΔZ) can be obtained from the law of sines, the law of cosines, and simple trigonometric functions, referencing three triangles in the lower half of Figure 3: triangle DTS (right triangle), triangle RPS (scalene triangle), and the triangle with ΔY and ΔZ as two of the sides and P as the remaining side (right triangle). From triangle DTS, we have Applying the law of sines to triangle RPS and solving for α yields The refraction correction can be considered a simple rotation and scaling of the in-water range vector, but it is most easily derived and intuited geometrically with reference to Figure 3. In the figure, θ 1 is the angle of incidence; θ 2 is the angle of refraction; D is the uncorrected depth, which is obtained by differencing the Z coordinates of the water surface model and uncorrected bottom return photons; S is the slant range to the uncorrected bottom return photon location; R is the corrected slant range; P is the distance between the uncorrected and corrected photon return points in the Y-Z plane, and other angles and distances are as labeled. It can be readily recognized that φ = θ 1 − θ 2 , and β = γ − α. From Snell's Law, the angle of refraction, θ 2 , is Due to the change in the speed of light that occurs at the air-water interface (which is not accounted for in the current ATL03 geolocation procedure), the relationship between R and S is The remaining relationships needed to compute (∆Y, ∆Z) can be obtained from the law of sines, the law of cosines, and simple trigonometric functions, referencing three triangles in the lower half of Figure 3: triangle DTS (right triangle), triangle RPS (scalene triangle), and the triangle with ∆Y and ∆Z as two of the sides and P as the remaining side (right triangle). From triangle DTS, we have Applying the law of sines to triangle RPS and solving for α yields Applying the law of cosines to triangle RPS gives From the right triangle with ∆Y, ∆Z, and P as sides And finally: A final step is to project the horizontal offset, ∆Y, onto the (E, N) axes using the azimuth of the laser pointing vector, κ: After alignment of local East and North directions to the map projection (e.g., Universal Transverse Mercator) Easting and Northing axes (using the convergence angle), ∆E and ∆N from Equations (10)- (11) serve as the corrections to the unrefracted measurement's horizontal coordinates. Simply stated, the refraction correction algorithm initiates with the input parameters n 1 , n 2 , W ij , and for each photon, E, N, Z, θ 1 , and κ. These variables systematically apply to Equations (1)- (11) to determine the corrections (∆E, ∆N, and ∆Z) to the unrefracted photon coordinates. The corrections are added to the original coordinates (E, N, and Z) to produce a corrected position, E' = E + ∆E, N' = N + ∆N, and Z' = Z + ∆Z, where prime denotes the corrected coordinate. From [21], the index of refraction of seawater (S = 35 PSU) at atmospheric pressure, a temperature of 20 • C, and a wavelength of 540 nm is 1.34116, while the corresponding value for fresh water (S = 0) is 1.33469. These values are taken to be the default settings for n 2 in seawater and fresh water, respectively. Meanwhile, the default value of the refractive index of air, n 1 , is 1.00029.
With the center beams at near-nadir, assuming a mean orbital altitude of 496 km and 3.3 km separation of the beams on the ground [22], the incidence angle for the outer beams is only 0.38 • , and the horizontal offset for the refraction correction is 0.003D. This equates to only 9 cm at a water depth of 30 m, which can be considered negligible with a beam footprint of 17 m. However, since ATLAS has the capability for up to 5º off-nadir pointing, the horizontal offset, due to refraction, cannot, in general, be considered negligible (This study used only data from a near-nadir orientation, as the off-nadir pointing has, to date, not been tested during science mode operations). In cases where a first-order approximation to the refraction correction is acceptable and it is determined that the planimetric component of the refraction correction can be safely ignored, the corrected elevation of a seafloor photon return can be approximated as Alternately, when greater accuracy in the refraction correction is needed, a correction for the curvature of the Earth over ATLAS's swath must also be included, since Earth's curvature contributes to the nonparallelism of ATLAS's nadir direction and the Earth surface's normal at the point of intersection of the incident photon(s). This Earth curvature correction, which is applied to θ 1 , is computed as follows: where H is the orbital altitude of ICESat-2 (e.g., 496 km as a mean value), and R e is the radius of the Earth (e.g., 6371 km as a mean value).
The ATLAS refraction correction developed in this work was implemented in MATLAB and validated using data acquired by Quantum Spatial, Inc. over Lake Tahoe, California, using a high-accuracy, airborne bathymetric lidar system: the Riegl VQ-880-G. Like other airborne bathymetric lidar systems (e.g., [23,24]), and unlike ATLAS, the VQ-880-G is a linear-mode lidar system, rather than a single photon. Despite this obvious difference, the VQ-880-G data provide a nearly-ideal validation data set for two reasons: (1) The VQ-880-G processing workflow is similar to the one described above, in which subaqueous returns are first geolocated as if they were subaerial returns (i.e., as if there were no water present), and then a refraction correction is applied; and (2) Riegl VQ-880-G bathymetric lidar data have been extensively investigated [24][25][26][27], such that new refraction correction techniques can be validated by comparing the results against those obtained by Quantum Spatial. Tested on the Lake Tahoe data set, the refraction correction algorithm described above yielded root mean square (RMS) differences from those obtained by Quantum Spatial of 0.02 m vertical and 0.14 m horizontal, even with a simplifying assumption of a flat water surface.

Study Site
The area selected for the ICESat-2 bathymetry validation encompasses the waters south of St. Thomas, U.S. Virgin Islands (USVI), including the West Gregerie Channel, Crown Bay, Krum Bay, as well as Druif Bay and Flamingo Bay on the western end of Water Island (Figure 4). The site, which intersects the St. Thomas Harbor Area of Particular Concern (APC) [28], is heavily used by commercial and cruise ships, as well as recreational boat traffic [29]. Interest in the site stems from both its economic importance to the island and threats to the health of coral reef ecosystems, which provide habitats for a number of species [28,[30][31][32].
Remote Sens. 2019, 11, x FOR PEER REVIEW 7 of 20 accuracy, airborne bathymetric lidar system: the Riegl VQ-880-G. Like other airborne bathymetric lidar systems (e.g., [23,24]), and unlike ATLAS, the VQ-880-G is a linear-mode lidar system, rather than a single photon. Despite this obvious difference, the VQ-880-G data provide a nearly-ideal validation data set for two reasons: (1) The VQ-880-G processing workflow is similar to the one described above, in which subaqueous returns are first geolocated as if they were subaerial returns (i.e., as if there were no water present), and then a refraction correction is applied; and (2) Riegl VQ-880-G bathymetric lidar data have been extensively investigated [24][25][26][27], such that new refraction correction techniques can be validated by comparing the results against those obtained by Quantum Spatial. Tested on the Lake Tahoe data set, the refraction correction algorithm described above yielded root mean square (RMS) differences from those obtained by Quantum Spatial of 0.02 m vertical and 0.14 m horizontal, even with a simplifying assumption of a flat water surface.

Study Site
The area selected for the ICESat-2 bathymetry validation encompasses the waters south of St. Thomas, U.S. Virgin Islands (USVI), including the West Gregerie Channel, Crown Bay, Krum Bay, as well as Druif Bay and Flamingo Bay on the western end of Water Island (Figure 4). The site, which intersects the St. Thomas Harbor Area of Particular Concern (APC) [28], is heavily used by commercial and cruise ships, as well as recreational boat traffic [29]. Interest in the site stems from both its economic importance to the island and threats to the health of coral reef ecosystems, which provide habitats for a number of species [28,[30][31][32].

ATLAS Data
ICESat-2 flew over the St. Thomas project site on November 22, 2018 at approximately 06:05 (UTC) on a descending track (North to South trajectory). The ATL03 geolocated photon returns for the strong beam in the westernmost ground track pair (GT3R) are shown in Figure 5. These data come from the initial data release (release 001) available at the National Snow and Ice Data Center (NSIDC: www.nsidc.org). To date, this mission is still considered within early on-orbit operations, as the calibration and validation studies contribute to continued improvement in the geolocation knowledge and pointing control capabilities. As such, there are still known biases in the range and position of the ATL03 product that can be identified and removed in those locations where reference data from other sources are available for comparison. Distinct geophysical or manmade features present in both datasets (truth and measured) create a basis for position alignment to determine geolocation (easting, northing) and range offsets within the satellite transects [11]. For the St. Thomas region, the bias identification and correction was possible using a high-resolution airborne topographic lidar data set acquired by Photo Science, Inc. (now part of Quantum Spatial, Inc.) for NOAA in November-December, 2013 with a Leica ALS70 lidar [35]. By fortunate coincidence, the track used in this study also passed directly over the Cyril E. King Airport (STT) in Charlotte Amalie West, St. Thomas, crossing Runway 10/28 and the parallel taxiway. The runway and other topographic features within the airborne collection, shown in Figure 5, were used to identify the ICESat-2 ATL03 geolocation and elevation biases prior to the bathymetric assessment.
The geolocation bias correction used the topographic (Leica) lidar data as a guide to determine the geolocation biases and vertical error in the satellite transect. Importantly, this Leica topographic lidar data set was independent of the EAARL-B bathymetric lidar data set used in the subsequent accuracy assessment of the refraction-corrected ATLAS bathymetry. The geolocation assessment methodology rasterizes the reference dataset at variable resolutions and iteratively determines the fit of the relevant ICESat-2 data by minimizing the error in the vertical direction. The initial offsets are typically determined using a coarse 8 m resolution with the final solution corresponding to the optimal fit between the ATL03 transect and the finer 1 m resolution reference grid. Although the analysis uses only the vertical error minimization, these cases with distinct geophysical or manmade features ensure that the horizontal offsets support the accuracy in the z-direction, as horizontal errors cause significant changes in the elevation if misrepresented. The realization of this process is the easting, northing and elevation offsets of the ICESat-2 data as the true geolocation and height knowledge of the satellite data product. Figure 6 shows the automated output from the analysis to visibly see the relationship between the vertical mean absolute errors and the easting and northing offset values between the 1.7 km coincident ATL03 transect to the Leica data as well as the best fit solution using biases of 2 m, 0 m, and -30 cm (easting, northing and vertical). These numbers specify that the current geolocation knowledge of the ATL03 data in this area is 2 m to the west and 31 cm above the 1 m resolution Leica surface model. The geolocation bias correction used the topographic (Leica) lidar data as a guide to determine the geolocation biases and vertical error in the satellite transect. Importantly, this Leica topographic lidar data set was independent of the EAARL-B bathymetric lidar data set used in the subsequent accuracy assessment of the refraction-corrected ATLAS bathymetry. The geolocation assessment methodology rasterizes the reference dataset at variable resolutions and iteratively determines the fit of the relevant ICESat-2 data by minimizing the error in the vertical direction. The initial offsets are typically determined using a coarse 8 m resolution with the final solution corresponding to the optimal fit between the ATL03 transect and the finer 1 m resolution reference grid. Although the analysis uses only the vertical error minimization, these cases with distinct geophysical or manmade features ensure that the horizontal offsets support the accuracy in the z-direction, as horizontal errors cause significant changes in the elevation if misrepresented. The realization of this process is the easting, northing and elevation offsets of the ICESat-2 data as the true geolocation and height knowledge of the satellite data product. Figure 6 shows the automated output from the analysis to visibly see the relationship between the vertical mean absolute errors and the easting and northing offset values between the 1.7 km coincident ATL03 transect to the Leica data as well as the best fit solution using biases of 2 m, 0 m, and -30 cm (easting, northing and vertical). These numbers specify that the current geolocation knowledge of the ATL03 data in this area is 2 m to the west and 31 cm above the 1 m resolution Leica surface model.

Empirical Depth Mapping Performance Determination
In addition to the St. Thomas site, ICESat-2 data were inspected over three additional sites to understand, globally, the bathymetric capabilities of the satellite in terms of depth penetration (Table  1). These sides include: Turks and Caicos (Figure 7), North West Australia (western Pilbara region) Figure 6. Determination of Easting, Northing, and Height correction to the St. Thomas ATLAS geolocated photon returns using the Leica topographic lidar data. It is important to note that the Experimental Advanced Airborne Research Lidar-B (EAARL-B) topobathymetric lidar data, which was used as the reference data set in assessing the accuracy of the ATLAS-derived bathymetry, was not used in the geolocation correction.

Empirical Depth Mapping Performance Determination
In addition to the St. Thomas site, ICESat-2 data were inspected over three additional sites to understand, globally, the bathymetric capabilities of the satellite in terms of depth penetration (Table 1). These sides include: Turks and Caicos (Figure 7), North West Australia (western Pilbara region) (Figure 8), and the Great Bahama Bank (Figure 9). These sites provide a diversity in location as well as feature content but, more importantly, were instances where the ICESat-2 ATL03 transects exhibited a clear extinction depth. These depths represent the point at which the downward-sloping seabed returns are no longer visible in the data, allowing for characterization of the satellite depth mapping performance.  Figure 6. Determination of Easting, Northing, and Height correction to the St. Thomas ATLAS geolocated photon returns using the Leica topographic lidar data. It is important to note that the Experimental Advanced Airborne Research Lidar-B (EAARL-B) topobathymetric lidar data, which was used as the reference data set in assessing the accuracy of the ATLAS-derived bathymetry, was not used in the geolocation correction.

Empirical Depth Mapping Performance Determination
In addition to the St. Thomas site, ICESat-2 data were inspected over three additional sites to understand, globally, the bathymetric capabilities of the satellite in terms of depth penetration (Table  1). These sides include: Turks and Caicos (Figure 7), North West Australia (western Pilbara region) (Figure 8), and the Great Bahama Bank (Figure 9). These sites provide a diversity in location as well as feature content but, more importantly, were instances where the ICESat-2 ATL03 transects exhibited a clear extinction depth. These depths represent the point at which the downward-sloping seabed returns are no longer visible in the data, allowing for characterization of the satellite depth mapping performance.

Results
The comparison of the z-errors (differences between EAARL-B reference bathymetry and ATLAS bathymetry) for the St. Thomas site allows for a quantitative performance assessment of the ICESat-2 depths with and without refraction corrections applied. Figure 10 provides an image of the ICESat-2 ground track over the St. Thomas region and the corresponding elevation profiles with and without the refraction correction. Similar to the previous comparison to the Leica (topographic) data raster surface, the EAARL-B data were gridded at both 10 m and 1 m resolutions to support the process of optimizing the offset determination and investigate the impact of refraction on the results.

Results
The comparison of the z-errors (differences between EAARL-B reference bathymetry and ATLAS bathymetry) for the St. Thomas site allows for a quantitative performance assessment of the ICESat-2 depths with and without refraction corrections applied. Figure 10 provides an image of the ICESat-2 ground track over the St. Thomas region and the corresponding elevation profiles with and without the refraction correction. Similar to the previous comparison to the Leica (topographic) data raster surface, the EAARL-B data were gridded at both 10 m and 1 m resolutions to support the process of optimizing the offset determination and investigate the impact of refraction on the results. The vertical errors associated with the bathymetric comparison both without and with the refraction correction are shown in Figures 11-12. The top plot in each figure provides the profiles of the measured data (ICESat-2) and the truth data (EAARL-B), over the full ~6 km length of the track beneath the water surface. The spatial variation of the errors the analysis is also shown for each 100 m segment along track. These errors indicate the increasing z-error with depth when the refraction correction is omitted, as would be anticipated. From comparing Figures 11 and 12, it is seen that by including the refraction correction, the mean absolute error decreased by 85%, and the RMSE reduced to 0.43 m. (As a side comment on Figures 10-12, the shore-adjacent data gap seen between ~0.5 and The vertical errors associated with the bathymetric comparison both without and with the refraction correction are shown in Figures 11 and 12. The top plot in each figure provides the profiles of the measured data (ICESat-2) and the truth data (EAARL-B), over the full~6 km length of the track beneath the water surface. The spatial variation of the errors the analysis is also shown for each 100 m segment along track. These errors indicate the increasing z-error with depth when the refraction correction is omitted, as would be anticipated. From comparing Figures 11 and 12, it is seen that by including the refraction correction, the mean absolute error decreased by 85%, and the RMSE reduced to 0.43 m. (As a side comment on Figures 10-12, the shore-adjacent data gap seen between~0.5 and 2.0 sec was only observed in the St. Thomas data set for the strong and weak beams of the right beam pair. From analyzing the data, we believe this data gap is due to a combination of two factors: 1) nearshore effects (e.g., breaking waves, bubbles, foam and/or suspended sediment close to shore); and 2) a darker (i.e., less reflective at 532 nm) substrate just south of the runway).
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 20 Figure 11. Elevation differences between EAARL-B ("truth") and ATLAS bathymetry at a 1 m grid resolution without the refraction correction in z. RMSE, root mean square error; MAE, mean absolute error. Figure 11. Elevation differences between EAARL-B ("truth") and ATLAS bathymetry at a 1 m grid resolution without the refraction correction in z. RMSE, root mean square error; MAE, mean absolute error.
We next repeated this analysis with the data for each ATLAS beam that intersected the available EAARL-B reference bathymetry, including the weak and strong beams of the center beam pair, and the weak beam of the right beam pair (in addition to the strong beam of the right beam pair, which was described above). The left beam pair was outside of the area for which EAARL-B data were available in this study. As anticipated, the density of seafloor points is significantly lower in the weak beam photon returns (by a factor of~4), which makes manual identification of the seafloor trace more difficult. However, we were able to identify seafloor points in the geolocated photon returns from all ATLAS beams that intersected the EAARL-B data, including the weak beams. The results ( Table 2) show similar accuracies between the weak and strong beam photon returns, after refraction correction. Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 20 Figure 12. Elevation differences between EAARL-B ("truth") and ATLAS bathymetry at a 1 m grid resolution with the refraction correction in z. Next, the depth-penetration capability of ATLAS was empirically assessed by analyzing the maximum depths of the refraction-corrected bottom return photons in all four sites (Turks and Caicos, North West Australia, and Great Bahama Bank, in addition to the St. Thomas site, which was also used in the comparison against the EAARL-B reference bathymetry described above). The maximum refraction-corrected depths for these four sites were considered in relation to Visible Infrared Imaging Radiometer Suite (VIIRS) Kd(490) data, obtained from the NOAA National Environmental Satellite, Data, and Information Service (NESDIS) Center for Satellite Applications and Research (STAR) in the form of netCDF (Network Common Data Form) files. The algorithm used by NOAA NESDIS to generate the Kd(490) product is described in [36] and summarized in [17]. Fortunately, VIIRS Kd(490) data were available for all sites from the same date as the ICESat-2 overpass. Although it would have been possible to convert from Kd(490) to Kd(532) following methods used by the airborne lidar bathymetry (ALB) company Teledyne Optech International [37], because  Next, the depth-penetration capability of ATLAS was empirically assessed by analyzing the maximum depths of the refraction-corrected bottom return photons in all four sites (Turks and Caicos, North West Australia, and Great Bahama Bank, in addition to the St. Thomas site, which was also used in the comparison against the EAARL-B reference bathymetry described above). The maximum refraction-corrected depths for these four sites were considered in relation to Visible Infrared Imaging Radiometer Suite (VIIRS) K d (490) data, obtained from the NOAA National Environmental Satellite, Data, and Information Service (NESDIS) Center for Satellite Applications and Research (STAR) in the form of netCDF (Network Common Data Form) files. The algorithm used by NOAA NESDIS to generate the K d (490) product is described in [36] and summarized in [17]. Fortunately, VIIRS K d (490) data were available for all sites from the same date as the ICESat-2 overpass. Although it would have been possible to convert from K d (490) to K d (532) following methods used by the airborne lidar bathymetry (ALB) company Teledyne Optech International [37], because the difference between the two is generally within a few percent, and the conversion requires additional information on water type that was not available for the sites used in this analysis, it was not performed in this study.
For consistency with the JALBTCX (Joint Airborne Lidar Bathymetry Technical Center of Expertise) coastal mapping and charting community, we chose to assess and report the bathymetric mapping performance of ATLAS following the same conventions used for commercial ALB systems. Specifically, the maximum depth measurement capability is generally specified as a function of water clarity, which, in turn, is sometimes quantified by the diffuse attenuation coefficient, K d , but perhaps more frequently by Secchi depth, Z SD (e.g., [23]). Seafloor reflectance at the laser wavelength, ρ 532 , is, of course, also a factor. While ρ 532 data were not available for the four sites, all four appear to contain mostly sandy substrates, with some coral (St. Thomas and Manihiki Island) and limestone substrates (Great Bahama Bank), which can be considered ideal for maximum depth mapping with lidar.
Secchi depths, which have been used in optical oceanography since 1865 [38], are measured by lowering a circular disk into the water and recording the depth at which it is no longer visible. (As an aside, there are a number of important considerations that must be taken into account to obtain reliable measurements for this seemingly simple parameter, including the type of disk used [39]). Notwithstanding the challenges of Secchi depth measurement, it is often preferred to quantify the maximum depth mapping performance of a lidar system in terms of Secchi depth, due to the intuitive nature of the metric.
Although no exact conversion between K d and Secchi depth exists, or is even theoretically possible [40], a number of empirical relationships, obtained from various regressions, have been developed and are in widespread use in the ALB community. One such relationship that is used in the field of ALB for coastal waters [41] is where Z SD is in meters, and K d is in inverse meters. This relationship is listed in [41] as being valid for the range: 4 m ≤ Z SD ≤ 43 m, or 0.06 m −1 ≤ K D ≤ 0.32 m −1 . A more general relationship [38,39] is In Equation (15), c is the beam attenuation coefficient, and R is the irradiance ratio, defined as the ratio of upward to downward irradiances at the surface [39]. Equation (15) is seemingly less useful than Equation (14), due to the dependence on two additional parameters, c and R. However, [39] demonstrate the relationship of both of these parameters to one additional parameter, the single scattering albedo, ω 0 , which is the dimensionless ratio of total scattering to total attenuation. By digitizing the data in Figure 7 of [39] using the online program "Web Plot Digitizer" [42] and performing a linear regression, we obtain the following relationship: Taking a typical value of ω 0 of 0.85 for coastal waters [39,41], Equation (16) further simplifies to: Equations (14) and (17) provide roughly consistent conversions between K d and Z SD , as shown in Figure 13. The mean of the two conversions, which is also depicted, was used in this work. Remote Sens. 2019, 11, x FOR PEER REVIEW 16 of 20 For very clear waters (Kd < 0.06 m -1 ), the relationships expressed in Equation 14 and 17 and shown graphically in Figure 8 cannot be assumed to hold. In this case, we use the following relationship, given in [41] and attributed to Poole-Atkins: The results of the depth penetration capability analysis for all four sites are shown in Table 3, with the means tabulated in the bottom row.

Discussion
The results show that, for the four areas tested, the maximum depth penetration for ATLAS was 0.96 Secchi depth, while the maximum optical depth [39], KdDmax, was 1.8, with reasonable consistency in values across the four sites (standard deviation of 0.10 and 0.17 for the maximum depth penetration expressed in Secchi depths and as KdDmax, respectively). Combined with the 0.43-0.60 RMSE from the comparison of refraction-corrected ATLAS bathymetry against high-accuracy airborne bathymetric lidar in the St. Thomas study site, these results provide a strong indication of the potential for generating nearshore bathymetry in a number of coastal regions around the world, with sufficient  Figure 8 cannot be assumed to hold. In this case, we use the following relationship, given in [41] and attributed to Poole-Atkins: The results of the depth penetration capability analysis for all four sites are shown in Table 3, with the means tabulated in the bottom row.

Discussion
The results show that, for the four areas tested, the maximum depth penetration for ATLAS was 0.96 Secchi depth, while the maximum optical depth [39], K d D max , was 1.8, with reasonable consistency in values across the four sites (standard deviation of 0.10 and 0.17 for the maximum depth penetration expressed in Secchi depths and as K d D max , respectively). Combined with the 0.43-0.60 RMSE from the comparison of refraction-corrected ATLAS bathymetry against high-accuracy airborne bathymetric lidar in the St. Thomas study site, these results provide a strong indication of the potential for generating nearshore bathymetry in a number of coastal regions around the world, with sufficient accuracy for a number of coastal science applications (e.g., benthic habitat mapping and the investigation of nearshore processes).
ATLAS ATL03 geolocation accuracy is expected to continue to improve, due to ongoing enhancements to the geolocation and calibration at NASA [43]. Therefore, it is recommended that ATLAS's bathymetric mapping accuracy continue to be assessed, using the latest ATL03 data products and the refraction correction developed in this work. Ideally, these follow-on analyses would include empirically assessing ATLAS's bathymetric mapping accuracy in additional geographic locations, using temporally coincident multibeam echosounder (MBES) data. While the EAARL-B and ATLAS are both lidar systems, their design characteristics and acquisition parameters are vastly different, with the EAARL-B being a dedicated airborne topobathymetric mapping system capable of meeting IHO S-44 hydrographic survey specifications for total vertical uncertainty (TVU) [33], which justifies the use of the EAARL-B data as a reference data set. However, the use of MBES data in future work would enable interesting comparisons between acoustic (sonar) data and spaceborne lidar data. In addition to the empirical accuracy assessments, it would also be valuable to develop a total propagated uncertainty (TPU) model for ATLAS bathymetry.
Another recommendation is to include additional sites in future investigation of the maximum depth mapping performance analysis to determine whether the results of this study hold over time. While our interest is primarily in coastal and estuarine regions, it would also be interesting to conduct similar analyses in inland freshwater sites. It should also be noted that there are a number of performance enhancements to the methods developed in this study, which could be investigated in follow-on research. The step of identifying seafloor returns from the ATLAS ATL03 georeferenced photons, which was performed manually in this work, could be automated, possibly leveraging recent work in detection of weak echoes in photon-counting lidar [44], and a variation of the random sample consensus (RANSAC) algorithm.

Conclusions
This study empirically investigated refraction-corrected bathymetry from NASA's ICESat-2 Advanced Topographic Laser Altimeter System (ATLAS) geolocated photon returns, through comparison against high-accuracy reference data, in the form of airborne bathymetric lidar collected by USGS and NOAA, for a study site in the U.S. Virgin Islands. Additionally, ATLAS's depth penetration capability, as a function of water clarity, was analyzed using Visible Infrared Imaging Radiometer Suite (VIIRS) K d (490) data for four project sites, from the Caribbean Sea to the Indian Ocean. The maximum depth mapping performance of nearly 1 Secchi depth, based on analysis for the four different locations, is noteworthy and, in fact, only slightly less than that of some dedicated airborne topobathymetric lidar systems (e.g., [23]). Additionally, the close visual agreement of the refraction corrected ATLAS bathymetry with USGS and NOAA Experimental Advanced Airborne Research Lidar-B (EAARL-B) bathymetry, along with the computed RMSE of 0.43-0.60 m, demonstrate that ATLAS can accurately measure bathymetry, if the change in direction and speed of the laser pulse at the air-water interface is accounted for via the refraction correction algorithm.
Future work for our project team will include investigating the fusion of ICESat-2 and Landsat 8 or Sentinel-2 for an active/passive bathymetric mapping approach [17]. A specific science goal to be addressed through this work is the long-term monitoring of coral reef ecosystems, which is currently hindered by the lack of nearshore bathymetry. Based on the anticipated value of ATLAS bathymetry for these scientific uses, we recommend that a new bathymetric data product from ICESat-2 be considered. Lastly, if the ICESat-2 bathymetry proves as useful as early indications suggest, there may be sufficient interest to consider recommending a future satellite mission dedicated specifically to global nearshore bathymetric mapping.