An Approach to Minimize Atmospheric Correction Error and Improve Physics-Based Satellite-Derived Bathymetry in a Coastal Environment

Physics-based radiative transfer model (RTM) inversion methods have been developed and implemented for satellite-derived bathymetry (SDB); however, precise atmospheric correction (AC) is required for robust bathymetry retrieval. In a previous study, we revealed that biases from AC may be related to imaging and environmental factors that are not considered sufficiently in all AC algorithms. Thus, the main aim of this study is to demonstrate how AC biases related to environmental factors can be minimized to improve SDB results. To achieve this, we first tested a physics-based inversion method to estimate bathymetry for a nearshore area in the Florida Keys, USA. Using a freely available water-based AC algorithm (ACOLITE), we used Landsat 8 (L8) images to derive per-pixel remote sensing reflectances, from which bathymetry was subsequently estimated. Then, we quantified known biases in the AC using a linear regression that estimated bias as a function of imaging and environmental factors and applied a correction to produce a new set of remote sensing reflectances. This correction improved bathymetry estimates for eight of the nine scenes we tested, with the resulting changes in bathymetry RMSE ranging from +0.09 m (worse) to −0.48 m (better) for a 1 to 25 m depth range, and from +0.07 m (worse) to −0.46 m (better) for an approximately 1 to 16 m depth range. In addition, we showed that an ensemble approach based on multiple images, with acquisitions ranging from optimal to sub-optimal conditions, can be used to estimate bathymetry with a result that is similar to what can be obtained from the best individual scene. This approach can reduce time spent on the pre-screening and filtering of scenes. The correction method implemented in this study is not a complete solution to the challenge of AC for satellite-derived bathymetry, but it can eliminate the effects of biases inherent to individual AC algorithms and thus improve bathymetry retrieval. It may also be beneficial for use with other AC algorithms and for the estimation of seafloor habitat and water quality products, although further validation in different nearshore waters is required.


Introduction
Bathymetric information from satellite data is of fundamental importance in optically shallow waters, where the seafloor is visible from space and the water-leaving radiance (L w ) is influenced by reflection off the seafloor. Such information, in the form of maps of water depth, is essential for a wide variety of purposes including offshore activities (e.g., pipeline laying), resource management (e.g., fishery), and defense operations (e.g., navigation). Traditional bathymetric charts are based on soundings obtained during hydrographic surveys. However, as ship-borne surveys are costly and an AC inter-comparison exercise. Warren et al. [21] evaluated the accuracy of a wide range of freely available AC processors by comparing them to reference R rs data from different coastal and inland waters. Similarly, in a more recent AC exercise, Zhang and Hu [22] also analyzed an AC algorithm, comparing its R rs images with those measured over a few sites from the AERONET-OC stations. Collectively, these studies demonstrated that accurate AC remains a challenge for OC remote sensing where precise R rs data are needed. Therefore, it is important to explore ways by which errors in AC outputs, and their effect on the products derived from them, can be minimized. One way to address some of the problems posed by imprecise AC is to assess and quantify the impacts of environmental variables on AC accuracy and then account for this in the atmospherically corrected image. In an earlier study [23], four publicly available AC processors (2 land-based and 2 water-based) for deriving the R rs in coastal waters were compared and validated with 54 R rs match-up datasets from AERONET-OC stations. The study revealed that biases from ACOLITE and SeaDAS, two of the state-of-the-art AC algorithms, are influenced by environmental variables. In this study, we demonstrated the potential of Landsat 8 (L8) data for SDB in US coastal waters and assessed the performance of a commonly used and publicly available water-based AC algorithm (ACOLITE [24]) for physics-based SDB. To minimize the effect of imperfect AC on the bathymetry retrieval, we further used a correction factor to improve the original atmospherically corrected image from ACOLITE. Using a set of 9 images, SDB estimates from these two AC procedures were then compared with LiDAR-derived bathymetry of the area. Lastly, we used an ensemble approach to produce SDB of the study area using all the corrected images.

Study Sites
The Florida Keys is a series of islands that extend from the southern end of Florida, USA, to the south-southwest. Their nearshore shallow waters include coral reef tracts, patch reefs, bank reefs, seagrass meadows, and unvegetated hard and soft bottom. This site was chosen because of its relatively clear waters, the good knowledge of seabed features, and availability of LiDAR-derived depth data for validating SDB estimates of water depth. The benthic environment of the section of the Florida Keys used in this study is dominated by extensive seagrass beds, with some patches of reef and unconsolidated sediments. Figure 1 shows this area with the distribution of bathymetric LiDAR data used for validating the SDB estimates in this study.

Satellite Data
Nine L8 images ( Figure 2) from the Florida Keys, acquired during both optimal and near-optimal conditions for SDB, were downloaded from the archive of the United States Geological Survey after visually inspecting all available images from May 2013 to May 2019. L8 OLI (Operational Land Imager) collects visible, Near Infrared (NIR) and Short-wave infrared (SWIR) spectral band imagery at 30 m spatial resolution. In addition to the improved positional accuracy of 14 m, compared to 50 m for its predecessors in the Landsat series, L8 includes coastal and aerosol (433-453 nm) and blue (450-515 nm) bands for coastal and bathymetric mapping [25,26].

LiDAR Data
To validate the SDB estimates, a bathymetry topographic digital elevation model (DEM) was acquired from the National Oceanic and Atmospheric Administration (NOAA) National Centres for Environment Information (NCEI) coastal LiDAR archive. The LiDAR data collection was conducted in December 2014 over South Florida and the Florida Keys as part of efforts by NOAA to study sea level rise and coastal flooding impacts on US coasts. Several LiDAR sources including topographic and bathymetric LiDAR sensors were used to develop and create a suite of tiled bathymetric-topographic DEMs for South Florida and the Florida Keys [27]. A portion of the DEM tiles covering the study site ( Figure 1) was retrieved from the Office of Coastal Managements Data Access Viewer [28] where all DEM data are archived. The DEMs, with a vertical accuracy of approximately 0.5 m, are referenced vertically to the North American Vertical Datum of 1988. Horizontal positions were provided in geographic coordinates and referenced to the North American Datum of 1983 [29,30]. A portion of the collection covering the Florida Keys coastal area was referenced to mean sea level and resampled from 0.3 m to 30 m to match the spatial resolution of L8.

Atmospheric Correction
We implemented two types of AC methods for water depth retrieval: (1) we used ACOLITE to process L8 images into R rs values (henceforth Rrs raw ) and (2) then applied a correction factor to reduce errors in the original ACOLITE output and create new corrected R rs values (henceforth Rrs corrected ). ACOLITE [24], specifically designed for AC over water surfaces, is an AC method that estimates L w by simulating contributions from molecular (Rayleigh) and particulate (aerosol) scattering using a 6SV-based LUT [31]. Based on Ruddick et al. [32], aerosol reflectance is estimated by determining a per-tile aerosol type (or epsilon) from the ratio of reflectances in two bands over water pixels where L w can be assumed to be zero. Then, the epsilon is used to extrapolate the observed aerosol reflectance to the visible bands to remove atmospheric contributions. ACOLITE was originally designed for processing L8 images, but it has been modified and updated to also process Sentinel-2 data [33]. Furthermore, the most recent version, which can be adapted to commercial sensors such as Pleiades, contains an additional AC scheme (now the default setting) called the dark spectrum fitting (DSF) algorithm, as well as a sun glint correction scheme [34]. In this study, ACOLITE (version 20170113.0) was used to produce all R rs images, which are the direct input into the bathymetry algorithm. The default SWIR option (1609 and 2201 nm band combination) was implemented for all images. This band combination takes advantage of the longest-wavelength SWIR band, where water absorption is the highest. In a previous study [23], in which a range of AC algorithms were compared and validated against in situ L w from 14 AERONET-OC stations, statistically significant relationships were demonstrated between errors in ACOLITE's R rs estimates for L8 s 443 nm and 482 nm bands and three environmental variables: Solar Zenith Angle (SZA), Aerosol Optical Thickness (AOT) at 865 nm (AOT 865 ), and wind speed (u 10 ); probable but statistically non-significant relationships were also demonstrated for the 561 nm and 655 nm bands. Using multiple linear regression, we therefore derived a set of coefficients that were used to estimate the error of ACOLITE's R rs estimates for each of those four bands in each image, as a function of SZA, AOT 865 , and wind speed. Then, each of the four bands used for depth retrieval in this study was corrected using Equation (1): where Rrs corrected and Rrs raw are the R rs images with and without correction, respectively; and a, b, c, and d are coefficients obtained through fitting a linear model to the data from Ilori et al. [23]. SZA was obtained from the metadata of each L8 scene. AOT 865 was processed and obtained using the l2gen processor in the SeaDAS software, and an average value used for each image was calculated by randomly sampling multiple pixels over the area of the study site. Wind speed data were obtained from the National Centers for Environmental Prediction Reanalysis project [35], where 6 h global wind speed estimates are archived. Table 1 presents the value of each environmental parameter for each image used in this study.

Sun Glint Correction
As sun glint correction is not inherently part of the ACOLITE version used in this study, we implemented the NIR method [36] to remove specular reflection off the sea surface for images where glint was visually obvious. This method assumes that for optically deep areas (where radiation reflected from the seafloor has a negligible influence on L w ), any remaining NIR signal after AC must be due to sea surface reflection. Thus, glint intensity and removal is performed by establishing a linear relationship between the NIR and visible bands over an optically deep area in the image, and that relationship is then used across all water pixels to reduce R rs for the visible bands to its assumed glint-free value.

Estimation of Noise Equivalent Reflectance
Bathymetry model inversion based on least squares optimization techniques is generally sensitive to environmental noise [37,38]; thus, high environmental noise may make images unsuitable for bathymetry extraction. The noise-equivalent difference in reflectance, NE∆R rs (sr −1 ), is a measure of image noise, with contributions from the sensor (e.g., instrument degradation) and the environment (e.g., variability in atmosphere and water surface state) [37,38]. The NE∆R rs can be used to assess the suitability of a satellite imagery for aquatic remote sensing applications. For example, it has been used to determine the suitability of the Compact Airborne Spectrographic Imager (CASI) for benthic mapping [11]. Therefore, following AC, we estimated the NE∆R rs (sr −1 ) [39] by calculating the band-wise standard deviation of R rs from a 33 × 33-pixel window over a homogeneous optically deep area using Equation (2) [40]. This approach assumes that any observed spectral variations in the selected area is due to noise; thus, selected pixels must be as homogenous as possible for a robust standard deviation estimate. Ideally, the NE∆R rs should be lower than 0.00025 sr −1 in each of the visible bands [41], which was the case for all nine images used in this study. Table 2 shows the per-band value obtained for each of the 9 images used in this study.
where σR rs is the standard deviation in each band over an as homogeneous as possible area of optically deep water within the image.

Parameterization of Environmental Properties
To implement the physics-based approach to SDB, values of optical properties and substratum spectral reflectance that are representative of the environment in question are needed. Water inherent optical properties (IOP) (P 440 , G 440 , and X 550 ) parameterization for forward modeling for each site was based on assessment from Level 3 OC products from the Visible Infrared Imaging Radiometer Suite Visible Infrared Imaging Radiometer Suite (VIIRS) Generalized Inherent Optical Property (GIOP) algorithms [42]. P 440 is the phytoplankton absorption coefficient at 440 nm, G 440 is the absorption of gelbstoff and detrital materials coefficient at 440 nm, and X 550 is the particulate backscattering of suspended particles coefficient at 550 nm. Using parameter values obtained from these OC products, ranges of values for each parameter were determined by observing the lowest and highest parameter values for all dates from GIOVANNI, which is an online visualization tool for OC products [43]. Then, values slightly lower and higher than the observed lowest and highest values, respectively, were chosen ( Table 3). As part of the inversion model, seafloor reflectance spectra are also needed. We used two seafloor spectra (Figure 3), based on the area's benthic description [44]. Depth (Z), which was also needed for forward modeling, was set to 0.1 and 25 m with the understanding that depth penetration greater than 25 m would be difficult.

Forward Modeling of Remote Sensing Reflectance
To derive water depth, we applied a modified version of the semi-analytical inversion model of Lee et al. [17,18] to the atmospherically corrected images. In this inversion scheme, the sub-surface remote sensing reflectance, r rs , (the ratio of upwelling radiance to downwelling irradiance just below the surface) is related to absorption (a) and backscattering properties (b b ) of the water column, the seafloor reflectance (ρ), and water depth (H). For nadir-viewing satellites, the model can be expressed as: where r rs , the subsurface remote-sensing reflectance, is expressed as: where θ w and θ v are the sub-surface solar zenith and sub-surface sensor viewing angles, respectively. Absorption (a) and backscattering coefficients (b b ) are functions of (1) the absorption coefficient of phytoplankton at 440 nm, P; (2) the absorption coefficient of colored dissolved materials at 440 nm, G; and (3) the backscattering coefficient of suspended particles at 550 nm, X. These are expressed as: where a w and b bw are the absorption and backscattering coefficients of pure water, respectively [45], a * phy is the specific absorption of coefficient of phytoplankton (normalized to a value of 1.0 at 440 nm), λ is the center wavelength, and Y is the spectral shape that depends on the particulate shape and size.
While Lee's inversion model uses the albedo of only one key benthic substrate (sand), our model includes a parameterization to set the seafloor reflectance as a linear mix of the two bottom types (i.e., sand and algae; [46]). To forward model the R rs as a function of water depth, water quality parameters, and the seafloor reflectance, the adaptive look-up table (ALUT) method [11,16] was implemented, which ensures efficient construction and search through the table. In this approach, an LUT consisting of the modeled R rs values of L8 bands 1-4, seafloor reflectance (Figure 3), water optical properties (absorption and scattering characteristics of water), and water depths of the optically shallow zone of the area in question (Table 3) is constructed. With realistic minimum and maximum values of all environmental parameters in the table, the LUT construction is optimized by using a hierarchical structure to efficiently cover the range of expected R rs values while minimizing the under-or over-sampling of spectrally similar regions of environmental space, which is common with discretization by regular intervals in conventional forward models. For example, a small change in depth in shallow water areas will cause a significant/large change in R rs, and will typically result in under-sampling if the depth parameter is discretized by regular intervals. Likewise, oversampling may occur in deep water areas where a small change in R rs is expected [11,16]. To identify the best parameter values to be included in the hierarchy in the ALUT technique, discretization is based on an evenly sampled spectral space and not on an evenly sampled parameter space (e.g., water depth). This approach requires bounded ranges for all the modeled parameters, for which we used the value ranges in Table 3.

Inversion of Remote Sensing Reflectance
Model inversion was subsequently performed using the binary space partitioning (BSP) approach [11,16], as described in Knudby et al. [12]. Briefly, this technique subdivides the LUT created during forward modeling into different nodes. First, the BSP splits the whole LUT into two (left and right child nodes) and subsequently subdivides the nodes into a partitioning tree, which facilitates the optimization of the per-pixel LUT search. After model inversion and depth retrieval, water depths were corrected for tidal height at the time of each image acquisition using tidal height estimates obtained from Oregon State University's tide prediction service [47].
Depth cannot be accurately estimated for optically deep pixels (i.e., where reflection off the seafloor is a negligible component of the radiance measured by the sensor); thus, we estimated a depth threshold value for each scene to distinguish between optically deep and optically shallow areas. These threshold values were obtained by calculating the mean minus 2 standard deviations of depth predictions within a 33 × 33 window of homogeneous optically deep water. Results are reported for each scene both a) for the full depth range, and b) for depths from the surface and down to these scene-specific thresholds.

Validation of Depth Estimates
Validation of depth estimates from the two AC procedures was performed by comparing the estimated depths to the LiDAR data. The number of depth estimates used for validation (Table 4) varied between the nine images due to differences in the number of pixels for which depth was successfully estimated, as pixels that did not pass the AC's internal quality checks (e.g., due to clouds), pixels with negative depths, and pixels that were visually impacted by boats, wake, or cloud shadows were eliminated prior to validation. Based on the remaining pixels, we used the coefficient of determination (R 2 ), RMSE (root-mean-squared-error) (Equation (9)) and bias (Equation (10)) to compare the accuracy of the uncorrected and corrected SDB estimates with the LiDAR datasets. The RMSE is used to measure the accuracy of the estimated depth values; and bias is used to indicate overestimation (positive value) or underestimation (negative value): where n is the number of observations, and x est and x obs are the estimated and measured depths, respectively. Values closer to zero for both error metrics indicate a better result. SDB obtained with Rrs raw and Rrs corrected are hereafter referred to as SDB raw and SDB corrected , respectively. These summary statistics (R 2 , RMSE, and bias) were calculated both for the full depth range and for depths ranging from the surface down to the per-scene depth threshold (Table 4).

Result and Discussion
Scatterplots showing water depth estimates produced from both Rrs raw and Rrs corrected images and the LiDAR depth measurements are shown in Figure 4, and summary statistics (R 2 , RMSE, and bias) are listed in Table 4.
The RMSE values for SDB raw /SDB corrected estimates range from 1.35/1.05 m to 2.21/2.04 m for the full depth range, and from 1.21/0.95 m to 2.19/1.89 m when applying the per-scene depth threshold. These values are broadly comparable with other SDB studies (e.g., [10,13,14,41,48]). For example, Dekker et al. [48], who compared one empirical and five physics-based approaches to bathymetry mapping using hyperspectral imagery, reported RMSE values between 0.86 (best) and 4.71 m (worst) for depths less than 13 m for two clear tropical waters in the Bahamas and eastern Australia, suggesting that our results are typical of what should be expected from a physics-based bathymetry method. It is worth mentioning here that impacts from recent hurricanes over parts of the Florida Keys, notably in 2016 and 2017, have resulted in an average increase of approximately +0.3 m in seafloor elevation over different habitat types [49,50]. Such changes were not accounted for in this study and may have had a small effect on the results, although it is worth noting that the best SDB estimate is from 2019 [ Figure 4h], after the hurricanes. Figure 4 shows that accuracy decreases with depth for both SDB raw and SDB corrected , particularly beyond approximately 15 m where the proportion of the measured signal originating from reflection at the seafloor becomes very small. In general, for depths shallower than 15 m, SDB corrected points cluster more tightly around the 1:1 line that do the SDB raw points.

Turbidity
Out of the nine Rrs images we applied the correction factor to, eight corrections resulted in negative RMSE changes when considering the two depth limits used in this study, with reductions ranging from 0.09 to 0.48 m (for the full depth range) and 0.07 to 0.46 m (for the per-scene depth threshold). For both depth limits, only one correction resulted in increased RMSE (i.e., the image from 01/05/2015, see Table 4) (RMSE values for SDB raw and SDB corrected will hereafter be referred to as RMSE raw and RMSE corrected , respectively). For this image, accurate depth estimates were not possible beyond approximately 14 m (Figure 4b), regardless of correction, and RMSE corrected increased marginally by 0.01 m for the full depth range and by 0.04 m for the per-scene depth threshold (i.e., 0-11.23 m). A visual inspection of this image shows sediment plumes in the study area (Figure 5b), which suggests that turbidity contributed to an underestimation of water depth for both SDB raw and SDB corrected [14], and the image is of marginal use for SDB, regardless of correction. Similarly, one of the two corrected images with the lowest RMSE reduction for the full depth range also has what looks to be a silt plume emerging from nearshore channels in the southwestern portion of the area for which depth was calculated (Figure 5a). For this image, RMSE value marginally decreased by 0.09 m for the full depth range and by 0.26 m for the per-scene depth threshold (i.e., 1 to 12.03 m).

Glint
With RMSE raw values of 2.21 m and 2.19 m for the full depth range and the per-scene depth threshold, respectively, the image from 13/01/2018 produces the poorest SDB raw and SDB corrected estimates out of the nine images. As shown in Table 4, this image has the highest RMSE and positive bias values. A visual inspection of this image (Figure 5c) indicates the presence of a moderate glint. Glint correction was not performed, as the image did not show any noticeable improvement after the initial testing. Likewise, for the image from 14/02/2018 (Figure 5d), the high negative bias values for both depth ranges, regardless of correction (Table 4), may be attributed to residual sun glint in addition to light turbidity. While an attempt was made to de-glint this image as described in Section 3.1.2, given that the NIR-based de-glinting method [36] implemented (1) relies on manual selection of deep-water pixels to estimate glint contribution, (2) assumes that there are glint-free pixels among those selected [51], and (3) assumes a homogenously low R rs (NIR) across all water pixels, failure to meet these conditions may have resulted in the observed residual glint. For example, R rs (NIR) may be non-negligible in glint-free but very shallow or turbid waters, or where reflective vegetation such as seagrass is close to the upper water column [52]. For these two images, the correction produces slightly reduced RMSE values (Figure 4e,f), substantially so for their respective depth thresholds (Table 4).

Effect of Wind Speed and SZA on SDB Performance
Out of the eight scenes whose SDB performance improved with the correction for the R rs images, greater corrections were done for five scenes (i.e., Figure 4a,c-f) acquired with SZA > 49 • (Table 1), and one scene (i.e., Figure 4g) acquired during high wind speed (3.21 m/s) ( Table 1), when considering the per-scene depth thresholds. Likewise, when considering the full depth range, greater corrections were also done for the same number of scenes under similar environmental conditions, with the exception of the image from 26 January 2017, whose change in RMSE is −0.09 m (Table 4). It should be noted that the most noticeable RMSE reduction for SDB corrected images for both depth limits considered in this study (i.e., Figure 4d) was observed for the scene with the highest SZA and highest wind speed. This is supporting evidence for the existence of a relationship between ACOLITE's overestimation of R rs in the first two bands of L8 and these two environmental variables [23], and it gives an idea of the magnitude of its impact on SDB performance. A recent study [53] also found a dependency between AC retrieval accuracy and wind speed in coastal waters. While there is strong evidence to conclude that the correction factor used in this study lowers RMSE values for images with high wind speed, it should be noted that wind speed data used in this study come from 6-h estimates in reanalysis model, and thus, they have their own uncertainty. For this reason, more testing may be needed for a firmer conclusion about the relationship between ACOLITE's error and wind speed. Figure 6 shows the performance of SDB raw and SDB corrected for each image, binned to 5 m depth increments for the 1-25 m depth range. The accuracy of SDB estimates decreases with increasing depth for both SDB raw and SDB corrected . While higher RMSE values should be expected at deeper depths due to the diminishing signal from seafloor reflectance, it should be noted that the number of LiDAR points for validation is also smaller at deeper depths, leading to increased uncertainty around the RMSE values reported at these depths.

SDB Estimates Using an Ensemble Approach
Most SDB studies are based on a single image for each study area, with researchers typically selecting the best available image using visual inspection [14]. Our results indicate that this may not be a robust approach. To illustrate the problem, we invite readers to visually inspect the nine images used in this study ( Figure 2) and identify the one that looks most suitable for SDB. Then, proceed to Table 4 to see if it was indeed the one that produced the best results, as measured by RMSE, bias, or R 2 . An informal test among our colleagues, all of whom work on OC remote sensing, suggests that it is not easy to identify the best scene. However, a unique advantage of optical remote sensing is the repetitive acquisition of images over the same area. This is especially important for SDB, where the suitability of a given image is determined by transient environmental factors, such as cloud and aerosols, sea surface state, and turbidity [54]. We explored one way of taking advantage of the multiple images available for the study area by testing an ensemble approach in which we calculated the per-pixel median depth value of all nine corrected images (i.e., SDB corrected ) used in this study. Then, we compared the resulting depth estimates with those obtained using the best individual image from the analysis in Section 4.1 (i.e., the image from 01/02/2019, see Table 4). Figure 7 shows that the results produced by the ensemble are very similar to those obtained with the best individual image. SDB estimates up to approximately 15 m are similar to the SDB corrected estimates from the 01/02/2019 image, as are the RMSE values for the 1-5, 5-10, and 10-15 m depth ranges ( Figure 8). Outliers are noticeably reduced in the ensemble result when compared to any of the sub-optimal images that were also included in its calculation, suggesting that the use of median depth is effective in eliminating noise in the ensemble. Thus, the ensemble approach eliminates the need for the selection of a single best image, while producing SDB results of similar accuracy. In this context, it is noteworthy that one of the two best images (i.e., image from 05/03/2019) is also the one with the lowest NE∆R rs in bands 1, 2, and 4, as well as the second-lowest NE∆R rs in band 3 ( Table 2). This suggests that one effective way to pre-screen images, either for a single best image approach or to determine which images should be included in an ensemble, could be to estimate NE∆R rs and select those images with the lowest values across the visible bands.

Conclusions
In this study, we demonstrated the use of Landsat 8 data for physics-based SDB in US coastal waters. A state-of-the-art AC method (ACOLITE) was used to convert per-pixel radiometric units to R rs , and an RTM was inverted to estimate water depth, which was compared to airborne LiDAR validation data. The results showed that ACOLITE can be used to produce SDB from imagery that is free of conditions such as clouds, glint, sediments plumes, boats and wakes, with an accuracy (RMSE 1.34 m for 1-25 m depth range and 1.21 m for approximately 1-15 m depth range) comparable to that reported from empirical and physics-based SDB elsewhere. To account for ACOLITE's known overestimation of R rs for Landsat 8 s coastal and blue bands, we applied a correction factor, which was calculated as a function of solar zenith angle, aerosol optical thickness, and wind speed, to obtain a corrected set of R rs images. This correction further improved bathymetry estimates for eight of the nine scenes we tested, with the resulting changes in bathymetry RMSE ranging from +0.01 m (worse) to −0.48 m (better) for a 1 to 25 m depth range, and from + 0.04 m (worse) to −0.46 m (better) for an approximately 1-16 m depth range. Using a total of nine Landsat 8 images, we showed that the correction factor improved SDB results, both on average (∆RMSE = −0.22 m) and for the best single image (∆RMSE = −0.30 m) for the 1-25 m depth range. SDB improvements from application of the correction factor were the greatest for images acquired at a high solar zenith angle and at high wind speeds, where ACOLITE is known to have the greatest bias. The correction method demonstrated in this study can be implemented with any appropriate AC algorithm. Finally, we demonstrated that an ensemble approach based on multiple images, with acquisitions ranging from optimal to sub-optimal conditions, can be used to derive bathymetry with a result that is similar to what can be obtained from the best individual image. This is important because it is rarely visually obvious which of several images is best for SDB, and the ensemble approach can be automated to reduce time spent on pre-screening and filtering scenes, and it can potentially also reduce the amount of missing pixels caused by clouds and cloud shadows encountered in any single image. Automating SDB will ultimately facilitate the efficient and operational use of the globally available L8 (and other multispectral) datasets.