Analysis of MABEL Bathymetry in Keweenaw Bay and Implications for ICESat-2 ATLAS

Abstract

In 2018, the National Aeronautics and Space Administration (NASA) is scheduled to launch the Ice, Cloud, and land Elevation Satellite-2 (ICESat-2), with a new six-beam, green-wavelength, photon-counting lidar system, Advanced Topographic Laser Altimeter System (ATLAS). The primary objectives of the ICESat-2 mission are to measure ice-sheet elevations, sea-ice thickness, and global biomass. However, if bathymetry can be reliably retrieved from ATLAS data, this could assist in addressing a key data need in many coastal and inland water body areas, including areas that are poorly-mapped and/or difficult to access. Additionally, ATLAS-derived bathymetry could be used to constrain bathymetry derived from complementary data, such as passive, multispectral imagery and synthetic aperture radar (SAR). As an important first step in evaluating the ability to map bathymetry from ATLAS, this study involves a detailed assessment of bathymetry from the Multiple Altimeter Beam Experimental Lidar (MABEL), NASA’s airborne ICESat-2 simulator, flown on the Earth Resources 2 (ER-2) high-altitude aircraft. An interactive, web interface, MABEL Viewer, was developed and used to identify bottom returns in Keweenaw Bay, Lake Superior. After applying corrections for refraction and channel-specific elevation biases, MABEL bathymetry was compared against National Oceanic and Atmospheric Administration (NOAA) data acquired two years earlier. The results indicate that MABEL reliably detected bathymetry in depths of up to 8 m, with a root mean square (RMS) difference of 0.7 m, with respect to the reference data. Additionally, a version of the lidar equation was developed for predicting bottom-return signal levels in MABEL and tested using the Keweenaw Bay data. Future work will entail extending these results to ATLAS, as the technical specifications of the sensor become available.

1. Introduction

NASA’s upcoming Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) mission, with a 91-day repeat period and near-polar orbit [1], provides a unique opportunity to assess ice sheet elevation change. The sole instrument aboard ICESat-2 will be the Advanced Topographic Laser Altimeter System (ATLAS), a micro-pulse, photon-counting lidar system operating at 532 nm using a frequency-doubled neodymium-doped yttrium aluminum garnet (Nd:YAG) laser [1,2]. Although the primary objectives of ICESat-2 will be to measure ice-sheet elevations, sea-ice thickness, and global biomass [1], the ATLAS design specifications suggest the feasibility of bathymetry retrieval from ATLAS data. If it can be demonstrated that reliable bathymetry estimation from ATLAS is, in fact, possible, this could greatly benefit studies of coastal and inland water bodies, which are often hindered by a dearth of shallow-water bathymetry [3]. The fact that ATLAS is an active remote sensing system makes it of particular interest for bathymetric mapping, since it may be possible to synergistically fuse water depth estimates from ATLAS with bathymetry retrieved from a variety of sources, including passive multispectral satellite imagery [4,5,6,7], such as from Landsat 8 and Sentinel-2, and synthetic aperture radar (SAR) data [8,9,10,11]. Using ATLAS to measure bathymetry and complement other bathymetric techniques is well aligned with NASA’s Applied Sciences Program mission to discover innovative and practical uses for Earth observations [12]. Additionally, the bathymetric mapping capability would represent an important advancement over ICESat-2’s predecessor, the original ICESat, which carried the Geoscience Laser Altimeter System (GLAS). While GLAS did contain a green laser, it was designed with a different laser architecture [13] for the secondary mission of measuring optically thin clouds and atmospheric aerosols [14,15]. GLAS used a 1064-nm laser, incapable of penetrating the water column to provide bathymetry, for elevation mapping.

In empirically evaluating the potential to derive bathymetry from ATLAS, a logical place to start is with NASA’s airborne photon-counting ATLAS simulator, the Multiple Altimeter Beam Experimental Lidar (MABEL), which has been flown on the NASA Earth Resources 2 (ER-2) high-altitude aircraft on over 50 missions since 2010. It is important to note that MABEL is not an exact model of what ATLAS will be [2], but its main purposes are to verify the ATLAS measurement concept and provide data similar to ATLAS data for algorithm development.

In keeping with the above considerations, our long-range research goals are to: (1) assess the ability to derive bathymetry from MABEL; (2) evaluate the accuracy and reliability of MABEL bathymetry; (3) extend the analysis of MABEL data to ATLAS; and (4) (pending the results of the previous phases) develop tools and techniques for bathymetric mapping using fusion-based approaches that leverage ATLAS bathymetry. The motivation for our long-range research program stems from the large number of nearshore, coastal areas around the world for which bathymetric data are entirely lacking, and the need for such data for applications ranging from modeling inundation, due to storm surge and sea level rise, to nautical charting [3]. The broad objective of obtaining coastal bathymetry worldwide is entirely aligned with initiatives to manage, expand, and add value to the suite of global bathymetry obtained through a variety of techniques by a number of national and international organizations [16,17,18,19,20]. Due to the challenges of mapping bathymetry in very shallow waters with any one particular technology [3], the synergistic fusion of spaceborne, active sensor bathymetric data with complementary data sets, including airborne bathymetric lidar, shipborne sonar, airborne and spaceborne multi-and hyperspectral imagery, and synthetic aperture radar (SAR), affords the best option for filling the voids in worldwide shallow bathymetric coverage.

This paper presents the results of the first two phases of the long-range research program outlined above. MABEL Viewer, a web-application developed in this study to facilitate interactive identification of MABEL bottom returns, was used to identify bathymetric returns in MABEL data acquired over Keweenaw Bay in Lake Superior. To quantitatively assess the quality of bathymetry from MABEL, the identified bottom returns were first corrected using National Oceanic and Atmospheric Administration (NOAA), Great Lakes Coastal Forecasting System (GLCFS) data and a first-order refraction correction. These corrected depths were compared against high-accuracy bathymetry acquired for NOAA with a Fugro Laser Airborne Depth Sounder (LADS) Mk II in 2010. The results show good agreement between MABEL bathymetry and the NOAA reference bathymetry, with a z-component root mean square error (RMSE_z) of 0.7 m in water depths of up to 8 m. Concurrently with this analysis, a version of the lidar equation was developed for predicting the number of bottom returns in the data set as a function of depth within the project site. The predictions were compared against the observed data to assess the general validity of our baseline model. We conclude with a discussion of the next steps, including the extension of this work to ATLAS, using simulations generated by project partners at the University of Houston, National Center for Airborne Laser Mapping (NCALM).

1.1. Bathymetric Lidar Background

Almost as old as the laser itself, airborne lidar (light detection and ranging) has its roots in the mid-1960s. The reader is directed to Guenther [21] and Petrie and Toth [22] for a detailed history of airborne lidar remote sensing, including the first practical demonstration of the capability of measuring depth in 1969 [23]. Lidar technology evolved rapidly over subsequent decades, leveraging developments in enabling technologies, including inertial and satellite navigation systems [21,22]. Conventionally, lasers operating in the near infrared (typically, 1064 and 1550 nm) have been preferred for topographic lidar, while bathymetric lidar has used water-penetrating 532 nm lasers, sometimes in combination with 1064 nm for obtaining water surface returns and simultaneously mapping topography. Industry-standard bathymetric lidar systems include the Optech Scanning Hydrographic Operational Airborne Lidar Survey (SHOALS) system and Coastal Zone Mapping and Imaging Lidar (CZMIL), Leica Airborne Hydrography AB (AHAB) Chiroptera II and HawkEye II and III, National Aeronautics and Space Administration (NASA) Experimental Advanced Airborne Research Lidar (EAARL), United States Geological Survey (USGS) EAARL-B, and Fugro LADS Mk II and Mk 3 [21,24,25]. Even more varied than the number of systems are the applications of these systems, including benthic habitat classification [26,27,28,29], shoreline delineation [30], inundation modeling [31], river morphodynamics [32,33], geomorphological mapping [34,35,36], and hydrography [21,37,38]. The conventional bathymetric lidar systems listed above are all waveform resolving, meaning they record and/or analyze in real time digitized, backscattered pulses containing up to thousands of photons.

In the topographic lidar community, an emerging technology is photon-counting lidar, which uses either photomuliplier tubes (PMTs) or avalanche photodiodes (APDs) operated well above the breakdown voltage (i.e., in Geiger mode) to record returns as low as a single photon as individual points. The low power consumption of photon-counting lidar systems, which also have low signal-to-noise ratios (SNR), make them ideal systems for spaceborne platforms with strict power budgets. However, photon-counting lidar systems, despite first being identified as a potential mapping technology in 1996 [39], have received considerably less attention in the ocean and coastal mapping community. Only with advances in photon-detector ranging resolution and dead time have photon-counting systems become viable mapping tools [40]. A number of recent publications target terrestrial and cryospheric applications, such as canopy extraction [41,42,43] and sea-ice and glacier profiling [44,45,46,47], and noise-filtering techniques [48,49,50], but research pertaining to using photon-counting lidar for bathymetric mapping remains limited.

Examples of photon-counting systems that have been shown to detect bathymetry are the Coastal Area Tactical-mapping System (CATS) [51] and the NASA Airborne Multi-kilohertz Microlaser Altimeter Instrument [40]. MABEL has also been shown to detect bathymetry, in the low-turbidity waters of Lake Mead and depths of up to approximately one Secchi depth [52]. Focused primarily on the retrieval of water surface height statistics, Jasinski et al. [52] includes a brief qualitative discussion of depth profiles observed in the photon-elevation data. However, the study did not account for index-of-refraction and vertical-datum corrections or include a quantitative comparison with existing bathymetry from dedicated bathymetric-mapping instruments. Whereas CATS and the NASA microaltimeter were operational at heights above ground level (AGLs) of 500 to 6700 m, MABEL is operational at an AGL of up to 20,000 m. A detailed investigation of MABEL, which is flown above 95% of Earth’s atmosphere, is, therefore, an important step toward evaluating the possibility of measuring bathymetry from photon-counting, spaceborne lidar.

1.2. Lidar Equation

In addition to the empirical analysis, it is also of interest to develop a model for predicting the bottom return signal level, given a specified set of environmental conditions and system parameters. Such a model is not only beneficial for evaluating the expected depth ranges in which bathymetry retrieval might be viable using MABEL, but, in future work, it could be extended to ATLAS, as the final system parameters are made available by NASA. In this study, the expected number of photoelectrons corresponding to bottom (seafloor) returns, $n_{bot}$, is evaluated using a modified version of the lidar equation in [53]:

$$n_{bot}={\eta }_q{\eta }_t{\eta }_r\frac{E_t}{hv}{\rho }_{\lambda }\cos \left({\alpha }_t\right)\frac{A_r}{\pi {\left(R_{air}+R_w\right)}^2}{\left(1-r_{int}\left({\alpha }_s\right)\right)}^2{T}_{\lambda }^2{e}^{-2c_{\lambda }{R}_w},$$

(1)

In Equation (1), h is Planck’s constant; v is the photon frequency; R_w is the pulse travel distance in water; $R_{air}$ is the pulse travel distance in air; ${\eta }_q$ is the detector quantum efficiency; ${\eta }_t$ is the transmitter optical efficiency; ${\eta }_r$ is the receiver optical efficiency; $E_t$ is the transmitted energy per channel per pulse; ${\rho }_{\lambda }$ is the bottom reflectance at the laser wavelength, λ; ${\alpha }_t$ is the incidence angle on the bottom; $A_r$ is the collecting area of receiver aperture; $r_{int}$ is the air-water interface reflectance; $T_{\lambda }^2$ is the two-way atmospheric transmittance; ${\alpha }_s$ is the local incidence angle on water surface; and $c_{\lambda }$ is the effective attenuation coefficient (for which the water column diffuse attenuation coefficient, K_d, establishes a lower bound). The two-way atmospheric transmittance can be further expressed as:

$$T_{\lambda }^2=e^{-2{\int }_0^H\sigma \left(z\right)dz},$$

(2)

where $\sigma \left(z\right)$ is the total atmospheric extinction coefficient, at the laser wavelength, as a function of altitude, z, and H is the altitude of the ER-2 aircraft at the time of the MABEL overflight.

2. Materials and Methods

2.1. MABEL Data

In missions spanning from December 2010 to August 2014, MABEL data have been acquired across North America, over a variety of landforms, marine environments, and vegetative regimes, including inland and coastal mountain ranges, Arctic sea ice and ice sheets, prairies, deserts, savannah, and coastal plains. This study focuses on two sites along the “Transit to KPMD” mission (where KPMD refers to the Palmdale Regional Airport, in Palmdale, CA, USA) identified using MABEL Viewer, an in-house prototype data-explorer web application built on Python and JavaScript. This custom software tool facilitates examining the existing MABEL data archive by providing a graphical user interface (Figure 1) to select a flight line of interest and interactively view and classify the photon elevation data for each channel. During future phases of the project, MABEL Viewer is intended to also serve as an algorithm-development environment and a template for distributing MABEL data.

The two study sites are in Lake Superior, at the eastern base of Michigan’s Keweenaw Peninsula (Figure 2). No local water-clarity information is available in the study area during the time of data acquisition, but in general, Lake Superior is an optically complex, oligotrophic case 2 water body [54]. Dominant coastal substrates in the region include high-reflectance white quartz sand derived from Jacobsville sandstone, lower-reflectance gray stamp sands, a product of the region’s historic copper mining, and localized cobble fields [55,56,57].

MABEL has 16 green (532 nm) channels and eight near infrared (1064 nm) channels, with a laser pulse length of 2 ns and an operational laser repetition rate of 5–20 kHz. The geometry of the channels is defined by selecting various fibers within the 215-fiber transmitter and receiver arrays located at the foci of the corresponding matched transmitter and receiver telescopes. This study examines data acquired during the “Transit to KPMD” mission, for which only 13 of the 16 green channels were configured. The corresponding beam angles, elevations, and relative energy levels for this mission are listed in

Table 1. MABEL channel specifications for the 13 green channels configured for the “Transmit to KPMD” mission.

Green Channel	Energy Level (Metadata)	Energy Level (Reclassified)	Angle (Mrad)	Elevation (Mrad)
1	Low	Low	5.0	−1.5
3	High	Low	−2.2	−0.5
4	High	Low	−0.2	−0.5
5	High	High	−5.0	−1.5
6	High	High	0.0	−1.5
7	Low	Low	1.8	−0.5
8	Low	High	−2.0	−1.5
9	Low	Low	2.2	−2.5
10	High	Low	0.2	−2.5
11	Low	High	2.0	−1.5
12	High	Low	−1.9	−2.0
14	Low	Low	−2.1	−1.0
15	Low	Low	−1.8	−2.5

. Only the green channels are considered in this study because ATLAS will operate only at the green wavelength and because it is the only one of the two wavelengths capable of penetrating the water column to provide bottom returns. At a nominal mission AGL of 20,000 m, the 13 green channels cover a swath of 200 m (Figure 3) and have a nominal footprint of 2 m. The channels corresponding to the anticipated low- and high-energy levels of ICESat-2, as per MABEL metadata, are also listed in Table 1; however, according to NASA, the documented power levels have not been rigorously verified and are affected by a number of environmental factors and configuration procedures [58]. To more closely represent the relative energies “as flown”, we reclassified the relative ICESat-2 energy level of each channel according to the along-track density of water-surface returns, which was used as a proxy indicator of signal strength. The energy level for each channel was then linearly scaled based on an assumed channel 6 energy level of 0.2 µJ, a value taken from a previous study examining MABEL data acquired 5 months prior under a similar energy regime [45,59]. The average low and high energy levels were calculated to be 0.04 and 0.2 µJ, up to two orders of magnitude lower than the 5–7 µJ reported in the original design specifications [2].

The MABEL data were directly georeferenced with a NovAtel GPS-aided inertial navigation system (INS) [60]. The MABEL trajectory and photon-elevation data, available from the online L2A product data archive, are relative to the WGS84 (G1674) ellipsoid, having been processed in NovAtel Inertial Explorer with a PPP (precise point positioning) post-processing paradigm, using precise ephemeris data obtained from the International GNSS (Global Navigation Satellite System) Service (IGS) [60]. The roll, pitch, and heading accuracies are documented to be, respectively, 0.007°, 0.007°, and 0.010°. Rotations between MABEL and the inertial measurement unit (IMU) reference frames and between the IMU and the aircraft reference frames are recorded in the MABEL data files.

Figure 4 shows example MABEL data with a distinctive bathymetric profile, along with discernable ground and vegetation returns. Raw photon elevations for each channel from the entire data file are shown in Figure 5. Photon classifications of water-surface and bottom were assigned manually using a classification tool in MABEL Viewer. Bottom returns were visually identified via expert knowledge of typical bathymetric profiles, in which the topographic surface can be seen to extend continuously below the water surface to the extinction depth of the lidar. The density of the data relative to ambient noise levels was also considered. The classification scheme shown in Figure 4 and Figure 5 is as follows: red = bottom; dark blue = water surface; light blue = all others (vegetation, land and noise classes).

2.2. Reference Data

The reference data set used in this study was a 2010 NOAA Coastal Services Center data set acquired with the Fugro LADS Mk II airborne bathymetric lidar system [61]. A high-power, full-waveform lidar, the LADS Mk II was operated at an altitude of 365 to 670 m, at 140 to 175 kts, with a pulse rate of 900 Hz and a scan rate of 18 Hz. Horizontal control, referenced to NAD83(CORS96), a realization of the North American Datum of 1983 based on Continuously Operating Reference Stations (CORS), was achieved with a control network consisting of two dedicated GPS base stations and seven CORS. Vertical control, referenced to the International Great Lakes Datum 1985 (IGLD85), was based on water-level data from four NOAA National Water Level Observation Network (NWLON) stations and a co-tidal model. The data set is documented in the distributed metadata as having a horizontal accuracy of 2.97 m and a vertical accuracy, after transformation to the Geodetic Reference System 1980 (GRS80) ellipsoid, of 0.29 m. The average point density was calculated to be 0.27 points per m². The data were obtained from NOAA Digital Coast referenced to WGS84 (G1674).

2.3. Predicted Number of Photoelectrons

The predicted number of photoelectrons per pulse as a function of depth for each channel was calculated by applying Equation (1). For the purposes of this study, the expected number of photoelectrons given by Equation (1) is assumed to also be the expected number of signal events, given a signal threshold of one photoelectron. The values used in Equation (1) are summarized in

Table 2. Parameter values used in a photon-counting version of the lidar equation to predict the number of photon returns from the bottom. Parameter values are purposefully stated with variable precisions, which reflect our relative level of knowledge of each.

Variable	Value	Unit	Description
${\eta }_q$	$0.18$	-	Detector quantum efficiency (Hamamatsu PMT H7260)
${\eta }_r$	$0.4$	-	Receiver optical efficiency
${\eta }_t$	$0.4$	-	Transmitter optical efficiency
$E_t$	0.04 and 0.2	µJ	Transmitted energy per channel pulse
${\rho }_{\lambda }$	$0.3$	-	Bottom reflectance at laser wavelength, λ
${\alpha }_t$	0.1	rad	Incidence angle on lake bottom
$A_r$	$0.01$	m2	Collecting area of receiver aperture
$R_{air}$	$20,000$	m	Pulse travel distance in air
$r_{int}$	$0.1$	-	Reflectance of air-water interface
$T_{\lambda }$	$0.8$	-	One-way atmospheric transmittance
$c_{\lambda }$	$0.25$	m−1	Effective total beam attenuation coefficient

. In order of preference, the values were obtained from: (1) NASA specifications, data sheets and correspondence (nominal energy level, collecting area of receiver aperture, pulse travel distance in air); (2) manufacturer specifications (detector quantum efficiency); (3) typical values for similar systems and environmental conditions obtained from published papers or reports (bottom reflectance); or (4) reasonable, assumed values (optical efficiencies, air-water reflectance, atmospheric transmittance, and effective total beam attenuation coefficient). The values stated in the table are purposefully listed with variable precisions, which reflect our relative level of knowledge of each. Project partners at the University of Houston are concurrently working on more rigorous simulations [62], which will be continually updated as new/improved system specifications for MABEL and ATLAS become available.

The expected numbers of photoelectrons were calculated in 2-D bins extending nominally 20 m along-track and 1-m vertically. Each bin corresponded to 0.1 s, or 500 laser shots, assuming an operational pulse rate of 5 kHz and speed of 200 m/s. Individual channel track lines were modeled by fitting a curve to the positions of the water-surface returns within each channel. The track lines were projected onto the bathymetric surface defined by the LADS Mk II reference bathymetry, and the average bottom elevation along each bin was used in Equation (1) to calculate the predicted number of photoelectrons for that bin. The average number of expected photoelectrons for each depth range for each channel was then calculated by averaging the populated bins at that depth range. High- and low-energy averages for each depth range were then calculated by averaging the corresponding channels.

2.4. Deriving MABEL Bathymetry

The next steps in our analysis focused on the quantitative comparison of MABEL bathymetry against the reference bathymetry. Before the identified MABEL bottom returns were analyzed for bathymetric accuracy and internal consistency, the data were corrected for water-surface refraction and vertically reduced to Lake Superior low water datum (LWD). The full workflow for generating bathymetry from MABEL data is illustrated graphically in Figure 6, and the steps which have not yet been described (i.e., those following the point classification in MABEL Viewer) are outlined below:

2.4.1. Refraction Correction

Because the native MABEL photon-elevation calculation assumes subaerial travel times only, a refraction correction was needed to account for the speed of light in the water column. In this study, a first-order depth correction factor was used, computed as the ratio of the index of refraction of air to that of fresh water, 0.7521. Planimetric corrections (which would require a water surface model for rigorous implementation) were not considered, because the maximum error associated with not applying them (~0.3 m for a 3° angle of incidence in 8 m water depth, assuming a flat water surface and no surface waves) was at the noise level for the purposes of this study. The raw depths were calculated by subtracting each bottom photon elevation from the mean elevation of all the water-surface photons in the corresponding data file.

2.4.2. Vertical Datum Transformation

The raw photon elevations were referenced to the WGS84 ellipsoid. Additional vertical corrections were necessary to account for the observation that the water level as inferred from the surface returns did not match the water level as portrayed by NOAA Great Lakes Environmental Research Laboratory (GLERL) data. The data reduction methodology used in this study incorporates both the traditional hydrographic practice of making water-level corrections based on local water-level observations [63] and the modern initiative to make “GPS-tide” corrections based on chart-datum—ellipsoid separation models [64].

Figure 7 shows the vertical differences between the water-surface as calculated from the surface-return elevations for each channel and the known water level elevation at the time of data acquisition. The average bias is 2.1 m, with a standard deviation of 0.2 m. The source of these channel-specific biases is unknown, but possible explanations include uncalibrated, channel-specific range error and any error in the z-component of the instrument installation lever-arms.

The channel-specific biases were calculated by differencing each channel’s average surface photon elevation and the elevation of the actual water level as modeled based on a regional Great Lakes Coastal Forecasting System (GLCFS) water-level value (Figure 8). The GLCFS water level, originally referenced to Lake Superior low water datum (LWD), was converted to a WGS84 (G1674) ellipsoid height based on a constant WGS84-IGLD85 separation of 35.38 m, which was derived from a model that was created using VDatum, a vertical datum transformation tool that uses a collection of transformation grids to transform vertically referenced data to and from a number of ellipsoidal, orthometric, and tidal datums [65,66]. Lake Superior LWD is defined as 183.2 m above IGLD85 [67]. The WGS84-IGLD85 separation model, shown in Figure 8, was generated by summing two component separation values at each node of a 100-m resolution 2-D grid: (1) the separation between WGS84 and NAD83 and (2) the separation between NAD83 and IGLD85. The single value used in the GLCFS water-level conversion was calculated by averaging the gridded separation values along the MABEL track line in the project site, the variation of which (σ = 3 cm) was deemed to be insignificant for the purposes of this study.

The vertical reduction methodology is summarized in Figure 9, and the corresponding variables are listed in

Table 3. Values used to reduce raw bottom elevations to WGS84 (G1674) elevations corrected for index-of-refraction and channel-specific water level biases.

Variable	Description
${\mathit{c}}_{\mathit{N}\mathit{C}}$	GLFCS water level (w.r.t. Lake Superior low water datum)
${\mathit{c}}_{\mathit{s}}$	Channel-specific water-level bias
${\Delta }_{\mathit{L}\mathit{W}\mathit{D}-\mathit{I}\mathit{G}\mathit{L}\mathit{D}}$	Lake Superior LWD–IGLD85 offset (183.2 m, by definition)
${\mathit{S}}_{\mathit{r}\mathit{a}\mathit{w}}$	Uncorrected channel water-surface ellipsoid height
${\mathit{D}}_{\mathit{r}\mathit{a}\mathit{w}}$	Raw depth (without index-of-refraction correction)
$\mathit{D}$	Depth (with index-of-refraction correction)
${\mathit{D}}^{\prime }$	D adjusted vertically for ${\mathit{c}}_{\mathit{s}}$
${\mathit{H}}_{\mathit{r}\mathit{a}\mathit{w}}$	Raw photon ellipsoid elevation
$\mathit{H}$	Final photon ellipsoid elevation
${\mathit{s}}_{\mathbf{1}}$	WGS84-NAD83 separation
${\mathit{s}}_{\mathbf{2}}$	NAD84-IGLD85 separation
${\mathit{n}}_{\mathit{a}\mathit{i}\mathit{r}}$	Index of refraction for air
${\mathit{n}}_{\mathit{w}\mathit{a}\mathit{t}\mathit{e}\mathit{r}}$	Index of refraction for fresh water

. The final ellipsoid elevation for each bottom return was calculated using Equations (3)–(6).

$$c_s={\Delta }_{LWD-IGLD}+c_{NC}-\left(s_1+s_2+s_{raw} \right),$$

(3)

$$D_{raw}=S_{raw}-H_{raw},$$

(4)

$$D^{\prime }=D_{raw}\left(\frac{n_{air}}{n_{water}}\right),$$

(5)

$$H=S+c_s-D^{\prime }$$

(6)

3. Results

3.1. Predicted Photon Returns

Figure 10 shows the predicted numbers of photoelectrons vs. depth for depths of 0–10 m, for each energy level, based on applying the values in Table 2 to Equation (1).

Figure 11 shows example results, from channel 6, comparing the expected and observed numbers of photoelectrons. The left panel displays the reference depth profile, with the along-track 2-D depth bins, and the derived MABEL depths. The right panel displays the average number of expected and observed photoelectrons per shot bin for each depth range.

Figure 12 summarizes the average corresponding numbers of expected and observed photoelectrons per pulse for each energy level.

Table 4. Average expected and observed numbers of photoelectrons per pulse for the reclassified low- and high-energy channels, which reflect the “as flown” energy regime.

Depth (m)	Low Energy Channels (0.04 µJ)			High Energy Channels (0.2 µJ)
	(1, 3, 4, 7, 9, 10, 12, 14, 15)			(5, 6, 8, 11)
	Expected	Observed	Difference	Expected	Observed	Difference
0–1	0.003	0.001	−0.003	0.017	0.005	−0.013
1–2	0.002	0.002	0.000	0.012	0.015	0.003
2–3	0.001	0.001	−0.001	0.007	0.003	−0.003
3–4	0.001	0.001	0.000	0.004	0.004	0.000
4–5	0.001	0.000	0.000	0.003	0.001	−0.001
5–6	0.000	0.000	0.000	0.002	0.002	0.000
6–7	0.000	0.000	0.000	0.001	0.001	0.000
7–8	0.000	0.000	0.000	0.001	0.000	0.000

lists the corresponding values. Overall, the number of observed photoelectrons decreases with increasing depth at a decay rate similar to the expected numbers of p.e. for each energy level. One noticeable exception is the relatively low number of observed photoelectrons in the 0–1 m depth range for both the high and low energy channels. This is not unexpected, as the shallow nearshore zone generally poses challenges for all bathymetric lidar systems [68,69,70]. The challenges associated with this depth range stem from a number of environmental factors in the surf zone, such as breaking waves, the associated bubbles, foam, and re-suspended sediment, which are often unavoidable. (As an aside, an additional challenge with waveform-resolving lidar systems in this very shallow depth range is that the convolution of the laser pulse with the water surface may result in a broadened pulse in which the bottom return is obscured.) Spatial variability of bottom reflectance, which was assumed to be constant in our model, is also a plausible contributing factor, because, as mentioned in Section 2.1, the region is characterized by high-reflectance natural white sand and lower-reflectance anthropogenic stamp sands that are transported along the coast.

The general (order-of-magnitude) agreement between the expected and observed signal levels both corroborates the observed results and provides an indication that our model produces reasonable estimates, which is important if the model is to be refined and extended to ATLAS in future work. Care should be taken to not over interpret the results of our prediction model, because the transmit energy levels have relatively high uncertainty that is difficult to quantify, given the limitations discussed in Section 2.1. The model suffices for the goals of the current phase of research: namely, it provides a means of assessing the general feasibility of deriving bathymetry from MABEL and developing a baseline model to be refined and extended in subsequent project phases.

3.2. Bathymetry Accuracy Assessment

The accuracy of the MABEL bathymetry was assessed by quantitatively comparing elevation differences between the reduced bottom-return ellipsoid heights with the reference bathymetry acquired in 2010 with the Fugro LADS MkII system. A detailed approach resolving uncertainties into vertical and horizontal components was not pursued because well-defined discrete features to use for planimetric positional comparisons were not discernable in the sparse photon data; however, we suspect that the MABEL photon positioning contains a horizontal or angular bias because the positioning of the shoreline relative to the channel track lines, which extend landward of the coast, suggests a nominal southwestward offset of the MABEL photon returns. The comparison was performed with the reduced ellipsoid heights and not the calculated depths to avoid introducing errors associated with additional vertical datum transformations. Differences were calculated for all photon returns classified as bottom, including the returns positioned up to 35 m outside of the extents of the reference data due to the suspected horizontal or angular bias.

The distribution of differences between the bottom-return ellipsoid heights and the ellipsoid heights of the nearest reference bathymetry data points, as found by a k-d tree nearest neighbor search, are shown in Figure 13. With an overall bias of −0.59 m, the two datasets agree to within an RMS difference of 0.74 m. Numerous factors could explain the general trend of the MABEL-derived depths being generally deeper than the reference surface, including depth and vertical-datum biases in the reference bathymetry and the suspected horizontal bias in the MABEL photon positioning. Conversely, the observed bias may also include a systematic underestimation of depth resulting from the possible penetration of “surface returns” into the upper portions of the water column [71]. A small, slope-induced high depth bias is also possible due to the LADS MKII system having a slightly larger footprint (nominally 2.5 m) than MABEL (2 m). A change in actual bathymetry over the intervening two years is also a likely factor, as the project area is characterized by dynamic bedform features, including nearshore sand bars and troughs and offshore shoals likely influenced by migrating sand [72]. The sedimentary regime of Keweenaw Bay has not been fully described, but a key regional sedimentary process occurring along the northern reaches of the Bay is the southwestward long-shore movement and deposition of stamp sands, which eventually migrate laterally to coalesce with deeper, larger bars [55,56]. Keweenaw Bay’s exposure to an approximately 230-km northeasterly fetch could also play a significant role in the sediment dynamics of the project area [57].

Figure 14 shows the spatial variability of the ellipsoid-height differences plotted over 1-m resolution NAIP imagery. Nearshore clusters of positive and negative differences are consistent with the combined effects of migrating sand bars and the suspected horizontal bias. Although the imagery does not temporally coincide with the MABEL and reference data, it offers insight into the observed spatial distribution of the ellipsoid height differences by revealing the dynamic geomorphological nature of the project area.

4. Discussion

As noted earlier, MABEL is not intended to be an exact replica of ATLAS.

Table 5. Comparison between MABEL and future Advanced Topographic Laser Altimeter System (ATLAS) system parameters. Note: the values in parentheses in the MABEL pulse energy field denote the reclassified, “as flown” energy levels, as opposed to the original design specifications.

Parameter	MABEL	ATLAS
Laser footprint ($1/e^2$)	2 m (100 µrad)	15 m (31 µrad)
Field of View	4.2 m (210 µrad)	41 m (83 µrad)
Laser pulse repetition freq.	5–20 kHz	10 kHz
Pulse energy	5–7 µJ (0.04–0.2 J)	41/160 µJ
Pulse pattern	16 532-nm beams, 8 1054-nm beams	6 beams (3 pairs of 2)
Swath width	2 km (max) (variable)	6 km
Wavelength	532 and 1064 nm	532 nm
Filter width	~150/~400 pm (532/1064 nm)	30 pm
Receiver aperture area	0.013 m2	0.79 m2
Operational altitude	20 km	500 km

lists some of the significant parameter differences between the two sensors.

Many of the parameter values needed to quantitatively assess ATLAS expected signal levels from the bottom in a particular area using Equation (1) are not currently available. However, given the significantly higher pulse energies and receiver aperture, it is reasonable to expect that the depth penetration capabilities of ATLAS will at least equal, if not exceed, those of MABEL. The much greater altitude of the ICESat-2 satellite (as compared with the ER-2 aircraft during the MABEL flights) impacts the reflected beam solid angle, but has relatively little impact on atmospheric transmittance, since the ER-2 operational altitude is already above ~95% of the Earth’s atmosphere. The fundamental spatial resolution of ATLAS will also be coarser, given its larger footprint.

If ICESat-2 ATLAS is found to have similar or greater bathymetric mapping capabilities than MABEL, the resulting bathymetric data may be useful for small-scale nautical charting in poorly mapped areas, as well as for inundation and storm surge modeling and related coastal zone management and science needs. Only one, small area was examined in this study, but a preliminary review of the suite of existing MABEL data reveals high-probability returns in other parts of Lake Superior, as well as a location along the Northern Icelandic coast. Additionally, Jasinski et al. [52] noted MABEL bathymetry in Lake Meade. (These data were qualitatively assessed, but not included in the quantitative analysis of bathymetry in this study, due to a lack of available reference bathymetry for the sites.)

Although the observed MABEL and anticipated ATLAS data densities and vertical uncertainties do not approach those of dedicated airborne bathymetric lidar using full-waveform systems, a satellite platform allows regular data acquisition in remote, poorly mapped areas, where high-resolution bathymetric lidar and shipborne acoustic data are typically scarce or altogether absent. Hence, there are multiple opportunities to acquire data in areas in which there are large seasonal—or even daily—variations in water clarity. Additionally, the repetitive coverage offered by repeat orbits may densify bottom returns and improve SNR.

Another potential application of ICESat-2 bathymetry data is to serve as a constraint in bathymetry retrieval from a variety of synergistic methodologies using various data sources, including passive, multispectral satellite imagery [6,7,73] and synthetic aperture radar (SAR) [9,11,74]. Imagery-derived bathymetry is inferred from the relative reflectance of different bands (e.g., blue and green), and SAR-derived bathymetry is inferred based on variations of surface roughness. Each of these methods has the advantage of providing dense depth estimates over large spatial extents; however, each method usually relies on “seed depths” typically soundings from nautical charts or lidar surveys [6] or shipboard sonar data [75], to calibrate derived values to meaningful local depths. By providing a more-direct form of depth measurement, MABEL (and, ultimately, ATLAS) bathymetry does not require existing reference soundings. On the other hand, bathymetry from MABEL or ATLAS has the disadvantage of being spatially sparse, as depths are constrained to lie along discrete track lines. MABEL (or ATLAS) bathymetry and bathymetry from multispectral satellite imagery and SAR data are, therefore, highly complementary, such that the fusion of each could potentially combine the strengths and overcome the limitations of each method individually. By providing a nearly seamless transition across the shoreline, ICESat-2 bathymetry is also a natural complement to coastal-mapping applications incorporating infrared, terrestrial lidar [76,77,78,79,80,81,82].

5. Conclusions

MABEL, the airborne simulator for NASA’s upcoming ICESat-2 mission, has been shown to detect bathymetry up to 8 m deep in the oligotrophic waters of Keweenaw Bay, Lake Superior. The depths were derived by applying index-of-refraction and water-level corrections to the raw MABEL photon elevations and were shown to be in good agreement with the NOAA reference data acquired two years earlier. Resolving positional uncertainty into vertical and horizontal components requires further research evaluating system and procedural parameters such as boresight calibration angles and georeferencing settings; and incorporating a thorough index-of-refraction model, accounting for the orientation of the aircraft, the corresponding angle of incidence of each beam relative to the water surface, and a water-surface model. Before MABEL data—and potentially ICESat-2 ATLAS data—can be reliably reduced to specific vertical data with sub-meter accuracy without directly relying on local water-level measurements at the time of data acquisition, vertical uncertainties resulting from 3-D georeferencing and channel-specific range calculations need to be refined.

In addition to assessing the depth accuracy of the derived MABEL bathymetry through a comparison with a high-accuracy reference dataset, a photon-counting version of the lidar equation was used to predict the number of signal events (taken to be one photoelectron) in 2-D bins along each channel track line. Anomalies in the 0–1 m depth range are consistent with known environmental limitations; and uncertainties in the assumed transmit-energy levels limit the applicability of our results, but the general agreement between the average observed and expected numbers of photoelectrons in the 1–8 m depth range for the low- and high-energy channels substantiates our approach and justifies developing our baseline model in future project phases.

Future research plans include developing an algorithm to automatically classify and process bottom returns from photon-counting lidar data and a corresponding uncertainty model. Research plans also include collaboration with the National Center for Airborne Laser Mapping (University of Houston) to extend this MABEL analysis to incorporate theoretical performance characteristics of the ATLAS sensor.