Evaluation of Satellite-Derived Bathymetry from High and Medium-Resolution Sensors Using Empirical Methods

: This study evaluates the accuracy of bathymetric maps generated from multispectral satellite datasets acquired from different multispectral sensors, namely the Worldview 2, PlanetScope, and the Sentinel 2, in the bay of Elounda in Crete. Image pre-processing steps were implemented before the use of the three empirical methods for estimating bathymetry. A dedicated correction and median ﬁlter have been applied to minimize noise from the sun glint and the sea waves. Due to the spectral complexity of the selected study area, statistical correlation with different numbers of bands was applied. The analysis indicated that blue and green bands obtained the best results with higher accuracy. Then, three empirical models, namely the Single Band Linear Algorithm, the Multiband Linear Algorithm, and the Ratio Transform Algorithm, were applied to the three multispectral images. Bathymetric and error distribution maps were created and used for the error assessment of results. The accuracy of the bathymetric maps estimated from different empirical models is compared with on-site Single beam Echo Sounder measurements. The most accurate bathymetric maps were obtained using the WorldView 2 and the empirical model of the Ratio Transform algorithm, with the RMSE reaching 1.01 m.


Introduction
Bathymetry measurements in oceans, rivers, or lakes are essential, especially in coastal areas with intense use of the coastal zone, heavy sea traffic, and vulnerable natural ecosystems. Monitoring of coastal areas is, thus, of great importance to implement sustainable coastal development and ecosystem protection strategies [1][2][3]. High spatiotemporal resolution and a vertical accuracy topographic and bathymetric data are also essential not only for understanding coastal systems evolution [2], but also for other environmental applications, such as benthic habitat mapping [4], seabed geomorphology [5], underwater archaeology [6], monitoring of coastal morphological changes, navigation, and fishing [7].
Traditional methods for estimating sea-bottom surface include Single Beam Echo Sound (SBES) and Multi-Beam Echo Sound (MBES) installed on boats, as well as active sensor systems-LIDAR installed on aerial platforms, remotely controlled vehicles (RVs) and autonomous underwater vehicles (AUVs) [4]. These are commonly used for highresolution bathymetry retrieval in nearshore areas [8,9]. Although these methods can provide high accuracy of bathymetrical data, they have specific limitations when it comes to cost and time constraints compared to the extensive coverage [10]. In contrast, satellite remote sensing techniques can provide a cost-effective solution to monitor large coastal areas, especially in remote areas, as they can provide a continuous and frequent updating set of data [11]. Indeed, the recent advancement of satellite remote sensors has generated new space-based methodologies for estimating bathymetry, the so-called Satellite-Derived Bathymetry (SDB) [12][13][14][15].
Nowadays, SDB is one of the leading research areas of Remote Sensing Science, studying the maritime environment. Associated products are used in many applications, such as the movement of deposited sediments, the bathymetry monitoring, etc. [11]. Both active and passive sensors can be used for estimating bathymetry. However, passive sensors are nowadays the most widespread ones concerning the use of multispectral and hyperspectral sensors. Over time, several studies were carried out using numerous sensors, techniques, and algorithms to estimate bathymetry. The estimation of bathymetry using passive sensors is based on the theory that as depth increases, while the intensity of electromagnetic energy gradually weakens due to the Inherent Optical Properties (IOPs) of the water column. The passive sensors provide a wide wavelength range of zones with predefined models to achieve the bathymetry estimation. However, the visible bands penetrate the water column differently, i.e., blue band up to 25 m, green band up to 15 m, and red band up to 5 m, based on respective studies using Landsat images [16].
SDB appeared in the late 1970s, where the extraction of the bathymetric substrate was achieved in clearwater, using a linear algorithm of multi bands, such as Lyzenga's approach [17,18]. This approach was followed by other researchers after modifications to refine the estimation of the bathymetry [19][20][21]. Moreover, semi-analytical approaches were applied using models and methods of deriving bathymetric data [22,23]. Stumpf extended the approach of Lyzenga, proposing an empirical method for estimating SDB known as the ratio transform algorithm, which contends that the ratio for two bands at a constant depth will be the same, independently from the substrate of the seabed [24]. Additionally, significant potentials of estimating bathymetry were presented by researchers with stereo techniques [25,26] and more recently with machine learning approaches [27][28][29].
Thus, fewer constraints are needed with empirical methods [30], where the associated models are mathematical equations studying the relationship between distance data, e.g., reflectance, radiance, and digital numbers, with shallow waters. Masita in [31] summarized twelve different empirical models for estimating bathymetry, including the principal component analysis (PCA), linear ratio, multiple linear regression polynomial of ratio transform, least squares boosting fitting ensemble, and support vector regression, among others. An additional recent review applied a comparative analysis for optical SDB within and between the methods with shallow water depth. This research indicated that some of the limitations of bathymetry that pose a significant challenge in assessing bathymetry, such as turbidity, chlorophyll, etc., were not resolved in recent years based on the SDB literature [32].
Data quality and pre-processing steps comprise a critical part of the SDB. Special attention is needed regarding the atmospheric effects, where studies have shown that a precise atmospheric correction increased the accuracy of the estimated depth using SDB [33,34]. An additional factor concerns the sun reflection and the ripple at the sea surface, confounding the properties of the water column and the sea bottoms in remote sensing. A further pre-processing step for SDB estimation is the sun glint corrections, where many studies evaluated techniques and bathymetry results, showing significant improvement in accuracy with sun glint correction [35,36].
The current paper focuses on estimating bathymetry concerning the bay of Elounda located at the NE part of Crete Island, covering an area of about 10 km 2 . It is a semi-enclosed bay with an eastward-facing entrance on the north side. The bathymetric conditions of the bay are considered challenging in estimating the bathymetry using satellite images, as the bay is shallow and covered with fine sediments. However, several studies were carried out to extract bathymetry by applying empirical methods with similar conditions. Casal G. et al. carried out a study of SDB in two bays of Ireland, the Dublin Bay and the bay of Mulroy. This study applied two empirical methods, the Linear band model and Ratio Transformed Algorithm using Sentinel 2 satellite images. The best bathymetrical results were estimated in Gulf of Mulroy using the Linear band model with the value correlation coefficient (R 2 ) set to 0.89 and the Root Mean Square Error (RMSE) 0.78 m, while using Ratio Transformed Algorithm the R 2 reached 0.83 with the RMSE at 0.98 m [37]. One more study was carried out in the bay Tralee of Ireland, where bathymetry was estimated using empirical models with three different sensors (Landsat 8, RapidEye, and Pleiades) [38].
The framework of this study is to estimate bathymetry using diverse spatial resolution of satellite images on a closed bay instead of an open sea. The Worldview 2, PlanetScope, and Sentinel 2 satellite data were used to extract bathymetrical data in the bay of Elounda in Crete. The present study is composed of six sections. After to the introduction, Section 2 mentions the in situ data of Single-Beam Echo Sounders and the satellite images. Then, the pre-processing steps, the three applied empirical models, and the methodology used in this study are described. The next section concentrates on the results of the applied methodology providing the bathymetrical maps, comparing the results for each method and satellite image. The extracted bathymetrical maps combined with the error distribution maps and the turbidity maps will help analyze the errors concerning turbidity. Finally, the paper ends with a discussion and the conclusions.

Study Area
The study site is the bay of Elounda, located northeast of the island of Crete in Greece. The site is located between 35 • 15 25 -35 • 17 7 N and 25 • 43 17 -25 • 44 35 E (WGS 84) as can be seen in Figure 1. The broader area was an important center and a strategic port of the island of Crete in ancient periods. According to the archaeological data, the area exhibited a continuous habitation spanning from the Minoan to Venetian period. The most relevant architectural constructions, however, were created in the Hellenistic, Roman, and Early Byzantine periods. In this area, the waters are characterized by suspended sediment due to fine grain sediments, hydrodynamic conditions, and marine traffic. Consequently, a high volume of suspended sediments regularly noted is an important factor that can alter the results in the estimation of satellite bathymetry.

In Situ Depth Data
In situ bathymetric data were obtained using hydroacoustic sounding and the GPS/RTK satellite positioning technique, with the use of a single beam echo sounder (Hitarget hd370) combined with two TOPCON RTK receivers using the Real Time Kinematic In this area, the waters are characterized by suspended sediment due to fine grain sediments, hydrodynamic conditions, and marine traffic. Consequently, a high volume of suspended sediments regularly noted is an important factor that can alter the results in the estimation of satellite bathymetry.

In Situ Depth Data
In situ bathymetric data were obtained using hydroacoustic sounding and the GPS/RTK satellite positioning technique, with the use of a single beam echo sounder (Hi-target hd370) combined with two TOPCON RTK receivers using the Real Time Kinematic technique. The local RTK/DGPS base station was set up only for the dedicated project. One served as base station and the other was in the vertical pole together with the echosounder. Before the initiation of bathymetric surveys, local geodetic control points were established that included an RTK/DGPS reference station. The error in the field measurements was less than 0.1 m for the position (x, y), while for the depth (z), the error ranged from 0.1 m in shallow to 0.2 m in deep waters, but without considering the pitch, roll, and yaw. Additionally, it is worth noting that wave height corrections were not made for this package of in situ data.
To increase the accuracy of the depth measurement conducted using single beam echo sounder, tests of water parameters were conducted to determine the sound propagation speed in water, which included conductivity and temperature with the use of an ONSER HOBO sensor. For the sea level corrections, a base station of water level measurements was established by deploying another ONSER HOBO sensor that was able to measure the water level, including tide oscillations and temperature per minute. Postprocessing of the measurements by considering all the above information formed the database for the development of the Elounda bay bathymetry.
For this study, 7500 in situ measurements ranging from 0 to 8.50 m depth were collected from the campaigns. After filtering and denoising the data, this number was estimated at 5665 points. Since the satellite sensors refer to a spatial resolution of 2 × 2 or 3 × 3 or 10 × 10 m areas, an extrapolation of the points, based on their mean value per pixel, was estimated. A similar approach has been already reported in the literature [34,38]. This has been processed for each bathymetric equations and method. These data were used for: (a) the calibration of bathymetric results with absolute depths and (b) the estimation of errors. It is worth noting that 20% of the points (i.e., 1145) were kept as validation points. Figure 2 shows the distribution of the calibration (green color) and validation points (red color) for the area of study. It is worth noting that 20% of the points (i.e., 1145) were kept as validation points. Figure 2 shows the distribution of the calibration (green color) and validation points (red color) for the area of study.

Optical Satellite Images
The current study aims to compare SDB maps of different optical sensor spatial resolutions. Three optical satellite images were used to analyze and estimate bathymetry. These images comprise three high-resolution Worldview 2 and PlanetScope images and one medium-resolution Sentinel 2 satellite image.

Optical Satellite Images
The current study aims to compare SDB maps of different optical sensor spatial resolutions. Three optical satellite images were used to analyze and estimate bathymetry.
These images comprise three high-resolution Worldview 2 and PlanetScope images and one medium-resolution Sentinel 2 satellite image.

Methodology
In contrast to the traditional depth measurement methods, satellite bathymetry is a "passive" technology and measures the intensity of solar radiation. The results from satellite bathymetry are influenced by many uncontrollable environmental factors that reinforce the uncertainties of the results. Limitations are that the accuracy of calculating the bathymetry with passive receivers is affected by the water quality (Turbidity), cloud cover, weather conditions, and sun glint.
Seven basic steps of estimating bathymetry were acknowledged to this challenging study area. The pre-processing steps are a crucial component of the methodology for estimating bathymetry. The first pre-processing step concerns the geometric, radiometric, and atmospheric correction of satellite images. The provider of the Worldview 2 and PlanetScope images gave calibrated the satellite images, while the medium resolution Sentinel 2 image was downloaded as the bottom of the atmosphere product. In addition, for the Sentinel-2 image, the resampling of all Sentinel 2 bands using a referenced band with a 10 m spatial resolution was applied in the Sentinel Application Platform (SNAP) software ( Figure 3-i).

Results
This section presents the overall results generated from the different methods and datasets described above. The Equations (4)-(6) were implemented to all different multispectral datasets. Due to the specificity of the study area, the combinations of all visible bands with all algorithms mentioned above were analyzed using the three optical images. This step was implemented because, as Vahtmäe and Kutser have stated, in turbid waters, the optimal bands that are shifted to higher wavelengths, such as the green-yellow spectral zones and not to the blue zone, have greater penetration in the water column of clear waters [44].
A table was created showing the correlation coefficients of in situ depth data between the Single Band Linear Algorithm (SBLA), Multi-Band Linear Algorithm (MBLA), and Ratio Transform Algorithm (RTA), using all visible bands for each optical sensor. Table 2 aims to select the best bands for each algorithm for the best derived bathymetrical results. The second step (Figure 3-ii) concerns the correct distribution of calibration/validation of in situ data. The in situ depth was divided into two packages, in situ depth data for calibration purposes (80%) and in situ depth data for validation purposes (20%).
The third step (Figure 3-iii) was the sun glint correction to all images. Sun glint can be directly reflected on the sensor due to the water's surface and the sun's position often displayed as scattered bright spots or white stripes along the edge of the wave. A significant obstacle to the estimation of satellite bathymetry is that the sun glint effect may present an error in SDB maps of up to 30% in high-resolution images [39]. The algorithm of Hedley was used in this study (Equation (1)). This technique is based on sun glint elimination, which assumes that the radiation of infrared bands (NIR) in deep water is negligible. Consequently, the NIR bands indicate the amount of glint in the received signal.
where R i is the sun glint-corrected pixel, R i is the pixel value in all bands, b i is the regression slope, and R N IR − Min N IR is the difference between the NIR pixel (infrared radiation) value and the minimum pixel value in the infrared. As this study focuses on the watery parts of the study area, masks were created. The masks focus mainly on land/sea separation, cloud removal, and cloud shading. The cloud removal and shading masks were not required and, therefore, not created. The Normalized Difference Water Index (NDWI) was used for land/sea separation (Figure 3-iv). This index uses the green and NIR band to enhance only the water features. Equation (2) shows the NDWI used in the fourth step: where B green is the green band and B N IR represents the band values of the reflected Near-Infrared Radiation. Turbidity can influence the accuracy of estimated bathymetry from satellite images, since it can produce suspended sediments that can contribute to the higher radiation in the water during the visible and near-infrared (NIR) bands. The depths in the shallow water are overestimated and underestimated in the deeper regions [37,40]. Lacaux in 2007 have developed the Normalized Difference Turbidity Index (NDTI), which can quantify the levels of turbidity in coastal waters [41]. Subsequently, using the ArcGIS Pro environment, nine error distribution maps were created by linearly interpolating the difference of in situ depth data and satellite bathymetry's depth values. Equation (3) shows the NDTI: where B green and B red represent the band values of the reflected green and red band, respectively. Then, a median filter (window) of 3 × 3 was applied (Figure 3-v) to reduce the noise. The median filter is often used to remove "sensing" noise, interrupted pixels and other false pixel features while maintaining the overall quality. The main disadvantage of other filters (low pass) instead of median filters is that only blur is corrected instead of removing noise [42].
This research used the most common methods (Ratio Transform Algorithm and Linear Method) to produce depth using satellite images (Figure 3-vi). Lyzenga followed the fundamental principle derived from the Beer-Lambert law. Lyzenga developed a linear method which assumes that the reflection at the bottom is a linear function of the reflectance of the seabed with the exponential function of the depth of water [5,6,30,43]. The estimation of the depth from a single band depends on the albedo and can be made using the following equation: where Z SDB is the estimated depth from optical images, α 0 and b o are constants, and L(λ i ) is the sun glint-corrected value of a band. However, by assuming that the bottom reflection ratio between two spectral bands is constant for all albedo types, Lyzenga showed that multi-bands could provide a correction for albedo in finding the depth, using the equation below: where Z SDB is the estimated depth from optical images, L(l ι ) is the sun glint-corrected value of a band (e.g., green), L∞(li) is the deepwater sun glint-corrected value of a band for spectral band l ι , α 0 and b 0 are constants, and N(i = 0, 1, . . . , N) is the number of spectral bands used.
Subsequently, Stumpf in 2003 [24] stated that the difference in solar attenuation for two zones could be used to extract water depth. This method argues that the ratio for two zones at a constant depth will be the same, regardless of the difference in the shade of the bottom and can be calculated using actual depths on clear shallow water. In this method, the approach of Beer's law was accepted, which considers that the attenuation of light in the water column increases exponentiation as the depth increases [30]. Thus, the ratio transform algorithm is related to the logarithm of the two-zone reflectance and the actual depths. The bands used in the transformation algorithm to estimate bathymetry are the blue, green, red, and infrared zones. Penetration into the blue and green spectrum is higher, while the absorption of electromagnetic radiation increases in the red part of the spectrum. The ratio conversion algorithm can be applied with zones that have different water absorption and can be applied to appropriate wavelengths of any sensor. Because the blue and green zones have a lower absorption, the ratio of the two zones remains the same despite the different shading of the bottom at a constant depth. The red and infrared channels are used to separate land and sea. The values of these zones were applied to the following equation to estimate bathymetry: where Z SDB is the depth of bathymetry, m0 is a constant displacement value of the results at a depth of 0 m, mi is the depths from in situ measurements for calibration purposes, n is a fixed constant, and R(li) and R(lj) are the zones li and lj, respectively. The following step (Figure 3-viii) is the creation of the error distribution maps using the validation points and the extracted bathymetry results. The error distribution maps evaluated the obtained SDB errors by comparing it with the NDTI maps. In this step, the absolute errors are calculated using the equation below: where E absolute is the absolute error, Z SDB is the satellite-derived bathymetry depth, and Z IDD is the in situ depth data. The final step (Figure 3-ix) is the comparison of the accuracy of the used algorithms, estimating the correlation coefficient R 2 and the assessment of the statistical parameter Root Mean Square Error (RMSE): where n is the number of the field points, Z SDB is the satellite-derived bathymetry depth, and Z IDD is the in situ depth data. Figure 3 presents the basic steps of SDB using optical sensors for this study.

Results
This section presents the overall results generated from the different methods and datasets described above. The Equations (4)-(6) were implemented to all different multispectral datasets. Due to the specificity of the study area, the combinations of all visible bands with all algorithms mentioned above were analyzed using the three optical images. This step was implemented because, as Vahtmäe and Kutser have stated, in turbid waters, the optimal bands that are shifted to higher wavelengths, such as the green-yellow spectral zones and not to the blue zone, have greater penetration in the water column of clear waters [44].
A table was created showing the correlation coefficients of in situ depth data between the Single Band Linear Algorithm (SBLA), Multi-Band Linear Algorithm (MBLA), and Ratio Transform Algorithm (RTA), using all visible bands for each optical sensor. Table 2 aims to select the best bands for each algorithm for the best derived bathymetrical results.   Table 2 states the correlation coefficients using the three algorithms for all satellite images. Comparing the correlation coefficients for each algorithm and satellite image noted the low values with red cells, mean levels with yellow cells, and high values with green cells. As observed in Table 2, the correlation coefficient with the highest value reaches up to 0.74. This observation is mainly due to the conditions of the study area (high suspended sediments), which is a critical component to the accuracy of SDBs. The highest correlation coefficients value in SBLA was observed in the green band for all three sensors. However, using the MBLA, the best fit values of the correlation coefficients differ for each image. The best fit value of R 2 was noted in the sum of the logarithm of the Green/Yellow band for WV2 sensor, Coastal/Green in Sentinel 2 image, and Green/Red in PlanetScope image using MBLA. Finally, using RTA resulted in the highest correlation coefficients value in Green/Blue bands for all three sensors. Subsequently, the bathymetric maps were created using the SBLA, MBLA, and RTA equations and bands with the highest correlation coefficients, noted in Table 2. Figure 4 presents the results of SDB. Table 2 states the correlation coefficients using the three algorithms for all satellite images. Comparing the correlation coefficients for each algorithm and satellite image noted the low values with red cells, mean levels with yellow cells, and high values with green cells. As observed in Table 2, the correlation coefficient with the highest value reaches up to 0.74. This observation is mainly due to the conditions of the study area (high suspended sediments), which is a critical component to the accuracy of SDBs. The highest correlation coefficients value in SBLA was observed in the green band for all three sensors. However, using the MBLA, the best fit values of the correlation coefficients differ for each image. The best fit value of R 2 was noted in the sum of the logarithm of the Green/Yellow band for WV2 sensor, Coastal/Green in Sentinel 2 image, and Green/Red in PlanetScope image using MBLA. Finally, using RTA resulted in the highest correlation coefficients value in Green/Blue bands for all three sensors. Subsequently, the bathymetric maps were created using the SBLA, MBLA, and RTA equations and bands with the highest correlation coefficients, noted in Table 2. Figure 4 presents the results of SDB.
A regression analysis was performed to estimate the correspondence of the satellite bathymetric outcomes against the in situ data. Each bathymetrical map, validation points, regression analyses, and the coefficient of determination (R 2 ) were calculated for the different satellite images. The overall results are presented in Figure 5. Nine scatterplots were designed with the x-axis representing the validation points, while the y-axis represented the estimated bathymetry points derived from the Worldview 2, Sentinel 2, and PlanetScope images. A regression analysis was performed to estimate the correspondence of the satellite bathymetric outcomes against the in situ data. Each bathymetrical map, validation points, regression analyses, and the coefficient of determination (R 2 ) were calculated for the different satellite images. The overall results are presented in Figure 5. Nine scatterplots were designed with the x-axis representing the validation points, while the y-axis represented the estimated bathymetry points derived from the Worldview 2, Sentinel 2, and Plan-etScope images.  In addition, the RMSE, using Equation 7, was calculated for each pair. The overall results are shown in Table 3. According to these results, the highest correlation is observed using the RTA for Worldview 2 and Sentinel 2 and MBLA for PlanetScope. The bathymetric map derived by WV2 using the Ratio Transform Algorithm marked the best results of RMSE to be 1.01 m. The highest value of RMSE can be found on Sentinel 2 images using the MBLA algorithm. For PlanetScope, the variation of the RMSE in all three empirical models is relatively small, and the RMSE ranges from 1.08 to 1.13 m for all three algorithms, with the best results observed in the MBLA.  In addition, the RMSE, using Equation (7), was calculated for each pair. The overall results are shown in Table 3. According to these results, the highest correlation is observed using the RTA for Worldview 2 and Sentinel 2 and MBLA for PlanetScope. The bathymetric map derived by WV2 using the Ratio Transform Algorithm marked the best results of RMSE to be 1.01 m. The highest value of RMSE can be found on Sentinel 2 images using the MBLA algorithm. For PlanetScope, the variation of the RMSE in all three empirical models is relatively small, and the RMSE ranges from 1.08 to 1.13 m for all three algorithms, with the best results observed in the MBLA. The following results present the areas of water turbidity and evaluate the bathymetrical results using the validation points and the error distribution maps. Consequently, the areas with turbidity for each satellite image were identified and classified using the NDTI indices. The water areas with extreme turbidity levels may have more radiation response of red light than green light. A prediction of an indicator of water quality was determined by classifying the NDTI values into five equal intervals categories. The results applying the NDTI indices for each satellite image are presented in Figure 6. The following results present the areas of water turbidity and evaluate the bathymetrical results using the validation points and the error distribution maps. Consequently, the areas with turbidity for each satellite image were identified and classified using the NDTI indices. The water areas with extreme turbidity levels may have more radiation response of red light than green light. A prediction of an indicator of water quality was determined by classifying the NDTI values into five equal intervals categories. The results applying the NDTI indices for each satellite image are presented in Figure 6.
The red color shows the areas with the highest NDTI values, the yellow areas with high NDTI values, the green areas for moderate NDTI values, light blue areas with low NDTI values, and blue areas with the lowest NDTI values. The results of the waters' turbidity in each satellite image differs as the sensing dates of the satellite images are different. For example, in the shallow waters in the satellite image of Sentinel 2, the high turbidity of the water was noted northeast of the bay compared to the other two, which had a medium to high level of turbidity. This may be due to the conditions prevailing in the bay of Elounda at the acquisition date.
Subsequently, based on the SDB results, error distribution maps were created to study the depth and areas with the highest errors presented on bathymetric maps. Subsequently, using the ArcGIS Pro environment, nine error distribution maps were created by linear interpolation of the difference of in situ depth data and satellite bathymetry's depth The red color shows the areas with the highest NDTI values, the yellow areas with high NDTI values, the green areas for moderate NDTI values, light blue areas with low NDTI values, and blue areas with the lowest NDTI values.
The results of the waters' turbidity in each satellite image differs as the sensing dates of the satellite images are different. For example, in the shallow waters in the satellite image of Sentinel 2, the high turbidity of the water was noted northeast of the bay compared to the other two, which had a medium to high level of turbidity. This may be due to the conditions prevailing in the bay of Elounda at the acquisition date.
Subsequently, based on the SDB results, error distribution maps were created to study the depth and areas with the highest errors presented on bathymetric maps. Subsequently, using the ArcGIS Pro environment, nine error distribution maps were created by linear interpolation of the difference of in situ depth data and satellite bathymetry's depth values, as shown in Figure 7. Each error distribution map refers to the use of a different algorithm/sensor.

Discussion
In the current study, three different satellite sensors, with the acquisition date being the same month, were used to estimate the bathymetry. The study area concerns the bay of Elounda, located northeast of the island of Crete in Greece. A relatively large number of calibration and validation points used in this study allowed us to assign the complexity

Discussion
In the current study, three different satellite sensors, with the acquisition date being the same month, were used to estimate the bathymetry. The study area concerns the bay of Elounda, located northeast of the island of Crete in Greece. A relatively large number of calibration and validation points used in this study allowed us to assign the complexity of the model better and be utilized in the regression models with a satisfactory assessment performance of the models [45]. During the analysis, the overall results show an error ranging from 1.01 m to 1.52 m. The best SDB accuracy was obtained using the RTA on the Worldview 2 sensor. The correlation coefficient reached 0.76 and the RMSE at 1.01 m.
However, the correlation coefficient (R 2 ) in all algorithms and satellite images is not high, and this can be due to several sources of errors. These possible sources of error may be led to the effect of the significant differences in the accuracy of the bathymetric maps. The source of error concerns the weather conditions and the activities carried out in the study area, as well as the specific date of the satellite image. The three algorithms operate by assuming that the water column and the substrate do not alter the radiance values. All three algorithms applied are based on the principle of exponential attenuation of depth, which considers that the depth in clear waters above a homogeneous substrate can be calculated using a single band by solving the Beer-Lambert law. However, high suspension of the fine-grained sediments of the sea floor and the turbidity observed seasonally depending on hydrodynamic conditions but also marine traffic increases the error in the estimation of satellite bathymetry. Furthermore, another possible source of error is the atmospheric corrections made to the satellite images used. Accurate atmospheric correction is essential and has a significant impact on the bathymetry results [39,46].
The error distribution maps show significant systematic errors appearing in the center, the deepest point of the south part of the bay, with 8.50 m depth. Additionally, systematic noise was observed in the shallow water in all models with all sensors. As identified in the bathymetric results with the RTA algorithm, we observed that the errors are visibly less in shallow and deeper areas of the bay. This can be linked with weather conditions or navigation prevailing in the study area at the specific time of the satellite sensing data or sensor resolution.
Concerning the study area's conditions of satellite images, images without clouds and winds are ideal to avoid phenomena of the sea surface reflections. These phenomena could not be avoided in satellite images where a median filter and sun glint corrections were used, which were necessary to estimate bathymetry. In addition to the intense activity of the study area, the bay waters are not clear because of turbidity and suspended fine-grain sediments. These phenomena reduce the accuracy of bathymetry results with diverse spatial resolutions.
Many factors could alter the bathymetric results, running into three different satellite images with different resolution and sensing dates. In this research, basic corrections were made to estimate bathymetry without considering the dynamic nature of the confounding variable, such as turbidity.
The next step of this research was the evaluation of the bathymetric results considering the three methods and turbidity. Subsequently, a comparison of the error about turbidity and depth for each method and satellite was achieved by creating graphs, as shown in Figure 8. The diagrams were created, distributed by two meters depth using the validation points and extracting the errors of Error Distribution Maps and the Turbidity maps. The x-axis of the diagram presents the type of satellite image, and the y-axis presents the values of mean errors by each method (meters) and the level of turbidity. Τhe diagrams of Figure 8 present the variation of turbidity levels in relation to the depth and the satellite images. This may be due to the activities that may take place at the sensing time, e.g., ship routes or even wind waves. It can be noticed that errors of all methods and satellite images have less than one meter for depths from four to six meters.
Comparing the diagrams separately, we noted importance results regarding the turbidity and errors for the depth per two meters. The first diagram in Figure 8 refers to the depth varying between 0-2 m. The turbidity level of the PlaneScope images has the lowest value, with the errors being smaller in contrast to the other two satellites and the three methods. Comparing Worldview 2 and Sentinel 2 results, Worldview 2 has a higher turbidity level with smaller errors using the methods of LDM, MBLA, and higher errors using the RTA method than Sentinel 2. The diagram with a depth of 0-2 m observed better results, with a lower turbidity level than the results derived with the highest resolution image. The diagrams of Figure 8 present the variation of turbidity levels in relation to the depth and the satellite images. This may be due to the activities that may take place at the sensing time, e.g., ship routes or even wind waves. It can be noticed that errors of all methods and satellite images have less than one meter for depths from four to six meters.
Comparing the diagrams separately, we noted importance results regarding the turbidity and errors for the depth per two meters. The first diagram in Figure 8 refers to the depth varying between 0-2 m. The turbidity level of the PlaneScope images has the lowest value, with the errors being smaller in contrast to the other two satellites and the three methods. Comparing Worldview 2 and Sentinel 2 results, Worldview 2 has a higher turbidity level with smaller errors using the methods of LDM, MBLA, and higher errors using the RTA method than Sentinel 2. The diagram with a depth of 0-2 m observed better results, with a lower turbidity level than the results derived with the highest resolution image.
The next diagram refers to the depth results between 2-4 m and is presented in Figure 8. This diagram shows smaller errors at the PlanetScope for all three methods and the smallest turbidity level. The other two satellite images results are noteworthy as the turbidity level is the same. The diagram shows smaller errors for all three methods on WordView 2 in contrast to Sentinel 2. This may mainly be due to the resolution of the satellite image.
The same turbidity level is also spotted for all three satellite images for a depth of 4-6 m (Figure 8). This diagram shows smaller errors in the bathymetric results for all methods with higher resolution images. More specifically, WorldView 2 errors are smaller than PlanetScope's, while the errors of PlanetScope are smaller than Sentinel 2 for all three methods.
In continuation, the results for depth 6-8 m are shown in diagram (d) of Figure 8. The lowest errors for all methods are observed in WorldView 2, having the lowest turbidity level. Comparing the PlanetScope and Sentinel 2 results, the same turbidity level is spotted. The LBM and MBLA methods have a smaller error derived with Planetscope's image, while the results derived from the method RTA are slightly better in Sentinel 2.
In the last diagram (Figure 8), that presents the results for a depth of 8-10 m, a significant error is observed for all the methods for all three satellite images. However, the lowest turbidity's level is spotted on WorldView 2, giving better results only for the methods of MBLA and RTA, as similar results are derived from the LBM method for all three satellite images.
To sum up, turbidity could affect bathymetric results regardless of the resolution of the satellite image. Additionally, smaller errors were noted in the bathymetric maps derived with higher resolution images, comparing areas with the same turbidity level. Moreover, concerning the sensor spatial resolution, similar studies show better bathymetry results with high-resolution sensors [47,48]. The main reason is that low-resolution satellite image causes mixed pixels, because of the decrease reducing the accuracy of bathymetric maps, mainly on the coastal waters. That explains the higher RMSE shown in the Sentinel 2 images than the WV 2 and PlanetScope satellite images.

Conclusions
This paper compares the derived bathymetric results obtained with three different satellite sensors using three different algorithms. The study was carried out expecting differences in the accuracy of SDB calculation due to the resolution of the satellite's sensors and the different sensing dates. The visible bands were evaluated, applying to the three empirical methods and all satellite images, due to the turbidity observing in the study area. This step was implemented to select the best bands for each algorithm for the best derived bathymetrical results, with the green and blue to better respond to all three methods. The pre-processing steps were followed in addition to the sun glint corrections, and median filters were applied. Subsequently, empirical methods were applied to satellite images to estimate bathymetry in the study area. Regarding Table 3 and the RMSE results, it was noted that the bathymetrical models of the study area have significant errors. Better accuracy was reached using the RTA for WV2 and Sentinel 2 images. The PlanetScope images noted better accuracy with the MBLA. Expressly, the lowest RMSE on the Worldview 2 image indicated 1.01 m (RTA), while this was 1.06 m for Sentinel 2 (RTA) and 1.08 m for PlanetScope (MBLA). As the Sentinel's 2 SDB errors were lower than PlanetScope's, a further analysis was performed, producing the error distribution and turbidity maps.
The creation of error distribution and turbidity maps helped statistically analyze the error in relation to depth and turbidity. The results show that satellite images with higher resolution have better accuracy in the bathymetric maps in areas with the same depth and the same level of turbidity. However, significant errors in bathymetric maps were observed in the same areas with different turbidity levels, regardless of the image's resolution for all three methods. This evaluation and comparison of SDB maps show that water quality plays an essential role in in-depth estimation using empirical methods and multispectral images of different spatial resolutions. Future work regarding the sensitivity of the spatial resolution of satellite sensors against in situ bathymetric campaigns is foreseen by the authors. Additionally, the authors will further elaborate on the impact of spatial resolution and turbidity with bathymetric accuracy. Funding: The compilation of this paper was supported by the project 'EXCELSIOR': ERATOS-THENES: EXcellence Research Centre for Earth Surveillance and Space-Based Monitoring of the Environment H2020 Widespread Teaming project (www.excelsior2020.eu, accessed on 18 January 2022). The 'EXCELSIOR' project has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No 857510 and from the Government of the Republic of Cyprus through the Directorate General for the European Programmes, Coordination and Development. Moreover, the compilation of this paper was also supported by the project "DIATOPO-Cretan cultural landscapes over the time: highlighting the marine and mountainous environment of Mirabello" with MIS 5028190. DIATOPO is part of the Priority Axis "Strengthening the competitiveness, innovation and entrepreneurship of Crete" of the Operational Program "Crete" 2014-2020 and co-financed by the European Regional Development Fund.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The spatial raw data are freely available to the academic and research community from the author upon reasonable request.