Shallow Bathymetry from Multiple Sentinel 2 Images via the Joint Estimation of Wave Celerity and Wavelength

: In this study, we present a new method called BathySent to retrieve shallow bathymetry from space that is based on the joint measurement of ocean wave celerity (c) and wavelength ( λ ). We developed the method to work with Sentinel 2 data, exploiting the time lag between two Sentinel 2 spectral bands, acquired quasi-simultaneously, from a single satellite dataset. Our method was based on the linear dispersion law, which related water depth to wave celerity and wavelength: when the water depth was less than about half the dominant wavelength, the wave celerity and wavelength decreased due to decreasing water depth (h) as the waves propagated towards the coast. Instead of using a best weighted ( c, λ ) ﬁt with the linear dispersion relation to retrieve h, we proposed solving the linear dispersion relation for each ( c, λ ) pair to ﬁnd multiple h-values within the same resolution cell. Then, we calculated the weighted averaged h-value for each resolution cell. To improve the precision of the ﬁnal bathymetric map, we stacked the bathymetry values from N -different datasets acquired from the same study area on different dates. We ﬁrst tested the algorithm on a set of images representing simulated ocean waves, then we applied it to a real set of Sentinel 2 data obtained of our study area, G â vres peninsula (France, 47 ◦ ,67 lat.; − 3 ◦ 35 lon.). A comparison with in situ bathymetry yielded good results from the synthetic images (r 2 = 0.9) and promising results with the Sentinel 2 images (r 2 = 0.7) in the 0–16 m depth zone.


Introduction
The retrieval of bathymetry from space is of major importance, particularly in sea areas with shallow depths, where access by conventional methods is tricky, and in areas where shallow morphology changes with time. Many methods based on multispectral imagery exist, which are appropriate in coastal areas where the sea floor can be directly imaged (e.g., [1][2][3][4][5][6] and references therein). In this study, we concentrated only on methods based on ocean wave motion from spaceborne optical imagery. In particular, we concentrated on measurements of wave celerity (c) and wavelength (λ) from quasi-simultaneous image acquisitions. These methods were promising in areas where the water was non-transparent, and they did not require ground calibration.
the dominant wavelengths and their associated celerities using the linear dispersion relation (also referred to simply as the dispersion relation) for free surface waves (e.g., [8]) [1] as: where h is the water depth and g is the universal gravitational constant. Equation (1) was first used in a practical sense during the Second World War [9], and works either with two sequential images acquired within a short time interval (i.e., shorter than the wave period), or with one image but with additional in situ measurements (e.g., [10]).

Related Works
In this study, we concentrate on spaceborne methods to measure wave celerity and wavelengths from space to retrieve bathymetry. If the reader is interested in other methods, using selective colors absorption for instance, we suggest to read Salameh et al., [11] and the bibliographytherein.
The retrieval of wave celerity from satellites is still challenging today because we need at least two consecutive images acquired with only a few s time delay. Such images can be obtained onboard specific satellites with push-broom sensors. With these sensors, the retrieval of c relies on the non-simultaneous stereo acquisition capability of given platforms, as shown by R. Abileah [12]. Or the information about c can be obtained when two multispectral bands cannot simultaneously exist in the focal plane of the sensor, yielding a de facto short time lag between two bands in the same dataset, allowing one to measure c directly, as proposed by de Michele et al. [13]. With these sensors and the joint measurement of c and λ, h becomes accessible.
In the literature, a few studies have proposed ways to access c and λ from space. For example, Abileah [12] suggested a method to measure c from optical satellite data using Ikonos stereo pairs that had an 11 s time lag. de Michele et al. [13] proposed a method that exploits the small time-lag (2.1 s) between two bands of SPOT-5 data. Danilo and Binet [14] have detailed the concept and limiting factors to access c and λ from space. Based on the aforementioned studies, Poupardin et al. [15,16] proposed a local spectral decomposition method that was employed prior to making celerity estimates to produce local (c, λ) pairs that were used to derive the water depth based on the dispersion relation (Equation (1)). Bergsma et al. [17] proposed a method based on multiple Sentinel 2 bands and the Fourier slice theorem. Almar et al. [18] proposed a method to extract bathymetry from the Pleiades satellite's "persistent" mode, based on (c, λ) measures.
Abileah [12], Bergsma et al. [17] and Almar et al. [18] make use of only one (c, λ) couple, the most energetic one, to extract one bathymetric value from a resolution cell. In this study, we furthered the work of de Michele et al. [13] and Poupardin et al. [15,16] to retrieve water depth by measuring multiple (c, λ) pairs in each resolution cell and we apply it to Sentinel 2 data. The advantages of using Sentinel 2 were manifold. First, the enhanced revisit time (five days) allowed one to improve the data acquisition frequency. Second, the data covered a wide swath, which could be potentially used to map large areas. Third, the data were provided free of charge.

Our Approach
Our method exploited multiple (c, λ) pairs within the same resolution cell. For the purpose of fast processing, we employed two steps: FFT (fast Fourier transform) and the weighted average of multiple h-values within the same resolution cell. First, we used FFT-based methods to produce local, sub-sampled, wave spectra and their associated celerities. To calculate the celerities, we exploited the time lag between two Sentinel-2 bands, acquired during a single Sentinel-2 pass. There are 1.05 s between acquisitions made in band 2 and band 4 (e.g., [19]). Band 2 records sunlight reflected from the Earth surface at the central wavelength of 0.49 µm (perceived as blue by the human eye) and band 4 at 0.68 µm (perceived as red by the human eye). First, we calculate the FFT spectra of these two images. Then, phases and spatial frequencies (i.e., wavenumbers) were combined (Section 2.2) to derive the multiple (c, λ) pairs.
The wave spectra were calculated from small sub-images on the basis of a grid defined by the sampling step, where the sampling step defined the resolution of the final bathymetric map. We calculated multiple (c, λ) pairs per resolution cell, but instead of using a best weighted (c, λ) fit with the linear dispersion relation to retrieve h, we proposed solving the linear dispersion relation for each (c, λ) couple. In this way, we found multiple h-values within the same resolution cell. Then, we calculated the weighted average h-value for each resolution cell per dataset, across N datasets. Then, we stacked the results from the N datasets. Because this method was conceived for Sentinel 2, we called it BathySent.
In this study, we tested our procedure on simulated images of waves using n = 6 images. First, from known in situ bathymetry and hydrodynamic modeling, we produced a set of simulated images of waves: a pair of images every 1.05 s, i.e., the time lag between Sentinel 2's band 2 and band 4. Second, we used the images of the simulated waves along with our BathySent procedure to retrieve the bathymetry. Third, we compared the in situ bathymetry with the ones retrieved by the BathySent method on the simulated dataset. Last, we applied the BathySent procedure to a real Sentinel 2 dataset of the same study area and compared the results.

Study Area, in situ Bathymetry, and Simulated Waves
To test our method, we first applied it to simulated images of propagating waves. We simulated the wave propagation using a hydrodynamic model, relying on in situ bathymetry.
Our study area was the Gâvres peninsula, French Atlantic coast, France (Figure 1), which was chosen for two reasons. First, because high-resolution in situ bathymetry, with 3-m spatial resolution, built using in situ surveys (for more details, see Idier et al., [20]) were available. Second, the SWASH model (acronym for simulating waves till shore [21]) was able to simulate waves on a wave-by-wave mode, which was already set-up by the authors of the study area [20]. The study area exhibited many different orientations and bathymetric gradients because of many rocky outcrops.   The model was set up with a horizontal resolution of 3 m and a 10-Hz sampling frequency. The model was previously validated for the study area by reproducing the wave overtopping and induced floods that occurred during the Johanna storm (10 March 2008). In this validation, we used a digital elevation model (DEM) representative of the topography and bathymetry of the area during 2008, built using in situ survey data. More details on the model set up, and the in situ bathymetry, can be found in Idier et al. [20]. We ran the model with the following conditions: still water level of 2 m (NGF-IGN69, the French national vertical datum), offshore 2D JONSWAP wave spectrum [22] characterized by a significant wave height (Hs) of 1 m, peak wave period (Tp) of 6 s, peak wave direction (Dp) of 260 • N (nautical convention), and a directional dispersion of 30 • . These offshore conditions were propagated to the open boundaries of the SWASH domain using the spectral wave model WW3 [23], such that the boundary conditions of the SWASH simulations were not uniform, and produced realistic local waves. WW3 is solving the random phase spectral action density balance equation for wavenumber-direction spectra (i.e., the wave spectrum), and thus allows estimating the wave spectrum in space and time. These spectra were used as forcing conditions of the non-hydrostatic phase-resolving SWASH model, which, based on these spectra and the still water level, generated a high frequency water level time series that was subsequently used to solve the nonlinear shallow water equations to provide instantaneous water levels in each grid cell (3 m) and at each time step (10 Hz). We created simulations every 1.05 s, each consisting of a set of six images of propagating waves. In Figure 2, we show the in situ bathymetry and two example simulated wave images.
step (10 Hz). We created simulations every 1.05 s, each consisting of a set of six images of propagating waves. In Figure 2, we show the in situ bathymetry and two example simulated wave images.

Sentinel 2 MultiSpectral Imager
The European Space Agency's (ESA) Copernicus Sentinel 2 mission includes two polar-orbiting satellites that monitor the Earth's surface with a wide swath width (290 km) and high revisit time (10 days at the equator with one satellite, and 5 days with both satellites). Approaching the poles, the revisit time of Sentinel 2 increases where the swaths overlap. When recording a single swath, the Sentinel 2 MultiSpectral Imager (MSI) uses a set of 12 focal plane modules (FPMs) [24]. The key element of being able to measure c relies on the fact that, for technical reasons, multiple CCD bands cannot coexist at the same place in the focal plane of an FPM. Therefore, there is a baseline between two bands due to the difference in the CCD positions, yielding a time lag in image acquisition. The spatial resolution of Sentinel 2's multispectral (MS) mode in the visible range of the electromagnetic spectrum is 10 m. Therefore, it is able to resolve ocean waves from 20-30 m wavelengths onwards. In this study, we used the time span between MSI bands 2 and 4; that is, 1.05 s. In theory, all resolvable waves and celerities from sentinel 2 can estimate water depth in the study area. We used six level 2C Sentinel 2 images provided orthorectified by ESA for the European Copernicus Program, in the Universal Transverse Mercator (UTM, zone 30 N) reference system. We have displayed the Sentinel 2 images and the acquisition dates, centered on the study area, in Figure 3.
The European Space Agency's (ESA) Copernicus Sentinel 2 mission includes two polar-orbiting satellites that monitor the Earth's surface with a wide swath width (290 km) and high revisit time (10 days at the equator with one satellite, and 5 days with both satellites). Approaching the poles, the revisit time of Sentinel 2 increases where the swaths overlap. When recording a single swath, the Sentinel 2 MultiSpectral Imager (MSI) uses a set of 12 focal plane modules (FPMs) [24]. The key element of being able to measure c relies on the fact that, for technical reasons, multiple CCD bands cannot coexist at the same place in the focal plane of an FPM. Therefore, there is a baseline between two bands due to the difference in the CCD positions, yielding a time lag in image acquisition. The spatial resolution of Sentinel 2's multispectral (MS) mode in the visible range of the electromagnetic spectrum is 10 m. Therefore, it is able to resolve ocean waves from 20-30 m wavelengths onwards. In this study, we used the time span between MSI bands 2 and 4; that is, 1.05 s. In theory, all resolvable waves and celerities from sentinel 2 can estimate water depth in the study area. We used six level 2C Sentinel 2 images provided orthorectified by ESA for the European Copernicus Program, in the Universal Transverse Mercator (UTM, zone 30 N) reference system. We have displayed the Sentinel 2 images and the acquisition dates, centered on the study area, in Figure 3.

Methodology
Let us consider a subsection of two MSI bands, normalized by the mean of each subsection amplitude, resulting in two matrices, [A] and [B]. We called this a «window» or resolution cell. Using i and j, we, respectively, identified the columns/lines indices in the window. Over a given window, we computed the two-dimensional discrete fast Fourier transform (DFFT, implemented in NumPy, e.g., [25]) per detector band (matrices [F A ] and [F B ]), where ν and µ were the columns/lines indices in the Fourier domain, respectively. Then, we computed the matrix [R], whose coefficients were given by: where '*' represents the complex conjugate and Am is the amplitude of the produced spectrum.
In the module of R ν,µ , we identified multiple peaks, which were considered representative of the most significant waves. Then, we derived their associated wavelengths, λ ν,µ , and phases, θ ν,µ .
[A] and [B], being separated in time by 1.05 s, were supposed to have the same spatial frequency content. For this reason, the positions of the peaks in |[R]|, |[F A ]|, and |[F B ]| coincided. Therefore, we identified significant wavelengths in [R]. Finally, we determined c by measuring the pixel offset in a wave's direction of motion (δr ν,µ ) between the two bands, which was directly related to the phase of R ν,µ by: Then, we obtained c by dividing Equation (3) by δt k,l , where δt k,l is the time span between band k and band l.
When the pixel offset was not wavelength-dependent (i.e., r ν,µ = r ∀(ν, µ)), Equation (3) was equivalent to the phase correlation algorithm used for earthquake-induced ground offset measurements (e.g., [26]). Now, knowing δt, we utilized the (λ ν,µ , θ ν,µ ) pairs, which were converted to (λ ν,µ , c ν,µ ). The spectral amplitude of R (Am in Equation (2)) was our criterion for selecting the most significant waves in each spectrum (Figure 4), and we kept Am > 0.5 in the normalized spectrum. resolution cell. Using i and j, we, respectively, identified the columns/lines indices in the window. Over a given window, we computed the two-dimensional discrete fast Fourier transform (DFFT, implemented in NumPy, e.g., [25]) per detector band (matrices [F ] and [F ]), where ν and μ were the columns/lines indices in the Fourier domain, respectively. Then, we computed the matrix [R], whose coefficients were given by: where '*' represents the complex conjugate and Am is the amplitude of the produced spectrum.
In the module of R , , we identified multiple peaks, which were considered representative of the most significant waves. Then, we derived their associated wavelengths, λ , , and phases, θ , .
[A] and [B], being separated in time by 1.05 s, were supposed to have the same spatial frequency content. For this reason, the positions of the peaks in |[R]|, | F |, and | F | coincided. Therefore, we identified significant wavelengths in [R]. Finally, we determined c by measuring the pixel offset in a wave's direction of motion (δr , ) between the two bands, which was directly related to the phase of R , by: Then, we obtained c by dividing Equation (3) by δt , , where δt , is the time span between band k and band l.
When the pixel offset was not wavelength-dependent (i.e., r , = ∀( , )), Equation (3) was equivalent to the phase correlation algorithm used for earthquake-induced ground offset measurements (e.g., [26]). Now, knowing δt, we utilized the (λ , , θ , ) pairs, which were converted to (λ , , c , ). The spectral amplitude of R (Am in Equation (2)) was our criterion for selecting the most significant waves in each spectrum (Figure 4), and we kept Am > 0.5 in the normalized spectrum.  . FFT spectra of the simulated waves (left) and Sentinel 2 waves (right). Size of the FFT window in this example: 106 pixels for the simulated waves and 32 pixels for Sentinel 2. (c, λ) couples are here represented as (f, Φ) couples; f is the FFT frequency-proportional to λand Φ is the FFT phase that we used to calculate c. Now, we wanted to extract (c, λ) pairs to find h. We defined small sampling intervals of equal wavelengths within [R]. Within these sampling intervals, we select values with an [R] amplitude above a given threshold. For each of these values, we extracted the correspondent [R] phases. Then, we estimated the average of these phases, weighted by Am. Proceeding in this manner, we obtained the phases associated with the selected wavelengths, and characterized their significance via their amplitude.
In other words, we defined the smaller possible circular sector of [R] containing the wavelength range λ − δλ 2 , λ + δλ 2 . We assumed a nearly constant water depth per FFT window. Therefore, all celerities associated with the circular sector (and thus their phases, as their wavelengths are the same) should be the same. Now that we had all the significant (c, λ) pairs, we estimated one value of h for each (c, λ) pair by solving the dispersion relation (Equation (1)). Afterwards, we now had one value of h for each (c, λ). Here, we took the weighted average of all h-values instead of the best fit to the dispersion relation. This provided us with the final water depth value for a given correlation window or resolution cell, per dataset. The significance of each solution was based on the corresponding amplitude value of [R]. Finally, we corrected each of the N bathymetries for tide offsets and stacked the results into a final bathymetry map. The stacking consists of taking the median value among the six bathymetric maps, for those pixels where there is more than one value. In areas where there is only one value, we keep that value. We propose this way to fill the gaps that might exist if we used one Sentinel 2 image only.

Correction for Tidal Offsets
A bathymetry is the measure of water depth over the study area at one given date. In this study, we used n = 6 images to estimate six bathymetric maps. Each of the six bathymetries, estimated from Sentinel 2 data, provided the water depth, h, at the time of each image acquisition; i.e., the absolute difference between the sea surface level and the sea bottom. Thus, each bathymetry value was provided with reference to the water level, which was a time-varying quantity. We had to correct this offset and estimated the bathymetry (h corrected ) in a fixed vertical reference frame for each Sentinel 2 acquisition date.
Firstly, we estimated the tide level relative to the mean sea level in our study area using the FES2014 tidal component database [27]. Secondly, we estimated the mean sea level relative to Z = 0 NGF-IGN9 using the data provided by the French Hydrografic Service [28]. Thirdly, we combined this information to estimate the tide level (Z tide ) relative to Z = 0 NGF-IGN69 (Table S1, Supplementary Materials). Finally, we subtracted Z tide to total water depth, h, and obtained h corrected . In the present work, we choose the Z = 0 of the NGF-IGN69 national vertical datum.

Results Obtained with the Simulated Dataset
We compiled bathymetry measurements from six simulated wave image couples where each couple had a time lag of 1.05 s and a 3-m pixel size. We used a window size of 32 × 32 pixels, with a sampling interval of 16 pixels. These parameters yielded bathymetry with a 48-m grid resolution. We then calculated the median values of the six bathymetric maps to retrieve an improved map of bathymetry. Figure 5 shows the in situ bathymetry and the results of the BathySent method applied to the simulated dataset. Our results fit the in situ bathymetry very well (Figure 5a,b). To assess the precision and accuracy of our results, we produced scatterplots of the in situ bathymetry versus the BathySent results, and calculated the coefficient of determination (r 2 ) and the related standard deviation (Figure 6). To do this, we down sampled the in situ bathymetry to the same spatial resolution as the BathySent bathymetry. In the scatterplots shown in Figure 6, we discretized Z into a number of intervals of 1-m width, which center takes the values from 0.5 m to 14.5 m. For each interval, we identified Z (in situ) within that interval, calculated the mean, and identified the geographic location of each point in the interval. At those points precisely, we extracted Z (BathySent), and we calculate the mean of Z (BathySent) and its standard deviation. Then, we plotted the mean Z (in situ) versus the mean Z (Bathysent) along with the root mean square error (RMSE) over the interval as a bar (+/− 0.5 RMSE). The total length of the bar was thus the RMSE. The quantitative comparison between the in situ bathymetry and the bathymetry obtained from the BathySent method yielded r 2 = 0.87 with a standard deviation of 2.1 m in the 0-14 m depth zone (Figure 6a). We observe that the BathySent algorithm systematically underestimated deeper values. We provided an explanation in the discussion section. Z (in situ) versus the mean Z (Bathysent) along with the root mean square error (RMSE) over the interval as a bar (+/− 0.5 RMSE). The total length of the bar was thus the RMSE. The quantitative comparison between the in situ bathymetry and the bathymetry obtained from the BathySent method yielded r 2 = 0.87 with a standard deviation of 2.1 m in the 0-14 m depth zone (Figure 6a). We observe that the BathySent algorithm systematically underestimated deeper values. We provided an explanation in the discussion section.  The BathySent bathymetry shown in Figure 5b was coarser with respect to the in situ bathymetry seen in Figure 2a. The sampling step of our methodology yielded a final map whose pixel size was 16 times larger than the original image pixel size. Nevertheless, we can see from Figure 5a,b that the shapes of bathymetric features was generally representative of the geomorphology of the in situ bathymetry. Therefore, we deduced that in extremely favorable conditions, such as the ones present in this exercise (where, for instance, The BathySent bathymetry shown in Figure 5b was coarser with respect to the in situ bathymetry seen in Figure 2a. The sampling step of our methodology yielded a final map whose pixel size was 16 times larger than the original image pixel size. Nevertheless, we can see from Figure 5a,b that the shapes of bathymetric features was generally representative of the geomorphology of the in situ bathymetry. Therefore, we deduced that in extremely favorable conditions, such as the ones present in this exercise (where, for instance, the still water level was known), our methodology is reliable.

Results with Sentinel 2
We chose to use the Sentinel 2, level 1C archive from 2020. Among all the data in the ESA SciHub archive, many contain clouds, and many others do not contain images of waves. As a result, we used six Sentinel 2 datasets (Figure 3) to produce six bathymetry maps with the BathySent method for the study area. For each image, we used a moving window of 32 × 32 pixels (320 m) with a sampling interval of 16 pixels (160 m). A too-large window would imply a decrease in spatial resolution of the results; hence, 320 m was well within the maximum wavelength expected in this area.
The results of applying the BathySent method to real Sentinel 2 data are shown in Figure 5c. We showed that our method estimated bathymetry for most of the study area. The measurements were quite scattered, but they represented a sound and valuable firstlevel information when bathymetry was unknown in a given area. Quantitatively, the comparison between the in situ bathymetry and the BathySent bathymetry yielded an r 2 = 0.7 and a standard deviation of 1.5 m (Figure 6b). Moreover, there is a tendency for the BathySent algorithm to underestimate deeper values. Qualitatively, we showed that the spatial patterns were preserved, but they were not always represented in high fidelity. We interpreted this as the result of the coarser resolution of the bathymetry values retrieved with Sentinel 2 (160 m grid) with respect to the ones retrieved from the synthetic images (48 m grid).

Discussion
We tested the BathySent algorithm to retrieve bathymetry from both simulated and Sentinel 2 images. Firstly, we chose a study area where bathymetry was known from in situ measurements. We started with testing the algorithm on a series of simulated images of waves, to retrieve the known input bathymetry, and the results of the tests were very good in terms of r 2 (>0.8). Next, we applied the algorithm to six real Sentinel 2 images and retrieved six bathymetry maps for the study area. Finally, we stacked the six bathymetry maps to improve the final precision. The process of stacking multiple data and the use of multiple Sentinel 2 images aimed at filling the potential gaps in the final bathymetry map. To evaluate the accuracy in this study, our strategy was not to look at each single bathymetry map but to compare the stacked bathymetry map to in situ bathymetry.
The results from the synthetic images proved that our method was able to retrieve bathymetry at least up to a 16-m depth in the study area. A quantitative comparison between our results and in situ bathymetry showed that our method was able to retrieve bathymetry throughout most of the study area, with a standard deviation of 1.5 m and an r 2 = 0.7. This r 2 result indicated higher dispersion with respect to the synthetic test, which was likely due to a tendency for the BathySent algorithm to underestimate deeper values. This systematic under-estimate of deeper value was expected. It is due principally to two factors. According to equation [1], to extract deeper bathymetry values we would need either very large wavelengths or short wavelengths but a very precise c (precise to the 0.01 m/s). Therefore, given that the maximum nominal precision that we can get on c is 1/10th of the image pixel size, deeper bathymetry values rely on the presence of large wavelengths in the scene (Figure 7). In fact, when we solve Equation (1) if c is overestimated h = nan (not a number). If c is underestimated, h = real but underestimated.  (1). It can be seen that the precision of the measured celerity was much more critical than the wavelength.
As highlighted by Poupardin et al. [15], the quality of c measurements is a key factor in analyses such as presented here. For instance, for a wavelength of 32 m and a celerity of 6 m/s, an error of 0.6 m/s (i.e., 10%) would contribute to a relative error of 42% in the bathymetry. Based on published studies of image correlations with FFT methods (e.g., [26]) we can assume that the precisions of the measured values of c vary between 1/4 and 1/10 the pixel size divided by the time span, depending on the image noise (i.e., between 1 and 2.5 m/s for Sentinel 2). With a smaller pixel size, we could obtain more precise values of c, which in turn would allow more accurate water depth estimations for the same λ. The minimum theoretical detectable λ is twice the pixel size (i.e., 20 m for Sentinel 2), limiting the possibility of investigating short wavelengths in small sub-images close to the shore. In addition, shorter λs are less sensitive to deeper bathymetry than larger λs (Figure 7). Therefore, the theoretical resolution of the final bathymetric map is expected to  (1). It can be seen that the precision of the measured celerity was much more critical than the wavelength.
As highlighted by Poupardin et al. [15], the quality of c measurements is a key factor in analyses such as presented here. For instance, for a wavelength of 32 m and a celerity of 6 m/s, an error of 0.6 m/s (i.e., 10%) would contribute to a relative error of 42% in the bathymetry. Based on published studies of image correlations with FFT methods (e.g., [26]) we can assume that the precisions of the measured values of c vary between 1/4 and 1/10 the pixel size divided by the time span, depending on the image noise (i.e., between 1 and 2.5 m/s for Sentinel 2). With a smaller pixel size, we could obtain more precise values of c, which in turn would allow more accurate water depth estimations for the same λ. The minimum theoretical detectable λ is twice the pixel size (i.e., 20 m for Sentinel 2), limiting the possibility of investigating short wavelengths in small sub-images close to the shore. In addition, shorter λs are less sensitive to deeper bathymetry than larger λs (Figure 7). Therefore, the theoretical resolution of the final bathymetric map is expected to decrease with depth as the detection of larger λs requires larger sub-images. In principle, an optimal algorithm-in terms of final spatial resolution and precision-should be able to simultaneously analyze short and long wavelengths in the image to address these issues.
In the present study, the effect of tidal currents on Sentinel 2 wave celerity is neglected, because it is very small with respect to Sentinel 2 pixel size. Indeed, in the study area, at the exact time of Sentinel 2 image acquisitions, the depth averaged currents-induced by the tide, atmospheric pressure and wind-are estimated by the Ifremer (Institut Français de Recherche pour l'Exploitation de la Mer) MARC service (Modélisation et Analyse pour la Recherche Côtière, [29]) to be lower than 0.1 m/s. The use of the high revisit time and fully available archive of Sentinel 2 was key for selecting multiple images that fulfilled all the requirements described above, which in turn encompassed a diversity of conditions in terms of wave conditions (wavelengths, directions, and visibility in the images) and sensors characteristics (resolution and time spans).

Conclusions
In conclusion, this study confirmed the potential use of multispectral images obtained by Sentinel 2 for shallow bathymetry mapping and it demonstrated the reliability of the BathySent concept.
The high revisit time of Sentinel 2 allowed us to use multiple images without clouds and with visible ocean waves. This point was key to fill the gaps in the bathymetry map that would appear inevitably if we used only one image.
The method has a tendency to underestimate deeper values. This is due to the limit in the precision of c for the range of λs in the study area when we solve the linear dispersion Equation (1). The precision on c depends on the pixel size and on the time lag. Therefore, we conclude that combination of multiple satellite missions, including the VENµS mission [30], Pléiades and SPOT, with a diverse range of pixel sizes and time lags, must be foreseen.