Towards Routine Mapping of Shallow Bathymetry in Environments with Variable Turbidity: Contribution of Sentinel-2A/B Satellites Mission

Satellite-Derived Bathymetry (SDB) has significant potential to enhance our knowledge of Earth’s coastal regions. However, SDB still has limitations when applied to the turbid, but optically shallow, nearshore regions that encompass large areas of the world’s coastal zone. Turbid water produces false shoaling in the imagery, constraining SDB for its routine application. This paper provides a framework that enables us to derive valid SDB over moderately turbid environments by using the high revisit time (5-day) of the Sentinel-2A/B twin mission from the Copernicus programme. The proposed methodology incorporates a robust atmospheric correction, a multi-scene compositing method to reduce the impact of turbidity, and a switching model to improve mapping in shallow water. Two study sites in United States are explored due to their varying water transparency conditions. Our results show that the approach yields accurate SDB, with median errors of under 0.5 m for depths 0–13 m when validated with lidar surveys, errors that favorably compare to uses of SDB in clear water. The approach allows for the semi-automated creation of bathymetric maps at 10 m spatial resolution, with manual intervention potentially limited only to the calibration to the absolute SDB. It also returns turbidity data to indicate areas that may still have residual shoaling bias. Because minimal in-situ information is required, this computationally-efficient technique has the potential for automated implementation, allowing rapid and repeated application in more environments than most existing methods, thereby helping with a range of issues in coastal research, management, and navigation.


Introduction
Seafloor mapping plays a pivotal role in using and managing the world's oceans in a way that is in accordance with the United Nations Sustainable Development Goal 14 ("Life below water -conserve and sustainably use the oceans, seas and marine resources") that aims to achieve a better and more sustainable future by 2030 [1]. While bathymetric information is key to the world's management of coastal environments, we still have sparse, outdated, and spatially limited coverage. According to the International Hydrographic Organization (IHO) and the Intergovernmental Oceanographic Commission (IOC), the global bathymetry baseline available to date is surprisingly incomplete: an estimated 70% of the world's coastal seafloor remains unmapped, unobserved or inadequately surveyed to modern standards, and is, therefore, poorly understood [2][3][4][5]. The coastal shallow water zone can be a challenging environment in which to acquire water depth information using conventional methods, such as the vessel-based multi-beam sonar or the active non-imaging airborne lidar. These This paper proposes the first multiple-image technique with Sentinel-2A and Sentinel-2B that allows semi-automated SDB estimation in areas of variable turbidity. Although some recent works have already obtained accurate results with Sentinel-2 and similar satellites, they were implemented in regions with transparent waters [18][19][20]29,30]. We focus on application of the log ratio model [31] because it was designed to support routine mapping in areas with extremely limited calibration data. Two different environments with varying water transparency conditions are examined and the main advantages and challenges of this new method are discussed.

Study Areas
The study areas were Saint Joseph Bay, Florida (29. scales [9,10,22]. This paper proposes the first multiple-image technique with Sentinel-2A and Sentinel-2B that allows semi-automated SDB estimation in areas of variable turbidity. Although some recent works have already obtained accurate results with Sentinel-2 and similar satellites, they were implemented in regions with transparent waters [18][19][20]29,30]. We focus on application of the log ratio model [31] because it was designed to support routine mapping in areas with extremely limited calibration data. Two different environments with varying water transparency conditions are examined and the main advantages and challenges of this new method are discussed.

Study Areas
The study areas were Saint Joseph Bay, Florida (29.75º N, 85.35º W, Figure 1a) and Cape Lookout, North Carolina (34.62º N, 76.54º W) (Figure 1b). Cape Lookout National Seashore preserves a 90-kmlong section of the Southern Outer Banks of North Carolina with three undeveloped barrier islands. St. Joseph Bay is located in Gulf County between Apalachicola and Panama City. The north end of the bay is a relatively narrow opening to the Gulf of Mexico, and is approximately 24 km-long north to south and 10 km-wide at its widest point. The waters of St. Joseph Bay contain the St. Joseph Bay State Buffer Preserve and the St. Joseph Bay Aquatic Preserve.
These areas were selected because they have variable turbidity [32][33][34][35] and the availability of recent airborne lidar bathymetry (ALB) data for final validation and error analysis. In addition, the National Oceanic and Atmospheric Administration (NOAA) was interested in evaluating hurricane induced depth changes: after Hurricane Matthew (October 2016) in Cape Lookout and before Hurricane Irma (September 2017) in Joseph Bay. These study sites also represent complex microtidal environments with variable bottom types.

Lidar Data for Validation
NOAA's National Geodetic Survey (NGS) collected Airborne Lidar Bathymetry (hereinafter ALB) for Cape Lookout in December 2016 and for Joseph Bay in November 2016. These point cloud data sets were both collected using the Riegl VQ-880-G sensor, which provided high-resolution bathymetric data in nearshore waters. The Riegl VQ-880-G used a green laser that operates in a These areas were selected because they have variable turbidity [32][33][34][35] and the availability of recent airborne lidar bathymetry (ALB) data for final validation and error analysis. In addition, the National Oceanic and Atmospheric Administration (NOAA) was interested in evaluating hurricane induced depth changes: after Hurricane Matthew (October 2016) in Cape Lookout and before Hurricane Irma (September 2017) in Joseph Bay. These study sites also represent complex microtidal environments with variable bottom types.

Lidar Data for Validation
NOAA's National Geodetic Survey (NGS) collected Airborne Lidar Bathymetry (hereinafter ALB) for Cape Lookout in December 2016 and for Joseph Bay in November 2016. These point cloud data sets were both collected using the Riegl VQ-880-G sensor, which provided high-resolution bathymetric data Remote Sens. 2020, 12, 451 4 of 23 in nearshore waters. The Riegl VQ-880-G used a green laser that operates in a circular scan pattern, which could penetrate shallow clear water to the seafloor. The high-density point samples were combined with GPS and other positional data to create precise 3D topobathy elevation models. NGS used coastal elevation data to map the mean high-water shoreline, which is considered the nation's official shoreline. These high-resolution observations at 1 m were selected as reference data set in the two study sites and compared to SDB products. The Mean Lower Low Water (MLLW), standard chart datum was used as the reference. The range of depths within this data set is 0-13 m. ALB data referenced to the MLLW was gridded at the Sentinel-2 resolution (10 m) via arithmetic averaging.

Sentinel-2A/B Imagery
Sentinel-2A and 2B twin polar-orbiting satellites, developed by the ESA to meet the operational needs of the Copernicus programme, were used. Both Multispectral Instruments (MSI) on-board are now operational: Sentinel-2A was launched on 23 June 2015 and Sentinel-2B followed on 7 March 2017. The radiometric, spectral and spatial characteristics of the bands used in this study are specified in the User Handbook [36]. Sentinel-2 Level-1C products were downloaded from the Sentinel's Scientific Data Hub [37] and images of zone 18 in Cape Lookout (sub-tile SUD) and of zone 16 in Joseph Bay (sub-tile RFT) were used. A temporal examination provided Level-1C Sentinel-2 images were typically geo-located within two pixels of each other (20 m) which is within the stated quality requirements for absolute geo-location [38]. The study period was selected based on lidar collection and the passage of hurricanes given that intense resuspension and currents may have modified shallow seabed morphology, confounding comparison with the lidar survey. For Cape Lookout, a one-year data series was evaluated from January to December 2017. A total number of 59 Sentinel-2A and 2B scenes were available, but only 15 optimal final images (25%) were further processed due to intense cloud coverage and sunglint effects (Table 1). In St. Joseph Bay, only seven usable Sentinel-2A scenes were available owing to cloud cover and surface reflectance (glint) for a study period of December 2016-March 2017. St. Joseph Bay is on the east side of the swath, as a result, surface glint was severe during spring and summer. Furthermore, no images were considered after Hurricane Irma (September 2017). During that period, only the Sentinel-2A satellite was operational, with a less than 10-day revisit over the Florida zone. Matlab R.2016a software and QGIS (version 3.6.0) were used for visualization and processing of satellite data. Table 1. List of Sentinel-2A and 2B images used in this study for the multi-temporal approach in Cape Lookout (15 scenes) and in Saint Joseph Bay (7 scenes). The percentage of pixels from each scene included in the final pSDBred and pSDBgreen (Equations (1) and (2), Section 2.5) after the multi-scene approach is indicated. The pixels evaluated (N) were only pixels with corresponding lidar data for final validation: N = 271635 in Cape Lookout and N = 678517 in Joseph Bay.

Atmospheric Correction
Sentinel-2 images were processed to Level-2A (L2A) with the ACOLITE processor developed by the Royal Belgian Institute of Natural Sciences (RBINS), which supports free processing, specifically for aquatic applications, of both Landsat-8 and Sentinel-2 [39][40][41][42]. ACOLITE products corresponded to Remote Sensing Reflectance (Rrs, 1/sr) in all visible and Near-Infrared (NIR) bands and chlorophyll concentration by the OC3 algorithm (Chl, mg/m 3 ) resampled to 10 m pixel size. We selected a combination of the NIR (865 nm) and Short Wave Infrared (SWIR) (1600 nm) bands in ACOLITE for the aerosol correction with a user defined epsilon value (maritime aerosol=1). This strategy has been shown to significantly improve the quality of the products by minimizing the influence of NIR/SWIR instrument noise [43]. A spatial filter (median filter 3x3) was conducted on the bands in order to remove some noise and inter-pixel variability. Recent experience with Sentinel-2 imagery indicated spatial filtering enhanced bathymetric products [10,22,44]. In this investigation, the spectral red-edge band at 704 nm (Rrs704) was used as a proxy for turbidity over optically shallow waters [22]. Several researchers have already indicated red-edge bands were appropriate for turbidity or suspended solids monitoring in optically shallow regions [43,45].

Satellite-Derived Bathymetry Model
In this study, the ratio model of log-transformed bands having different water absorption was applied (Equations (1) and (2)). The model was designed to support routine mapping in clear waters with extremely limited calibration data [31]. The model uses the Remote Sensing Reflectance (Rrs, units of sr −1 ) of the blue (490 nm), green (560 nm) and red (664 nm) bands for each satellite image corrected for atmospheric effects (Equations (1) and (2)), and the log-transform addresses the exponential decrease in light with depth [31]. In this case, we used the ratio of blue (λ i ) to either green or red (λ j ) bands to produce the SDB (hereinafter called SDBgreen and SDBred, respectively). SDBgreen performs better in deep areas, and SDBred performs better in shallow water [22][23][24]44]: where, pSDB = ln(n πRrs(λ i )) ln n πRrs λ j pSDB is the relative or "pseudo" depth from satellite (dimensionless), SDB is the satellite-derived depth (meters), Rrs is the Remote Sensing Reflectance, m 1 and m 0 are tunable constants to linearly transform the model results to actual depth, π has units of sr, and n=1000 is a fixed constant to assure that both logarithms will be positive under any condition and that a residual non-linearity in the ratio is removed from depths that are retrievable from satellite [31].

Multi-Scene Approach
Thanks to the 5-day revisit of the Sentinel-2 twin satellites, the mission can offer imagery that would identify transient turbidity features that produce false shoals. Accordingly, in this study, we used Remote Sens. 2020, 12, 451 6 of 23 a multi-scene approach (see Figure 2). The pSDBgreen and pSDBred were determined, per Equation (2), for all available scenes (Table 1). It was reasonable to assume that the contribution from the temporal variability in water-column turbidity to SDB models was greater than the temporal variability of the seafloor [46]. As turbidity produces a false shoaling of depth, all scenes were compared to identify the maximum pSDB at each pixel from all the scenes, and two resultant composite images were created: one of the maximum pSDBred, and one of the maximum pSDBgreen, with assumption being maximum pSDB values would most likely be without any shoal effect. Composite images of Rrs704 and chlorophyll from the OC3 algorithm were also returned, where the turbidity value at each pixel came from the same scene as the depth for that pixel. Both Rrs704 and chl were surrogates for addressing the false shallowing, not intended to describe actual turbidity or chl concentration.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 23 the temporal variability in water-column turbidity to SDB models was greater than the temporal variability of the seafloor [46]. As turbidity produces a false shoaling of depth, all scenes were compared to identify the maximum pSDB at each pixel from all the scenes, and two resultant composite images were created: one of the maximum pSDBred, and one of the maximum pSDBgreen, with assumption being maximum pSDB values would most likely be without any shoal effect. Composite images of Rrs704 and chlorophyll from the OC3 algorithm were also returned, where the turbidity value at each pixel came from the same scene as the depth for that pixel. Both Rrs704 and chl were surrogates for addressing the false shallowing, not intended to describe actual turbidity or chl concentration.

Vertical Referencing with Chart Soundings
After generating the final pSDBred and pSDBgreen maps, the parameters m1 and m0 (Equation 1) were tuned by linear regression with about ten control points from charts in each of the study

Vertical Referencing with Chart Soundings
After generating the final pSDBred and pSDBgreen maps, the parameters m 1 and m 0 (Equation (1)) were tuned by linear regression with about ten control points from charts in each of the study regions and for pSDBred and pSDBgreen separately. Depth measurements (also known, as soundings) from NOAA charts [47] were used to select the control points over Cape Lookout (11545) and Joseph Bay (11389). Using chart points keeps the calibration independent of the lidar used for validation, and corresponds to methods that would be utilized in typical application [31]. The points were chosen as areas of uniform depth that were less likely to have changed over time (regions away from inlets and sand waves). The selected calibration points (pixels) for each of the four composite images were found to originate in different input scenes (3 and 4 different scenes respectively for SDBred and SDBgreen for St-Joseph Bay, and 5 for both Cape Lookout composites). The vertical calibration inherently corrects for the reference data by shifting the depths to the tidal datum, which in the case of USA is the MLLW, through the tuning of the coefficient m 0 . Sea level depends on satellite observation time; in this case, both regions are microtidal environments. For this approach, our hypothesis was that the influence of the tide for a pSDB relative to each scene is small compared to the influence of the false shallowing generated by turbidity from multiple satellite images. The validation of this hypothesis and the improvement of accuracy by compensating for tide effects in meso and macrotidal regimes are subjects for future study.
In addition, to test the robustness of the multi-scene approach and the final SDB result in terms of a common calibration for remote unsurveyed locations, we interchanged the coefficients for pSDBred and pSDBgreen between the two regions, thus using the coefficients of Cape Lookout in St. Joseph Bay, and vice versa.

Switching Model
The final step for a corrected bathymetry mapping (SDB) corresponded to a switching model implementation between SDBred and SDBgreen due to the sensitivity of each model: while SDBred performs better over shallow regions [48], SDBgreen performs better over deeper regions [21,31,44,49]. Moreover, SDBgreen frequently yields severe overestimation in shallow regions [22,44,50], especially with dense seagrass, highly undesirable for navigational purposes. The switch between models used the conditions detailed below: SDBred < 2 m, SDB=SDBred SDBred > 2 m and SDBgreen > 3.5 m, SDB=SDBgreen SDBred >= 2 m and SDBgreen <= 3.5 m, SDB=SDBweighted SDBweighted was determined by a simple linear weighting calculation to account for a smooth transition (Equation (3)): where the depth weighting (for 3.5 m and 2 m) is determined by: The final SDB map was generated in each region then compared to the lidar surveys for validation. Assessment of the discrepancy between SDB and lidar used the mean difference as the bias metric, the Median Absolute Error (MedAE) providing the typical total error, and the interquartile range (IQR) as a measure of statistical dispersion. These are robust metrics that do not require an assumption of a Gaussian error distribution.

Submerged and Floating Aquatic Vegetation Masking
An additional step was carried out in St. Joseph Bay due to the existence of floating aquatic vegetation (FAV) and shallowed submerged aquatic vegetation (SAV). The complex coastal waters of Saint Joseph Bay are characterized by several benthic types, including high density submerged and floating aquatic vegetation at depths shallower than 2 m [51]. In order to correct this issue, we established a masking strategy for the FAV and SAV area. The MCI (Maximum Chlorophyll Index) first designed for the MEdium Resolution Imaging Spectrometer (MERIS), measuring the radiance peak at the red-edge band (709 nm) in water leaving radiance, indicates the presence of a high surface concentration of chlorophyll against a scattering background [52]. This index was formed with three MERIS narrow channels centered near 681, 709 and 754 nm used to define a linear baseline. The MCI has been used to map floating vegetation or benthic [53,54]. In this study, given that Sentinel-2 has three bands (665, 704, and 740 nm) similar to the MERIS bands used for MCI, we utilized Equation (4).
FAV and SAV areas were established for the pixels with positive values of MCI (MCI>0). This MCI mask was applied to the final SDB map after the multi-scene approach and switching model in St. Joseph Bay in order to remove possible overestimation for shallow depths where vegetation was present. Therefore, the MCI was not used to map seagrass, but to locate and remove an error source with SDB.

Cape Lookout in North Carolina
Prior to implementing the multi-temporal approach with Sentinel-2A/B images in Cape Lookout (Figure 1b), we estimated bathymetry for several single scenes in order to evaluate the impact of varying turbidity.  Figure 3h). These outcomes demonstrated the impact of variable water quality conditions on the accuracy of the predicted SDB. Although the 23 January scene had lower turbidity overall-and might be selected as the optimal scene for analysis-there are still areas of elevated turbidity evident offshore (see the lower right of 3e and 3f). Over complex turbid environments such as Cape Lookout [32,34], using only a single scene for SDB has definite limitations. October 2017, the color bars indicate pixel frequency at the point; residual errors (SDBgreen -ALB) for (c) 23 January (standard deviation=2.07m, percentile 5%/95% =-4.06m/2.93m), and (d) 30 October (standard deviation=3.17m, percentile 5%/95% =-6.88m/3.16m), the bins of the histograms correspond to the number of pixels; maps of SDBgreen for (e) 23 January and (g) 30 October; and maps of the turbidity proxy used in this study, Remote Sensing Reflectance at 704 nm (Rrs704), for (f) 23 January and (h) 30 October. Gray color represents land mask and white color at the left is the limit of Sentinel-2 tile for each date.
The compositing selected the pixels of the multiple scenes that were slightly affected by turbidity. Compositing multiple scenes to identify the pixels in order to correct for water turbidity substantially improved the accuracy, as can be observed in Figure 4 and Figure 5. The assumption of the approach is based on the submarine terrain remaining slightly unchanged during the study period (1 year) to accumulate a time series of scenes to use in the compositing, while turbidity (and noise: waves, cloud shadows, ships, and bubbles) affect the accuracy of bathymetric inversion. Using the multi-scene pixel compositing approach with the sixteen images (Table 1), and then applying the switching model between SDBred and SDBgreen, gave better performance from shallow to deeper regions. The composite solution scheme efficiently described the nearshore depth with minimum scatter and MedAE=0.41 m over the range of ALB data (up to 12.5 m) compared to the individual scenes (Figure 4a  The compositing selected the pixels of the multiple scenes that were slightly affected by turbidity. Compositing multiple scenes to identify the pixels in order to correct for water turbidity substantially improved the accuracy, as can be observed in Figures 4 and 5. The assumption of the approach is based on the submarine terrain remaining slightly unchanged during the study period (1 year) to accumulate a time series of scenes to use in the compositing, while turbidity (and noise: waves, cloud shadows, ships, and bubbles) affect the accuracy of bathymetric inversion. Using the multi-scene pixel compositing approach with the sixteen images (Table 1)       In addition, the switching strategy reduced error over using only SDBred (Figure 6a,c) or SDBgreen (Figure 6b,d). Recently, the switching algorithm has been suggested to be an opportune method for mapping regions from very shallow to deeper waters [22,44]. The quality of the combined models is shown in the histogram of the residual errors, which has a symmetric and narrow distribution (Figure 4b). The vertical calibration results are presented in Figure 7. Compared to the current standard approaches focused on the selection of the optimal scenes; in this case, the clearest image was acquired on 27 November 2017, compositing reduced the mean error by 50% and the IQR from 1.6 m to 0.9 m (Figure 4a,b vs. Figure 4c,d).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 23 In addition, the switching strategy reduced error over using only SDBred (Figure 6a and 6c) or SDBgreen (Figure 6b and 6d). Recently, the switching algorithm has been suggested to be an opportune method for mapping regions from very shallow to deeper waters [22,44]. The quality of the combined models is shown in the histogram of the residual errors, which has a symmetric and narrow distribution (Figure 4b). The vertical calibration results are presented in Figure 7. Compared to the current standard approaches focused on the selection of the optimal scenes; in this case, the clearest image was acquired on 27 November 2017, compositing reduced the mean error by 50% and the IQR from 1.6 m to 0.9 m (Figure 4a,b vs. Figure 4c,d).   The method was evaluated by taking into account the selection of images, confirming that pixels from multiple scenes with high turbidity conditions (Rrs704) and chlorophyll (chl) were minimally incorporated into the model (<1%). Table 1 detailed the percentage of each scene incorporated into the final model. Turbidity (defined as either Rrs704 or chl) produces a false shoaling because the relative or "pseudo" depth pSDB (before the vertical calibration) decreases as turbidity increases (see an example for Cape Lookout in Figure 8). There is a relationship between the increase in the two water quality parameters and the decrease in relative depths for both models, pSDBred and pSDBgreen, thus confirming that intense false shoaling patterns were associated with highly turbid water. Similar results have been found in other sites along the Florida coastal waters using imagery from Sentinel-2A and Sentinel-3A [22]. Figure 5d represents the turbidity proxy (Rrs704) associated with the final SDB map, providing an indication of areas where residual turbidity may exhibit potential residual biases. The method was evaluated by taking into account the selection of images, confirming that pixels from multiple scenes with high turbidity conditions (Rrs704) and chlorophyll (chl) were minimally incorporated into the model (<1%). Table 1 detailed the percentage of each scene incorporated into the final model. Turbidity (defined as either Rrs704 or chl) produces a false shoaling because the relative or "pseudo" depth pSDB (before the vertical calibration) decreases as turbidity increases (see an example for Cape Lookout in Figure 8). There is a relationship between the increase in the two water quality parameters and the decrease in relative depths for both models, pSDBred and pSDBgreen, thus confirming that intense false shoaling patterns were associated with highly turbid water. Similar results have been found in other sites along the Florida coastal waters using imagery from Sentinel-2A and Sentinel-3A [22]. Figure 5d represents the turbidity proxy (Rrs704) associated with the final SDB map, providing an indication of areas where residual turbidity may exhibit potential residual biases.  (Table S1). The red and green circles indicate examples of the maximum pSDBred and pSDBgreen, respectively, selected for the multi-scene approach.
The method also exhibited some advantages over the in-situ data. The lidar survey was limited by water clarity as well, and some water areas nearshore were not retrieved (Figure 5c). In contrast, by using multiple dates of imagery, features could be identified by the composited satellite data, along the entire coastal fringe and within the inner banks (SDB on Figure 5b compared to ALB on Figure 5c). We also examined the switching model used in this study, so the final SDB map ( Figure  5b) was split into the SDBred and SDBweighted (0-3.5m, Figure 5e), and the SDBgreen (>3.5m, Figure  5f). SDBred and SDBweighted mapped the shallow regions within the banks and barrier islands, while the SDBgreen mapped the offshore and the deep channels in the banks. In offshore areas where the bottom was not visible, the attributed depth was removed with a cut-off depth of 17 m. These results were support by visual comparison with the RGB composite ( Figure 5a) and standard chart soundings collected by the National Oceanic and Atmospheric Administration [47], given that there was not high-resolution information available for validation of those deeper water areas.

Saint Joseph Bay in Florida
Previous studies have characterized the Saint Joseph Bay and the adjacent area as a moderately turbid region under influence of the Apalachicola Bay turbid waters [33,35]. The region also has areas  (Table 1). The red and green circles indicate examples of the maximum pSDBred and pSDBgreen, respectively, selected for the multi-scene approach.
The method also exhibited some advantages over the in-situ data. The lidar survey was limited by water clarity as well, and some water areas nearshore were not retrieved (Figure 5c). In contrast, by using multiple dates of imagery, features could be identified by the composited satellite data, along the entire coastal fringe and within the inner banks (SDB on Figure 5b compared to ALB on Figure 5c). We also examined the switching model used in this study, so the final SDB map (Figure 5b) was split into the SDBred and SDBweighted (0-3.5 m, Figure 5e), and the SDBgreen (>3.5 m, Figure 5f). SDBred and SDBweighted mapped the shallow regions within the banks and barrier islands, while the SDBgreen mapped the offshore and the deep channels in the banks. In offshore areas where the bottom was not visible, the attributed depth was removed with a cut-off depth of 17 m. These results were support by visual comparison with the RGB composite ( Figure 5a) and standard chart soundings collected by the National Oceanic and Atmospheric Administration [47], given that there was not high-resolution information available for validation of those deeper water areas.

Saint Joseph Bay in Florida
Previous studies have characterized the Saint Joseph Bay and the adjacent area as a moderately turbid region under influence of the Apalachicola Bay turbid waters [33,35]. The region also has areas of dense shallow aquatic vegetation [51], so a single scene was used to evaluate model performance over these areas.  In order to correct the overestimation issue, a masking procedure to locate and remove an error source based on a common algae index-the Maximum Chlorophyll Index MCI [53,54]-was developed. In [54], it is shown that slightly submerged vegetation will cause an MCI (false positive against phytoplankton). This approach allowed us to use the satellite without requiring any additional in-situ information. Pixels with positive MCI values (MCI>0) were identified as vegetation ( Figure 5d) and then masked out. The scatterplots of ALB against SDBred and SDBgreen after masking vegetation (Figures 9f and 9c, respectively) show that the overestimation was eliminated, thus MCI can be applied to remove the overestimation. Very shallow SDBred pixels with good accuracy (lying on the 1:1 line) were also eliminated, because MCI does register on bright sand in very shallow water (owing to the signal from Rrs704 in <1 m). As a conservative solution, the FAV and SAV masking procedure would be appropriate owing to the severe consequences of overestimating very shallow bathymetry, as identified by the Chart Standards Groups (CSG) at NOAA's Office of Coast Survey/Marine Chart Division. In order to correct the overestimation issue, a masking procedure to locate and remove an error source based on a common algae index-the Maximum Chlorophyll Index MCI [53,54]-was developed. In [54], it is shown that slightly submerged vegetation will cause an MCI (false positive against phytoplankton). This approach allowed us to use the satellite without requiring any additional in-situ information. Pixels with positive MCI values (MCI>0) were identified as vegetation (Figure 5d) and then masked out. The scatterplots of ALB against SDBred and SDBgreen after masking vegetation (Figure 9c,f, respectively) show that the overestimation was eliminated, thus MCI can be applied to remove the overestimation. Very shallow SDBred pixels with good accuracy (lying on the 1:1 line) were also eliminated, because MCI does register on bright sand in very shallow water (owing to the signal from Rrs704 in <1 m). As a conservative solution, the FAV and SAV masking procedure would be appropriate owing to the severe consequences of overestimating very shallow bathymetry, as identified by the Chart Standards Groups (CSG) at NOAA's Office of Coast Survey/Marine Chart Division.
The same multi-scene approach and switching model as used for Cape Lookout was carried out in this region, although only seven cloud-free images were available during the study period (Table 1).
In this case, we tested the additional MCI masking to remove issues associated with FAV and SAV. The validation against lidar data confirmed the high accuracy obtained for the entire strategy with (Figure 10a) or without (Figure 10c) MCI masking. There was low scatter (MedAE=0.3 m) for depths up to 11 m (limit of ALB data set). The vertical calibration results are presented in Figure 7. Even without MCI correction, the switching model alone substantially reduced intense overestimation of SDBgreen for depths <1 m (Figure 10c). The importance of the switching strategy is evident in the validation of the final multi-scene SDBred and SDBgreen, where an acute overestimation occurred at depths <1 m for SDBgreen, similar to Cape Lookout (Figure 6a-d). In terms of error assessment, the distribution of the residuals indicated minimum discrepancies were achieved (Figure 10b,d). Furthermore, the compositing algorithm produced accurate SDB results compared to the single "optimal scene" approach (picking the scene with the lowest overall turbidity) combined with the most widely used model in the literature, SDBgreen (Figure 10e,f). The composite map of the proxy for turbidity (Rrs704) associated to the final SDB model shows minimum turbid conditions throughout (Figure 11d).
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 The same multi-scene approach and switching model as used for Cape Lookout was carried out in this region, although only seven cloud-free images were available during the study period (Table  1). In this case, we tested the additional MCI masking to remove issues associated with FAV and SAV. The validation against lidar data confirmed the high accuracy obtained for the entire strategy with (Figure 10a) or without (Figure 10c) MCI masking. There was low scatter (MedAE=0.3 m) for depths up to 11 m (limit of ALB data set). The vertical calibration results are presented in Figure 7. Even without MCI correction, the switching model alone substantially reduced intense overestimation of SDBgreen for depths <1 m (Figure 10c). The importance of the switching strategy is evident in the validation of the final multi-scene SDBred and SDBgreen, where an acute overestimation occurred at depths <1 m for SDBgreen, similar to Cape Lookout (Figure 6a-d). In terms of error assessment, the distribution of the residuals indicated minimum discrepancies were achieved (Figure 10b,d). Furthermore, the compositing algorithm produced accurate SDB results compared to the single "optimal scene" approach (picking the scene with the lowest overall turbidity) combined with the most widely used model in the literature, SDBgreen (Figure 10e-f). The composite map of the proxy for turbidity (Rrs704) associated to the final SDB model shows minimum turbid conditions throughout (Figure 11d).   The geographic distribution of features in the SDB (Figure 11b) corresponded to those identifiable in the lidar survey (Figure 11c), although much more spatial information is available in the SDB product. Similar to Cape Lookout region, the method offered a complete map of the St. Joseph coastal area compared to the restricted map provided by lidar surveys. In this case, we did not apply the MCI masking to add data over the vegetated area (depths <1 m), thereby allowing maximum coverage of the final SDB product. Shallow areas and shoals < 3.5 m were accurately described with SDBred and SDBweighted (Fig 11e) whereas depths >3.5 m were mapped with SDBgreen (Figure 11f).

Discussion
In the present investigation, we successfully applied both Sentinel-2A and 2B satellites to estimate water depth in two areas with complex bathymetry and water clarity, Cape Lookout, North Carolina, and Saint Joseph Bay, Florida. The multi-scene compositing approach addressed limitations The geographic distribution of features in the SDB (Figure 11b) corresponded to those identifiable in the lidar survey (Figure 11c), although much more spatial information is available in the SDB product. Similar to Cape Lookout region, the method offered a complete map of the St. Joseph coastal area compared to the restricted map provided by lidar surveys. In this case, we did not apply the MCI masking to add data over the vegetated area (depths <1 m), thereby allowing maximum coverage of the final SDB product. Shallow areas and shoals < 3.5 m were accurately described with SDBred and SDBweighted (Figure 11e) whereas depths >3.5 m were mapped with SDBgreen (Figure 11f).

Discussion
In the present investigation, we successfully applied both Sentinel-2A and 2B satellites to estimate water depth in two areas with complex bathymetry and water clarity, Cape Lookout, North Carolina, and Saint Joseph Bay, Florida. The multi-scene compositing approach addressed limitations inherent in conventional methods and reduced the impact of turbidity, performing better than the standard "pick the best scene" method that relies on a single image (Figures 3 and 4c vs. Figure 4a in Cape Lookout; and Figure 10e vs. Figure 10a,c in Saint Joseph Bay). The final corrected SDB produced robust depths up to the limit of the lidar surveys, with typical errors ≤0.4 m. These excellent results from Sentinel-2 compared favorably with those produced in relatively low turbidity water in south Florida [22,44], and in regions with transparent waters [10,18,20,30]. Whereas some researchers suggested there is still work to be performed regarding the identification of the optimal period throughout the year where bathymetric errors are minimized [18,29], others asked for novel strategies to allow seabed mapping without the laborious analysis per image and the visual inspection of the "clearest scene" [19]. Recent studies have already indicated the potential of multi-scene approaches in order to select the optimal scene or eliminate noise over clear waters [16,17,19,20]. However, our temporal compositing strategy successfully reduced the turbidity impact without requirement of visual inspection, thereby enhancing SDB performance in an easy way. The high temporal resolution and interchangeability of the Sentinel-2 twin mission may rapidly overcome SDB anomalies introduced by highly heterogeneous water transparency regimes [22].
The multi-scene strategy applied here did not require any screening or manual adjustment of the imagery prior to compositing. It automatically picked the pixels least impacted by turbidity (e.g., Figure 8) from the set of scenes provided (Table 1), substantially simplifying the effort compared with other studies where the selection of optimal images with variation in the water quality conditions were essential in the extraction of SDB [20,22,29,55]. Manual selection of the optimal scene is not only highly subjective, but requires considerable time and effort, and may still include regions having patches of turbidity. The consistency in the pSDB products, which were not yet calibrated to true depth, indicated that ACOLITE produced an effective and robust atmospheric correction across scenes (Figure 8), as previously demonstrated [22,44]. Inconsistencies or errors in ACOLITE (or any other atmospheric correction) would force manual intervention and make compositing impractical [10].
We chose chart values for the calibration to depth in order to demonstrate the benefit of the method for remote areas under likely applied conditions, rather than under the optimal (and extremely unlikely) condition of a contemporaneous lidar survey. Using limited chart soundings would have only introduced some of the error, most likely as bias in the depth. Some depth errors may also be due to the differences in resolution, i.e., 10 m from the satellite, and the 1 m spot size for lidar-a consideration in areas with steep gradients, such as channel edges. Using a median filter on the SDB reduced other artifacts that typically lead to random noise [49]. One potential error factor that needs to be considered in the future is tide range. The tide range in these areas is relatively small (<1 m). The influence of tide on the accuracy of this approach for areas with tidal ranges greater than 1 m requires further investigation. For meso-tidal areas, one option would be to constrain the input images to a common range of water level, such as 1 m maximum difference in water level. Macro-tidal areas are likely to be too turbid for optical SDB.
An interesting consideration is that the calibration to depth for this method may be dependent on water clarity rather than on geography. Applying the coefficients for composited pSDBred and pSDBgreen of Cape Lookout to St. Joseph Bay (Figure 7), and vice versa, generated precise SDB for both areas with errors ≤0.65 m (see results in Figure 12). This comparison suggests that an existing calibration after the multi-scene approach may be applicable to other remote and unsurveyed areas, assuming the compositing can retrieve the same turbidity. These results suggest the potential suitability of SDB estimates for charting applications using the IHO S-57 standard [56], which defines CATZOC levels containing depth accuracy specifications for depth ranges up to 30 m in order to allow for their incorporation into nautical products.  Figure 8.
The use of the SDBred for shallow water retrieval addressed the worst of the overestimation issues with the SDBgreen in both study regions. The use of different bands to treat overestimation in shallow water (or underestimation by SDBred in deep water) has been pointed out recently by some researchers [19,30,57]. Switching to SDBred has been previously used to provide better discrimination of depths in very shallow water over bright targets like carbonate sand [48], and other studies have applied each model over different depth ranges [22,44]. Here, we established an automated switching method between SDBred and SDBgreen that performed accurately over the two study sites. In addition, SDBgreen overestimated on dense seagrass in extremely shallow water using either the single scene (Figure 9b, Figure 10e) or the multi-scene approach (dense seagrass has a very low albedo, so in deeper water there is often no detectable signal to be retrieved from satellite). By switching to the SDBred, the combined product reduced most of the impact of shallow seagrass on SDBgreen (as seen in Figure 9b against Figure 9e or Figure 10e against Figure 10c). If more rigorous masking is required, the MCI can identify the presence of dense seagrass in extremely shallow water (<1-2 m, Figure 10a), without the need for external data sets on bottom type. Combined with calibration using only a few points (Figure 7), the approach detailed here offers an effective means of assessing bathymetry in regions lacking any data.
While the accuracy of SDB cannot match ALB, the spatial coverage of satellite imagery shown here surpasses ALB (Figure 5c and Figure 11c vs. Figure 5b and Figure 11b). This capability represents one of the great benefits of satellite monitoring; whereas in-situ surveys (lidar, sonar) are expensive and extremely limited due to technical and deployment cost, SDB can give us a broad picture of the entire study site with wide swaths, low-cost repeated coverage, and easy access to remote locations [8]. Lidar data is also constrained by availability and water quality, so information can be lost. We characterized post Hurricane Matthew bathymetry in Cape Lookout, extending the limited amount of lidar data to a much larger and more complete area.
Sentinel-2 may provide a way to use other satellites more effectively. Accordingly, Sentinel-2 might offer a low turbidity reference data set that can be used to inform mapping conducted with the various commercial very high-resolution sensors such as the World-View fleet [14,26,58] or the novel CubeSats from Planet [19,29]. On the other scale, Landsat has a four-decade history using Thematic Mapper, and has been shown to be useful for some SDB in clear water [10,21,59,60], even with 30 m The use of the SDBred for shallow water retrieval addressed the worst of the overestimation issues with the SDBgreen in both study regions. The use of different bands to treat overestimation in shallow water (or underestimation by SDBred in deep water) has been pointed out recently by some researchers [19,30,57]. Switching to SDBred has been previously used to provide better discrimination of depths in very shallow water over bright targets like carbonate sand [48], and other studies have applied each model over different depth ranges [22,44]. Here, we established an automated switching method between SDBred and SDBgreen that performed accurately over the two study sites. In addition, SDBgreen overestimated on dense seagrass in extremely shallow water using either the single scene ( Figure 9b, Figure 10e) or the multi-scene approach (dense seagrass has a very low albedo, so in deeper water there is often no detectable signal to be retrieved from satellite). By switching to the SDBred, the combined product reduced most of the impact of shallow seagrass on SDBgreen (as seen in Figure 9b against Figure 9e or Figure 10e against Figure 10c). If more rigorous masking is required, the MCI can identify the presence of dense seagrass in extremely shallow water (<1-2 m, Figure 10a), without the need for external data sets on bottom type. Combined with calibration using only a few points (Figure 7), the approach detailed here offers an effective means of assessing bathymetry in regions lacking any data.
While the accuracy of SDB cannot match ALB, the spatial coverage of satellite imagery shown here surpasses ALB (Figures 5c and 11c vs. Figures 5b and 11b). This capability represents one of the great benefits of satellite monitoring; whereas in-situ surveys (lidar, sonar) are expensive and extremely limited due to technical and deployment cost, SDB can give us a broad picture of the entire study site with wide swaths, low-cost repeated coverage, and easy access to remote locations [8]. Lidar data is also constrained by availability and water quality, so information can be lost. We characterized post Hurricane Matthew bathymetry in Cape Lookout, extending the limited amount of lidar data to a much larger and more complete area.
Sentinel-2 may provide a way to use other satellites more effectively. Accordingly, Sentinel-2 might offer a low turbidity reference data set that can be used to inform mapping conducted with the various commercial very high-resolution sensors such as the World-View fleet [14,26,58] or the novel CubeSats from Planet [19,29]. On the other scale, Landsat has a four-decade history using Thematic Mapper, and has been shown to be useful for some SDB in clear water [10,21,59,60], even with 30 m pixels. Comparisons of Sentinel-2 with Landsat-8 may ultimately lead to results that could expand the utility of the entire Landsat data record for change detection.
In this study, interpretation of error assessment and uncertainty can be achieved by means of the water quality composites of turbidity (Rrs704) associated with the final SDB maps (Figures 5d and 11d). These products may indicate areas that still have residual bias. They can also provide areas where the water is chronically optically deep due to turbidity, which can lead to better allocation of resources. Previous research has detailed the impact of water quality parameters on SDB results using the ratio model over several regions in Florida [22]. Recent studies using the ratio model showed that bottom reflectance has little influence on the accuracy of SDB estimates, while chlorophyll concentration has a strong influence [22,61]. This problem would also be resolved with the compositing method.
In accordance with our findings, the portability and reliability application of the proposed approach using minimum in-situ measurements is a distinct advantage in effectiveness and may especially benefit developing countries. In addition, using the semi-automated solution and a computer cloud-based system may allow exploitation of the Sentinel-2 data for regional to global scale coastal SDB [16]. The average time required for processing a Sentinel-2 image for bathymetry estimation is 1 hour. For further avenues of research with Sentinel-2, we intend to upscale the study results to more environments with different water and atmospheric conditions in order to evaluate repeatability as well as implement the multi-temporal approach in cloud-based computing platforms such as the ESA Coastal Thematic Exploitation Platform-Coastal-TEP [62]. The promising transferability of this temporal technique exploiting the open and free archive of the Sentinel-2 mission will allow advancement of SDB applications and optimized scientific coastal mapping worldwide, especially in data-poor regions [63].

Concluding Remarks
This paper presents a framework to obtain SDB while potentially automated the image processing. Precise SDB was derived over moderately turbid environments by using the high revisit time (5-day) of the Sentinel-2A/B twin mission from the Copernicus programme. The proposed methodology incorporates a robust atmospheric correction with ACOLITE, a multi-scene compositing method to reduce the impact of turbidity, and a switching model to improve mapping in shallow water. Two study sites in United States were explored due to their varying water transparency conditions. The method establishes a plausible high-quality strategy for SDB with sufficient sensitivity to support interannual and pre-post hazards (e.g., hurricanes or tropical storms), which may benefit change detection analysis from multiple images to discriminate variability. The open data policy and long-term mission commitment of Sentinel-2 opens future promising time series evaluation over years and even decades that can be an important tool to provide crucial missing information on the bathymetry distribution, especially in data-poor or remote areas with large gaps in a retrospective, rapid and non-intrusive manner. It would enable immediate commencement of coastal mapping using archived and/or forthcoming imagery from the Copernicus programme at regional to global spatial scales, and any other source from which repeated scenes may be collected over the world's coastal systems. There is a requirement to promote the development of technology that provides low-cost solutions with enhanced data quality to address the needs of a wide range of user groups. The semi-automated framework demonstrated here may address that need, and help mitigate the long-standing gap of depth information for the majority of coastal regions on Earth. Such a capability can aid scientists, managers, and policymakers around the globe in assessing vulnerabilities to ecosystems, infrastructure, navigation and many other coastal concerns.
Author Contributions: I.C. and R.P.S. conceived and designed the research. I.C. collected the remote sensing data, processed the data, conducted the analysis, and prepared the figures. R.P.S. supervised the research. I.C. led the writing of the manuscript with contribution from R.P.S. All authors have read and agreed to the published version of the manuscript.