Tested in the Gulf of California Suggests Field Data are Not Needed to Derive Satellite Bathymetry

Satellite derived bathymetry methods over coastal areas were born to deliver basic and useful information like bathymetry. However, the process is not straightforward, the main limitation being the need of field data. The Self-calibrated Spectral Supervised Shallow-water Modeler (4SM) method was tested to obtain coastal bathymetry without the use of any field data. Using LANDSAT-8 multispectral images from 2013 to 2016, a bathymetric time series was produced. Groundtruthed depths and an alternative method, Stumpf’s Band Ratio Algorithm, were used to verify the results. Retrieved (4SM) vs groundtruthed depths scored an average r2 (0.90), and a low error (RMSE = 1.47 m). Also 4SM showed, over the whole time series, the same average accuracy of the control method (40%). Advantages, limitations and operability under complex atmosphere and water column conditions, and high and low-albedo bottom processing capabilities of 4SM are discussed. In conclusion, the findings suggest that 4SM is equally accurate as the commonly used Stumpf’s method, the only difference being the independence of 4SM to previous field data, and the potential to deliver bottom spectral characteristics for further modelling. 4SM represents a significative advance in coastal remote sensing potential to obtain bathymetry and optical properties of the marine bottom.


Introduction
Accurate bathymetry in the coastal area is one of the most basic yet fundamental sources of information for many purposes like navigation, military, resource exploitation, fisheries and tourism.Water depth is also a major factor in community zonation and it is therefore relevant to map and model benthic habitats.In marine sciences, acoustic, radar and optical systems have been used to assess bathymetry [1].However, they require substantial time and effort to be deployed, operated and interpreted.Such an effort can be challenging, especially in developing countries.Furthermore, vast regions are too costly or difficult to reach for these conventional methods: this creates gaps in the knowledge of underwater habitats.Recently, satellite technologies have developed rapidly and the NASA-USGS Landsat Data Continuity Mission [2] offers a unique opportunity to get open-source multispectral images (https://earthexplorer.usgs.gov).LANDSAT-8 is the newest multispectral satellite developed from the Landsat series which have been used to identify benthic objects in shallow waters [3][4][5][6][7][8][9].LANDSAT-8 features a higher spectral resolution and quality than its predecessors.The availability of images makes LANDSAT-8 the main open source alternative for coastal monitoring.However, the challenge is to obtain the underwater optical properties which are required to model bathymetry and bottom reflectance.
Among the several factors that determine the reflectance spectrum recorded by the satellite, there are four main issues which the practitioner has to work around to obtain acceptable results: 1) the influence of the atmosphere through Rayleigh scattering; 2) the sun-glint effect where photons are reflected by the air-sea interface; 3) the interaction between refracted photons, water molecules and other constituents in the water column; and, finally 4) the interaction of photons with the bottom [10].The way to account for these influences is by applying "corrections" to the images.The water column correction algorithms published to date have been extensively reviewed [10,11].They differ in the way of estimating the several contributions to the water leaving spectral signature, and all of them require field data for calibration purposes.They are based on the knowledge of the physics involved (analytical) or based on field data (empirical), or both (semi-analytical).The most popular among empirical methods are band combination algorithms, which are simple and easy to use; they assume that the effects of the bottom albedo in multispectral bands are exponential functions of depth.The water leaving signal over a shallow bottom is a function of the optical properties of the water column, like the reflectance over deep water (Lw) and the diffuse attenuation coefficient (Kd).These are essential parameters for water column correction methods [10], and have also been used to classify water types [12], estimate chlorophyll concentrations [13], predict light intensity at depth [14], estimate the rate of solar heating in surface waters [15].Since the depth at each pixel is constant for each multispectral band of the satellite image, the relationship between radiance in two distinct bands can be used to estimate Kd.Estimating Kd from field data is difficult; it can be approximated using various band ratio methods [16][17][18][19].The first and most popular [16], assumes that differences of a group of pixels on the same substrate are due to depth variation and not another factor.This can be true in small areas, but if the area is large, or the bottom type is heterogeneous, the only solution is a calibration with ground-truthed Kd measurements [20].Moreover, Lyzenga's algorithm does not retrieve substrate reflectance, and the performance depends on the wavelength of the bands applied.While the algorithm removes the depth effect, it generates a depth-invariant bottom index which is not related to the physics behind [10].Lyzenga et al. [21] further enhanced the methodology using multiple linear regressions which, however, still required knowledge of water's optical properties at the time of image acquisition, which is not always available nor easy to obtain.Stumpf et al. [22] presented a Log Ratio method to account for variable bottom types, which has been used since then by NOAA for coastal mapping and is available in an image processing software (ENVI, and SPEAR toolbox).This method still requires existing depth measurements for calibration and does not produce a water column corrected image.The method developed by Sagawa et al. [19] for low transparency waters requires depth and attenuation coefficient observations measured in the field.It may be valid only within small areas and authors emphasize the need for an accurate sea-truthed map to obtain results.Therefore, it does not produce an estimation of retrieved bottom reflectance.
In contrast with all other methods in use to-date, the Self-calibrated Spectral Supervised Shallow water Modeler (4SM) [23], does not rely on preliminary sets of field data for the calibration of the model parameters.4SM was first presented by Morel and Lindell [24] and is explained in detail by Morel and Favoretto [23] and at http://www.watercolumncorrection.com/.We tested the 4SM approach to process a time series of nine LANDSAT-8 images without the use of any field data for the calibration of the simplified radiative transfer equation (RTE).Furthermore, we compare the results with those obtained using the Log Ratio method and with a depth sounding dataset.The potential of 4SM to retrieve depths from multispectral images without the need for field data is a significant advance in the exploration and monitoring potentials in coastal areas.

Study area
The San Lorenzo Channel, is a shallow area that connects the Bay of La Paz to the Gulf of California (Figure 1).The communication between the Bay and the Gulf of California happens also through the "Boca Norte" (Northern mouth): this large and deep opening reaches 350 m deep.San Lorenzo Channel was chosen for its morphology and for the variety of shallow habitats [25,26].The region has two contrasting seasons: a dry season (October-June) and a humid season (July-September); the latter is referred to as "Mexican monsoon" [27] and is associated with occasional significant precipitation, often through intense events like hurricanes.Hydrographic conditions in the bay, particularly in deeper areas, are strongly influenced by the seasonal variability of main thermohaline circulation of the southern Gulf of California [28].Salinity ranges between 34 and 36 depending on the time of the year [29].In particular, the deeper layers, are more influenced by the Gulf of California, while shallower ones are dominated by local processes like insolation, evaporation, precipitation and mixing [30].Because of the high solar radiation during summer, and weak wind regimes, in the upper 20 m an intense thermocline can form, and the intense evaporation cause the surface layer to be more saline.This situation can be reverted with strong monsoonal winds, or during the end of the humid season characterized by northerly colder winds that promotes vertical mixing of the water column.In particular, San Lorenzo channel can show a vertical and horizontal variability of the water column caused by a combinations of events that generate a highly dynamical habitat.Among them, upwelling of colder waters from the Gulf of California can cause localized phytoplanktonic blooms [31, and references therein], and strong tidal currents [32] can transport sediments from the bottom and enhance turbidity of the bottom layers of the water column (Author, FF, pers.obs.); finally, during the monsoon season, sediment runoff increases due to intense episodic precipitations, decreasing water quality, while atmospheric conditions become more complex (see Appendix A).

Radiative transfer equation
4SM relies on the simplified radiative transfer equation (RTE) for optically shallow water (equation 1).The RTE has been discussed in detail by Philpot [33] and Maritorena et al. [34].In this algorithm, the irradiance reflectance of shallow waters below the surface is equal to the deep water reflectance, plus bottom substrate contrast attenuated by the depth effect, using the exponential decrease term e (2KZ) in equation 1, to give the inverted RTE in equation 2. Following Maritorena et al. [34] for nadir radiances at the water-air interface (Bottom Of Atmosphere, BOA) the RTE is: it is inverted as BOA: or at the sensor TOA: = + ( − ) ( ) .
(3) Where L is the spectral BOA radiance where the bottom substrate is at depth Z; Lw is the spectral BOA radiance of the sea where the bottom is optically deep, also defined as backscatter or water volume reflectance; LB is the spectral BOA radiance of the bottom substrate at null depth.By adding the atmospheric path radiance (La) to all BOA radiance terms, one gets the Top Of the Atmosphere (TOA) spectral radiances at the sensor LsB, Ls and Lsw in equation (3).2K is an operational two-way diffuse attenuation coefficient, in units of m-1, equals to the down-welling and up-welling term together (2K = Kd + Ku) and Z is the actual depth of the seabed.The optical depth 2KZ, also called optical path length, is dimensionless.Wavelength dependency subscripts of all L and K terms are omitted for brevity.This will finally allow for the linearization of the bottom contrast: X=log(Ls-Lsw).
(4) 4SM is a "ratio method".For a pair of bands i and j, with 2Ki<2Kj and suitable bottom contrast in band j, water column correction is achieved by increasing Z in (2) until the ratio LBi/LBj matches the Soil Line (see the optical calibration section below); LBi is the average of all bands with 2Ki<2Kj and suitable bottom contrast.This process operates ratios among relative numbers, and therefore does not require radiance terms to be converted into calibrated reflectance.All L radiance terms in the above equations represent the signal as measured by the sensor, and the terms digital numbers (DN), radiance and reflectance are used interchangeably and remain dimensionless.This is the privilege of a "ratio method": the 4SM approach does not require formal atmospheric correction.

Image resources preparation
It was soon realized that the hydrological conditions of the La Paz area are complex and highly variable.Therefore, a time-series of LANDSAT-8 images have been downloaded from https://earthexplorer.usgs.gov/;all scenes were searched with a less than 10% cloud cover condition.Preference was given to images from different months and seasons, avoiding duplicates from the same month, since such a time-detailed study was not the aim of the work.Nine images were selected to be processed and considered reasonably safe to be modeled (Table 1).We processed them using the 4SM method, then using the ENVI/SPEAR tool method.As regards ENVI/SPEARS, the image preparation and preprocessing started with the creation of vector polygons to be used as masks for clouds and boats.Then, if necessary, the images are deglinted following Hedley et al. [35].The de-glinting assumes a null penetration in water of the NIR-SWIR bands of the spectral images.Areas that are affected by glint are selected over deep water and then are used to calculate linear models of all bands against the NIR or SWIR band.The amount of glint is evidenced and then subtracted using the linear models.The images are then pansharpened using a modified version of the OrfeoToolbox algorithm in R-project in the Raster Package [36].In this process, the generation of a high resolution image (PXS) is achieved through the fusion of the panchromatic (PAN) with the multispectral (XS) data.First, the two images need to be orthorectified by oversampling of the XS image, then a low pass filter is applied to the PAN band to equalize (in the Fourier domain) the spectral content with XS data.Once the spectral content is similar, XS are normalized with the low-pass filtered (PAN) and results are multiplied with the original panchromatic band.The whole process can be represented as: As regards 4SM, the image preparation only consists in importing the raw data as digital numbers in a working database, and creating a few shapefiles for sampling glint, optically deep water, vegetation, and bright sandy beach.Deglinting and pansharpening in 4SM uses the same protocols as above; however, they are applied on the fly rather than as image pre-processing.

Optical calibration of the simplified RTE in 4SM
As presented by Morel and Favoretto [23] and at www.watercolumncorrection.com, it is possible to retrieve from the image all the information needed for the calibration of the simplified RTE.Solving for Z and LB in equation (2) requires four spectral parameters which are estimated from the image itself: 1) the deep water reflectance Lsw; 2) the Soil Line; 3) the water volume reflectance Lw; 4) the Brightest Pixels Line.
Estimating the TOA deep water reflectance Lsw requires the presence of an area in the image which is free of any glint (or de-glinted), and is optically deep in all wavebands.
To estimate the Soil Line [37], the non-vegetated areas on the land part of the image are analyzed in order to provide the operational values of the ratios Lbi/LBj of the bottom substrate at null depth for all possible pairs of wavebands i and j.These ratios are then assumed to apply to spectrally nondifferentiated shallow water pixels.
Estimating the water volume reflectance Lw is done by visual inspection of the Soil Line in the optical calibration diagram, and using some commonsense knowledge.Therefore, the atmospheric path radiance La may be estimated as the difference La=Lsw-Lw.Please note that, by subtracting La from all Ls terms in equation ( 3), a first order atmospheric correction of the imagery is achieved.
The Brightest Pixel Line is a radiometric model of the brightest shallow bottom over the whole deep range in the image.It assumes that the shallow water part of the image contains at least some stretches of the brightest bottom substrate that exists at the scene at various depths; it represents the exponential decrease of the radiances as the bottom depth increases.The BPL is used to estimate the ratio Kblue/Kgreen of the attenuation coefficients in the blue and green bands.Then this ratio is used to interpolate spectral K for all visible bands, using table XXVII of Jerlov [12]; this dataset reports the diffuse attenuation coefficients for downwelling irradiance in oceanic and coastal waters worldwide.
For two examples of straightforward optical calibration diagrams over apparently homogeneous waters (one Oceanic water type, and one Coastal water type), please refer to diagrams in Appendix B. For most scenes, the calibration diagrams exhibit a complex hydrologic situation, which we tentatively explain as follows: the clear deep oceanic waters are usually covered by a layer of slightly turbid coastal waters; the thickness of this upper layer varies across the scene, from just a very few meters to in excess of 30 m (see Appendix B).This vertical structure is illustrated in Figure 2 for a Landsat 8 scene acquired on January 1rst 2015: X2, X3 and X5 are linearized Blue, Green and Red bands respectively (equation 4).The backdrop in shades of gray represents the bi-dimensional histogram for all image pixels (white ROI in Figure 1).Blue dots are the actual pixels used for BPL calibration, each referenced to its row/line position in the image.Small white circles represent the scatter of averaged bare land pixels, as a proxy for the Soil Line.This diagram exhibits two layered water types: • 0-5 m depth range: on that day, the BPL pixels in Figure 2a  Finally, through a series of interactive calibration steps, the practitioner now must progressively fine tune manually all these variables in the command line script, until satisfied that a sound consistency has been achieved: the calibration diagrams allow to take control of the calibration process.This completes the set of parameters that are necessary for operating the inverted RTE (equation 2) without the need or use of any field data.

Depth retrieval in 4SM
4SM is a ratio method, meaning that a band ratio is used to derive both Z and a spectral bottom signature.Let Lbj be the water column corrected radiance in band j with suitable bottom contrast Ls-Lsw; Let Lbi be the average water column corrected radiance of all available wavebands with Ki<Kj and suitable bottom contrast.The inversion of the model is achieved by increasing Z term in equation 2 until the ratio LBi/LBj matches that of the Soil Line.As detailed by Morel and Favoretto [23], for a Landsat-8 or WorldView-2 image, several solutions are available: the NIR solution, the Red solution, the Green solution, and the PAN solution.In practice, the PAN solution is preferred, as it carries several distinct advantages.The outputs are (i) a raster DTM of the shallow water area in units of meter, and (ii) several rasters of water column corrected wavebands in units of relative radiance, ready for bottom typing (a "low-tide" view of the scene).Some seatruth data can be used for fine tuning the estimation of Z and for a tide correction.

Combining depths in 4SM
Single scenes are inevitably affected by alien features that impact the result (e.g.boats and their wakes, small clouds and their shadows, algal blooms, local variations of water optical properties, etc.).Therefore, a combined depth image can be produced by averaging all the retrieved depths that satisfy some quality criteria: for each shallow pixel, first an average depth is calculated; then depths that are outside the standard deviation range are excluded and a new average is calculated.This eliminates discrete artifacts without any smoothing.See for example in Figure 3 that, on November 4th 2013 (purple profile), a major depth under-estimation affects the southern half of the deep San Lorenzo Channel between 4 and 5 km on Profile A; an increase in turbidity possibly caused this artifact.These profiles also show that the depth combination reduces much of the noise without any smoothing.

Groundtruthed DTM
The collection of ground truth depth soundings in the San Lorenzo Channel was completed in two days by the use of a boat equipped with a Lowrance fishfinder Elite5x HDI.The Transducer (Airmar TM260 1kW with 50kHz wide 19o beam / 200kHz narrow 6o beam), was set to read from the waterline.The offset between the antenna and the transducer position was automatically calibrated by the instrument, and the instrument flagged data output in case of excessive wave action with an internal gyroscope.However, during the whole groundtruthing process, sea conditions were excellent.Data have been collected at a constant speed of 10 kn on parallel transects.Up to 8 depth readings were averaged for each UTM navigation point, then tide height was subtracted using open source tide data (http://predmar.cicese.mx/calendarios/).Finally, depth data were converted in a shapefile (n= 5004, red points in Figure 1), and a Digital Terrain Model (DTM) was created by spline interpolation of groundtruthed depths (Figure 4A).

Depth retrieval in ENVI-SPEAR tool
Eight scenes were processed with the ENVI-SPEAR tool over the same Region Of Interest (ROI) as 4SM (Figure 1).This tool operates the Log Ratio Transform algorithm developed by Stumpf et al., [22].It is integrated in the ENVI software (Exelis Visual Information Solutions, Boulder, Colorado, v.5.3) as a Spectral Processing Exploitation and Analysis Resource (SPEAR) tool, from now ENVI-SPEAR tool is capable of retrieving relative depths from multispectral images which then are calibrated into meters using existing depth sounding dataset.To this end, the groundtruthed depths described above were used.Eight images were processed.An example of the output is shown in Figure 4C.Alien features were masked.Then an approximate atmospheric correction was applied on the whole image by Dark Pixel Subtraction.Images were pansharpened, then the Log Ratio Transform algorithm was applied.A median filter 3x3 was chosen to remove high frequency noise.Then all groundtruthed depth points (n = 5004) were used to calibrate relative depths into absolute depths.In the calibration tool (Figure 5), it is possible to select a model that best fits data based on r2 results (Linear, Exponential, Polynomial model).To compare different results, a linear fitting (Figure 5A) and the model with the highest r2 were chosen for each scene (e.g. Figure 5B).

Groundtruthing regressions and comparisons
As regards the accuracy of 4SM depths, pixels with both groundtruthed depths and satellite derived depths are selected and these depth pairs are assigned to bins in decimeters [23].Bins with few data pairs are excluded as outliers and plotted as red points in the scatterplots of Figure 6: these outliers are generated by various alien features and must be excluded for a fair regression calculation.Then all the remaining bins are averaged into bins in meters.In Figure 6, all one-meter bins which are represented by a star symbol account for a total of 99 % of all accepted depth pairs, while the remaining one-meter bins are represented by a white circle.Inside the two parallel gray lines in the scatterplot, satellite retrieved depths within ± 1 m of groundtruthed depths are included.A correlation coefficient (r2) and a Root Mean Squared Error (RMSE) were calculated for each linear model: see Figure 6 and Table 1.Tide heights were adjusted manually and added to seatruth depths in order to minimize the RMSE result.Two scenes stand out with unrealistic tide heights (1.6 and 1.7 m): this is possibly the result of a bad estimation of the Soil Line for these two scenes.In order to compare 4SM and ENVI-SPEAR tool results, an accuracy index was estimated by calculating Depth_residual = (Depth_retrieved -Depth_measured).So a positive depth residual signals an overestimated depth, and vice versa.Depth residuals have been classified.Figure 7 shows 5 classes of accuracy: blue tones for two depth overestimation classes, red tones for two depth underestimation classes, and turquoise tone where the residual is less than 1 meter.A barplots in Figure 7 shows the percentage of pixels that belong to each of these classes.

Results
This section first presents the results of the seatruth regressions (Figures 5-6), then the results of the accuracy assessment for retrieved depth (Figure 7), and finally the depth residuals (Figure 8).

Seatruth regressions
As regards 4SM results, linear relationships between retrieved depths and groundtruthed depths are shown for each scene in Figure 6.The regression lines do not account for all depth pairs; red outliers are excluded from the statistics.Linear regression model explained in excess of 80% of the depth variations.The October 27th 2016 scene reports the highest r2 (0.99) although it also reports the highest RMSE (2.11 m).The October 22th 2014 scene scored the lowest RMSE.From Figure 6 it is also possible to see how combining depths improved linear fitting results and reduced the noise.
For ENVI-SPEAR tool results, please refer to Figure 5.In Figure 5A a linear model fitting is applied to the scatterplot, while in Figure 5B the best model that showed the higher r2 is plotted.In particular, for eight ENVI-SPEAR scenes regression plots, the higher correlation coefficient r2 was always achieved with a polynomial model.

Accuracy assessment
Accuracy assessments are presented in Figure 7. A) is a LANDSAT RGB color composite view of October 22th 2014, white arrow indicates an example of bright shallow bottom, while black arrow represents darker shallow bottoms.A.1) is 4SM accuracy index; A.2) is the accuracy index calculated on ENVI-SPEAR with line model calibration; A.3) is the accuracy index calculated on ENVI-SPEAR with polynomial model calibration.Topright barplot report % of pixels belongings to each index classes, <-5 m differences class is due to the detection limit constraint of remote sensing in coastal water, since it is represented by pixels in deeper waters.B) is a LANDSAT RGB color composite view of October 19th 2013, B.1) is 4SM accuracy index; B.2) is the accuracy index calculated on ENVI-SPEAR with line model calibration; A.3) is the accuracy index calculated on ENVI-SPEAR with polynomial model calibration.Topright barplot report % of pixels belongings to each index classes calculated on the scenes.The accuracy index allowed a visual estimation of overall accuracy of the retrieving methods applied for each scene.
In terms of % of accuracy between ±1 m from groundtruthed depths, 4SM performed better in 4 out of 9 scenes; in the 3 remaining scenes, it scored similar results with ENVI-SPEAR with the exception of one scene: October 22th 2014 (Figure 7 A.1 and black bars in linked barplot chart) that is reported to be further discussed.In Figure 7 B a scene from October 19th 2013 is reported because 4SM scored its best results (50 % of pixel between ±1 m from groundtruthed depths, Figure 7 B.1 and black bars in linked barplot chart).Over shallow bright bottoms, 4SM performed better than ENVI-SPEAR (Figure 7 A, B white arrows); in particular, with a linear fitting model, ENVI-SPEAR tends to underestimate depth by more than 3 m over bright bottom areas (Figure 7 A.2, B.2).It is also observed, in Figure 5A, that linear fitting of ENVI-SPEAR exhibits a detection limit: deeper than ~15 m, most retrieved depths are not correlated to real depths anymore; polynomial model exhibits a better fit because of the curved shape of scatterplots above 0.8 Log Ratio results (Figure 5B).Moreover, it is observed, in Figure 7 A.1, A.2, A.3, that both 4SM and ENVI-SPEAR tend to overestimate depths over dark bottoms at shallow depths.Finally, except for October 22th 2014 and October 27th 2016 which both exhibit complex atmospheric and hydrologic conditions (Appendix A4), accuracy of the methods is similar.Therefore, in spite of the adverse conditions of the San Lorenzo Channel, and without the need of any field data, the 4SM method is comparable to the ENVI-SPEAR method, and equally valid.

Depth residuals
Depth residuals are defined by the absolute difference between all of 4SM retrieved depths from all the nine scenes investigated.Depth residuals are represented graphically in Figure 8.Most of the depth residuals were ≤ 0.5 m, especially over bright shallow and deeper (≈ 15 m) substrate, where 4SM showed more consistency.Conversely, higher residuals were found in two areas corresponding to darker bottoms (see also Figure 7 A, B black arrows).Over time this area showed a change in bottom properties, when in October 22th 2014 a large portion of a low albedo patch on the northern side disappeared and wasn't recorded anymore.Finally, it is evident how areas below 20 m were the ones that showed higher residuals, where depth retrieval reaches its optical limit.

San Lorenzo Channel's conditions
The results suggest that 4SM can be safely used on single scenes under reasonable conditions of low cloud cover, low glint effects, and homogeneous vertical and horizontal water optical properties.The San Lorenzo Channel did not show any relevant variations in depth along the 4 years analyzed.However, the area was characterized by high seasonal variations in vertical and horizontal stratification of the waters; the same can be said for variations of spectral characteristic of the bottom, that are currently under investigation (Favoretto et al., in prep.).The application of 4SM not only allowed the retrieval of the bathymetry of the San Lorenzo Channel that showed to be consistent with groundtruthed depths, but it generated a bathymetry model for the whole shallow area in the Bay of La Paz (Figure 9).In particular, the San Lorenzo Channel showed a fair amount of layered waters with seasonal and spatial complexity, and 4SM does account to a limited extent for such layered waters while calibrating the simplified RTE.Overall, waters in the channel were vertically stratified, with common algal blooms events (red tides) in the complex hydrodynamic conditions in the Bay of La Paz ( [31] and references therein).In the October 22th 2013 scene for example, deep oceanic waters are covered by an upper layer of Coastal waters: in the San Gabriel bay (Figure 1, red square), this upper layer is ~5 m thick, but it is ~9 m thick in the San Lorenzo Channel.Such heterogeneity is commonly found through the seasonal cycle the scenes investigated.All of this considered, it is possible to obtain useful information even from low quality areas in the scenes, accounting for their limitations.

Models comparison
In Figure 7 A, B, a RGB color composites present the worst (A) and the best (B) scenes in terms of 4SM accuracy.It is important to note a major change in bottom composition from 2013 that clearly show a dark patch (black arrow figure 7A) of what probably was macroalgae (Favoretto et al., in prep.) that disappeared in 2014, and up to now it has not been detected again.This change may have been caused by Hurricane Odile in September 2014, that displaced significant amounts of sediments (Favoretto et al., in prep.).This event allowed for a comparison between low-albedo bottom processing capabilities of 4SM and ENVI-SPEAR.The bottom albedo independent nature of ENVI-SPEAR (Stumpf's) algorithm should guarantee in some limits that dark bottoms are shown at the same depth as bright sand, when field data report so.However, the method applied in the San Lorenzo Channel shows how depths over dark bottoms are overestimated within 5 -10 m depth range (Figure 7 B.2, B.3).In particular, the linear model calibration of ENVI-SPEAR algorithm yields higher depth residuals than the polynomial model calibration.(Figure 7 A.2, A.3).This limitation that was already pointed out in the introduction of this article about Stumpf's algorithm.As the results suggests, 4SM perform better over very bright shallow areas, but has the same problems as ENVI-SPEAR over very dark bottoms (Figure 7  In the 3-10 m range, retrieved depths are over-estimated by 4SM: this means that bottom substrates are actually less green than assumed by default in 4SM.In this depth range, dense red coralline algae formations, called rhodolith beds, characterize the area [25,26].Coralline algae absorb additional wavelength over the green range for their photosynthetic activities, and their common presence in the area could in part explain a lower than expected green reflectance.This hypothesis is currently under investigation (Favoretto et al., in prep).
Over the 10-20 m range, most retrieved depths are under-estimated by 4SM: this means that bottom substrates are actually greener than assumed by default.This depth range is present in the San Lorenzo Channel exclusively in its central portion, where tidal currents are supposed to merge in stronger bidirectional deeper currents.In this portion of the channel, where strong currents shape the sandy bottom: cyanobacterial biofilms can cover large portions of the seabed (Favoretto et al., in press).These organisms, stretching for meters, can have a significant production and affect the green reflectivity of the substrate in the green band (see Appendix D).

The 4SM method
In 4SM the detection of the shallow bottom in more than two wavebands allows for a better estimation of Z, and therefore yields a richer spectral bottom reflectance associated to the pixel.To this extent, and from many tests made in the development of 4SM (see www.watercolumncorrection.com , and [23]) it can be said that 4SM is not site-specific, and is capable of delivering both a bathymetric map and a water column corrected image in units of image's DN over many different water types, in different location of the world, when the high cost of collecting similar data in the field is simply not affordable.While an empirical method (e.g.[21]) explores the statistical relationship between image pixel values and field measured depths, 4SM operates a radiance inversion approach which does not need any field measured water depths for calibration.Using 4SM there is no need for a formal atmospheric correction since the spectral path radiance and the water volume reflectance are estimated through the calibration process.All of these advantages reduce artifacts and enhance increase operability simplifying image processing and interpretation.
The RTE results are an estimate of water depth, and an acceptable fit of the water column corrected spectral bottom signature of each shallow pixel with the Soil Line, although contrasted bottom signatures entail potentially severe bias on retrieved depth.Another advantage of the 4SM method is the pansharpening adaptation, that generates 15 m resolution images to be modelled.Using this modified pansharpening algorithm, pixels which are brighter than their surrounding in the PAN band get a boost in the multispectral band and vice versa.Other empirical methods use only multispectral bands, and this is applicable as long as there is enough contrast in colors (e.g.Kblue<Kgreen<Kred) like in clear oceanic waters.However, in coastal waters, this is not the case and there is hardly any color separation between Kblue and Kgreen.This actually is a fundamental limitation which precludes a reliable estimation of the model parameters for Coastal water types in methods like Lyzenga's and Stumpf's.The PAN solution offers a valuable alternative, as it will always generate enough color contrast, even with Coastal type waters.In the end, 4SM uses the image metadata to convert spectral water column corrected bottom signatures into calibrated reflectance (0-1), ready for bottom typing and time series studies.The method produces a low tide view of the shallow areas along with a DTM; this is a unique feature of 4SM, apart from semi-analytical methods [10].
Moreover, the application of combined depths increased model fitting results; however, it did not improve accuracy compared to other scenes.Strong depth residuals (>2m) can hardly be interpreted in this case as caused by natural or anthropic events; they are likely noise in depth retrievals due to local perturbations like atmospherics or hydrological effects, but also phytoplankton blooms, turbid waters, wind or even boats wakes.Therefore, any pixel that shows this kind of variation should not be considered in the overall depth model.Combining depths is an important tool in reducing these types of noise (that are inevitable for the practitioner).The fact that most of the pixels showed little or no variation (0 -0.5 m, Figure 8) in the surface to 20 m depth range is supporting 4SM as a robust depth retrieval method that can cope, to a limited extent, with a layered structure of the water column.

Conclusions
It is possible to conclude about 4SM that: • is independent on field data to achieve the optical calibration; • offers a simpler operational method compared to previous methods (e.g.[16,19,22]); • its accuracy is equally valid, and in some case better, compared to the accepted and most used Stumpf's algorithm; • provides a priori important insight on ROI water column and hydrological condition; • use all bands with significant bottom detection and delivers computed depths, spectral K (i.e. water quality) and spectral water column corrected bottom reflectance; • provides all valuable information that allows to explore and monitor large coastal areas in a more efficient and cheaper way in terms of resources and time; These results highlight the progress made by 4SM as a band ratio method for satellite derived bathymetry, especially if considered the independence from any field calibration, and show how the 4SM algorithm can be a relevant alternative to traditional ratio-methods.

Figure 1 .
Landsat scene: LC80340432013292LGN00.(a) whole image with study area (orange rectangle), and Espiritu Santo island and the San Lorenzo Channel (white rectangle); (b) red square is the area used for parameter calibration, and green rectangle is the San Lorenzo Channel clearly display along a straight line over the 0-5 m depth range.The ratio Kblue/Kgreen for this surficial layer is estimated at 0.95; this denotes a water type C1+0.17 of Jerlov.The BPL pixels in Figure2bdisplay along two straight lines for the pairs Blue/Red and Green/Red.These two straight lines have virtually the same slope.Attenuation coefficients in units of m-1 are estimated at Kblue=0.272, Kgreen=0.285,and Kred=0.774.Please note that 0.272/0.285=0.95;• 5-10 m depth range: Figure 2a seems to exhibit a progressive change in water quality over the 5-10 m depth range; • 10~25 m depth range: on that day, the the BPL pixels in Figure 2a clearly display along a straight line over the 10-~25 m depth range.The ratio Kblue/Kgreen for this deeper layer is estimated at 0.75; this denotes a water type OII+0.53 of Jerlov.Attenuation coefficients in units of m-1 are estimated at Kblue=0.173, Kgreen=0.232.

Figure 2 .
Figure 2. Example of linearized calibration diagram for the white ROI of Figure 1

Figure 3 .
Figure 3. Plot of the depth combination process along four profiles across San Lorenzo Channel.All profiles run from North to South.
a) is the interpolated DTM of groundtruthed depth points; b) is 4SM satellite retrieved depth layer from October 19th 2013; c) is satellite retrieved depth layer processed with ENVI-SPEAR tool on the same scene.

Figure 6 .
Figure 6.Scatterplots of 4SM retrieved depths vs tide corrected groundtruthed depths.Please note the different scale in the combined depths plot.

Figure 7 .
Figure 7. Results of the accuracy assessment on retrieved depth for 4SM and ENVI/SPEAR tool methods.Please refer to text for details.

Figure 8 .
Figure 8. absolute difference between retrieved depths in all of the 4SM retrieved depths scenes

Figure 9 .
Figure 9. 4SM is capable of calculating bathymetry values without field data and this is the results of satellite derived bathymetry for the whole La Paz bay.

Figure A2 .
Figure A2.True color composite of the October 27th 2016 image.This scene exhibits a complex atmosphere and distinct blooms of discolored waters

Figure A3 .
Figure A3.True color composite of the October 27th 2016 image.This scene exhibits distinct blooms of discolored waters (algal blooms).

Table 1 .
Landsat 8OLI scenes codes processed in this study, with seatruth regression for 4SM retrieved depths.