Analysis of Retrackers ’ Performances and Water Level Retrieval over the Ebro River Basin Using Sentinel-3

Satellite altimeters have been used to monitor river and reservoir water levels, from which water storage estimates can be derived. Inland water altimetry can, therefore, play an important role in continental water resource management. Traditionally, satellite altimeters were designed to monitor homogeneous surfaces such as oceans or ice sheets, resulting in poor performance over small inland water bodies due to the contribution from land contamination in the returned waveforms. The advent of synthetic aperture radar (SAR) altimetry (with its improved along-track spatial resolution) has enabled the measurement of inland water levels with a better accuracy and an increased spatial resolution. This study aimed to retrieve water levels from Level-1B Sentinel-3 data with focus on the minimization of the land contamination over smallto middle-sized water bodies (130 m to 4.5 km), where continuous clean waveforms rarely exist. Three specialized algorithms or retrackers, together with a new waveform portion selection method, were evaluated to minimize land contamination in the waveforms and to select the nadir return associated with the water body being overflown. The waveform portion selection method, with consideration of the Digital Elevation Model (DEM), was used to fit the multipeak waveforms that arise when overflying the continental water bodies, exploiting a subwaveform-based approach to pick up the one corresponding to the nadir. The performances of the proposed waveform portion selection method with three retrackers, namely, the threshold retracker, Offset Center of Gravity (OCOG) retracker and two-step SAR physical-based retracker, were compared. No significant difference was found in the results of the three retrackers. However, waveform portion selection using DEM information great improved the results. Time series of water levels were retrieved for water bodies in the Ebro River basin (Spain). The results show good agreement with in situ measurements from the Ebro Reservoir (approximately 1.8 km wide) and Ribarroja Reservoir (approximately 400 m wide), with unbiased root-mean-square errors (RMSEs) down to 0.28 m and 0.16 m, respectively, depending on the retracker.


Introduction
Inland water systems constitute crucial resources of fresh water necessary for human survival [1].Water within rivers and reservoirs represent the primary supply of drinking water, agricultural irrigation and industrial water usage globally [2], and also used to produce renewable hydrological energy.In addition, floods that periodically occur in every region of the world represent threats to crops, settlements, infrastructure, and most importantly life.Therefore, it is important to monitor the water levels of inland water bodies to provide early-warning alerts for water shortages or flood predictions.Unfortunately, in situ gauging stations are not always available in many parts of the world, or not publicly available, being maintained by local authorities [3][4][5].An alternative to this is the use of satellite radar altimeters.Satellite radar altimeters are essential tools for monitoring the oceans, a task they have performed for over 25 years [6].Satellite radar altimetry has also proven to be valuable tools for monitoring the water levels within inland water systems, including lakes and rivers [3,7,8].Satellite radar altimetry can help overcome the lack of data in many parts of the world, and contribute to monitoring the water levels both within inland water systems where no in situ data are available, and also within transboundary river basins.
A list of satellite radar altimetry missions and their resolutions is shown in Table 1.Until CryoSat-2 was launched, the satellite radar altimetry was pulse limited to low-resolution mode (LRM).CryoSat-2 is the first altimetry mission with the synthetic aperture radar (SAR) mode available, followed by Sentinel-3.Unlike classical pulse-limited altimetry, SAR altimetry exploits coherent processing of groups of transmitted pulses to make the most efficient use of the power reflected from the surface [9], therefore significantly improving along-track resolution [9][10][11], as shown in Table 1.Sentinel-3, the newest in-orbit satellite, has a temporal resolution of approximately 27 days and an inclination angle of approximately 98.6 degrees, meaning it can cover almost all Earth latitudes except the poles.Operating in high-resolution SAR mode [10], the ground track separation of Sentinel-3 at the equator is approximately 104 km.In contrast, the CryoSat-2 ground track separation at the equator reaches 7.7 km at the expense of its temporal resolution, at 369 days.Sentinel-3 exhibits the best along-track resolution at 300 m, close to the resolution of CryoSat-2 in SAR and interferometric SAR (SARIn) modes.However, CryoSat-2 has a focus on the cryosphere, and as a result covers most continental surfaces in LRM.Therefore, Sentinel-3 constitutes the first altimetry mission that covers the globe completely in SAR mode, and is thus the most ideal tool for inland water level monitoring thanks to its good global coverage and sufficient temporal and spatial resolutions. 1 Geodetic Satellite. 2European remote sensing satellite. 3Topography Experiment/Poseidon. 4 GEOSAT Follow-On satellite. 5Environmental Satellite. 6Earth Explorer Opportunity Mission. 7Haiyang-2/Ocean-2. 8Satellite with Argos (Data Collection System) and ALtiKa (Altimeter in Ka-band). 9Surface Water Ocean Topography Mission. 10 Jason Continuity of Service Mission.
Numerous studies [33][34][35][36][37][38][39][40][41][42][43] have also focused on middle-sized water bodies.For example, Birkinshaw et al. [33] measured the water level of the Mekong River, the width of which varies from 400 m to 1700 m; their results show RMSEs of approximately 44-65 cm for Envisat and 46-76 cm for European Remote-Sensing Satellite-2 (ERS-2).Da Silva et al. [34] studied water level time series over the Amazon River basin, which exhibits widths ranging from several kilometers to less than a hundred meters, using both Envisat and ERS-2; they showed a RMSE varying from 12 cm in the best cases, 40 cm in most cases, to several meters in the worst cases.Furthermore, Michailovsky et al. [35] employed Envisat to monitor the Zambezi River basin and reported river widths reaching 80 m with RMSEs of 32-72 cm relative to in situ measurements at different locations.In addition, Maillard et al. [36] measured the water levels of medium-sized rivers with widths between 100 and 1000 m over the São Francisco River, Brazil, with RMSEs lower than 60 cm, and better than 30 cm in some cases using Envisat and Satellite with ARgos and ALtiKa (SARAL) data.Most studies of middle-sized water bodies used mainly Envisat and ERS.In addition, previous results referred to LRM altimeters, which provide accuracies on the order of tens of centimeters.RMSEs over different water targets vary with size and morphology, with smaller size and complexity of the shape bring more uncertainties.With the launch of Sentinel-3, equipped with a SAR Altimeter (SRAL), there is a need to evaluate its performance.Because of its higher spatial along-track resolution, SRAL is naturally better suited to eliminate land contaminations found within LRM altimetry footprint [44].Our study objective was to retrieve water levels in challenging environments exploiting SAR altimeter data capabilities.We focused on middle-sized and small-sized water bodies whose width varies from hundreds of meters (close to along-track resolution) to a few kilometers, and/or surrounded by rapid changing topography.
As the backscattered waveform depends on the surface reflecting the signal, altimeter data are usually applied by different retrackers adapted to different surfaces to better locate the height of the reflective surface.A lot of research has been carried out regarding this, and is still ongoing in the retracking modeling framework, specifically considering the evolution from the LRM to SAR altimetric operation [45].There are two types of retrackers, based on either empirical observations and practical experience, or theoretical knowledge of microwave scattering at nadir, namely empirical retracker and physical-based retracker, respectively.Most notable empirical retrackers include the threshold [46], beta-5 [47], and the Offset Center of Gravity (OCOG) [46,48], with the OCOG being the core of ice-1 retracker [49].The physical-based retrackers include the LRM ocean retrackers based on the well-known Brown ocean waveform model [50] (e.g., ice-2 retracker [51,52]), and the altimetry SAR mode ocean retrackers based on the original work of Ray et al. [53] (e.g., the SAMOSA retracker [53,54]).
Whilst well-developed retrackers have been defined for operation over the oceans and ice surfaces, called ocean retrackers and ice retrackers, respectively, no retracker has been adapted for inland water bodies as of yet.Jarihani et al. [15] compared the results from different satellite missions for inland water bodies, where most previous studies used ice retrackers [15][16][17][18]24,34,36], and some others included ocean retrackers [15,16,24].Nevertheless, different retrackers need to be developed and compared for inland water bodies.Due to the relatively recent advance towards SAR altimetry operation, started by the CryoSat-2 mission [55], and continued with the Sentinel-3 [56] mission, most of the studies over inland water are quite limited to the classical LRM operation.When trying to retrack waveforms over continental water bodies, the radar altimetric community has widely applied different retrackers including threshold, ice-1 and ice-2 retrackers.The robust OCOG retracker shows a consistent behavior over the different types of waveforms found in [57] with SARAL/Altika mission over different rivers and reservoirs in India.The SAMOSA retracker and its further adaptations have proven well fitting with CryoSat-2 data not only over the open-ocean and coastal zones [58], but also over ice sheets and inland waters where waveforms show specular characteristics [25,59].
Studies also show the possibilities of retrieving water levels over very small water bodies using techniques that allow precisely monitoring water surfaces.Investigation of individual echoes of Envisat LRM altimetry [60][61][62][63] has shown the possibility of retrieving water levels over water bodies with widths down to 50 m [63].The SARvatore service provides 80 Hz data dedicated to inland water bodies.Another retracker, "fully focused SAR", has shown the potential of monitoring 40 m wide water surfaces by reducing the along-track resolution down to the theoretical limit, equaling half the antenna length [64].The fully focused processing is a very promising processing technique that can provide very high resolution Level-1B products, and has the ability to open a new paradigm in SAR altimetry when it becomes fully operational.
The objective of our study was to develop an operational methodology for water level retrieval from altimetry SAR mode Level-1B waveforms.We focused on small water bodies where continuous clean waveforms rarely exist, and efforts were made to better retrack the contaminated waveforms to improve the accuracy of the water level.First, we modeled three different retrackers including the most widely used threshold retracker, the OCOG retracker, and our in-house SAR physical-based retracker [65] (isardSAT SAR retracker in [65] is an implementation based on the original physical-based model developed by Ray et.al [53], which is also the origin of the sometimes-called SAMOSA model).The latter was integrated in a two-step fitting procedure, exploiting both its normal operation over ocean-like waveforms, and the adaptation to peaky-like waveforms, expected from inland waters.For the second typology, the mean squared slope (MSS) describing the sea surface roughness [66] was used as a fitting parameter instead of the significant wave height (SWH), as proposed and initially evaluated in [59] over Tibetan lakes when using SAMOSA model.(In isardSAT retracker, a rough surface of Gaussian distributions is assumed to define the scattering model, such that the normalized radar cross section (NRCS) as a function of incidence angle is modeled as a Gaussian [67], where the MSS is a parameter defining the directivity of surface radiation pattern, i.e. small values of MSS lead to very narrow and directive radiation patterns related to very specular surfaces.Therefore, considering this parameter (MSS) as fitting term rather than a constant one would allow better fitting specular waveforms.)Additionally, to limit land contamination within the received waveforms, we included Digital Elevation Model (DEM) information to select the waveform portion from nadir, which was performed as a pre-processing in the whole retracking chain.The remainder of this paper is organized as follows.In Section 2, the studied area and the database are presented.In Section 3, the methodologies are introduced, including the selection of the waveform portion, different retrackers and waveform filtering approach.In Section 4, the performances of the different retrackers and the results of the water level time series are shown in comparison with the Sentinel-3 Level-2 ocean retracked results from the European Space Agency (ESA) and in situ measurements.In Section 5, the method and results are discussed.In Section 6, the conclusions are presented.

Study Area
Our study area encompassed the reservoirs and rivers in the Ebro River basin (Figure 1).The Ebro River, which has a length of approximately 910 km and a drainage basin with an area of approximately 85,362 km 2 , is one of the most important rivers on the Iberian Peninsula [68].The river flow is irregular throughout the year, with low levels at the end of summer and high levels during the spring due to melt runoff in the Pyrenees, leading to a danger of flooding.The annual variation of the river flow is more than 800 m 3 s −1 at Ribarroja dam, which is the last dam controlling the water flow entering the Mediterranean Sea [69].The Ebro River is of great importance for agriculture in summer, during which drought often occurs due to the continental Mediterranean climate.Nevertheless, the mean annual flow has decreased by approximately 29% during the 20th century due to many causes: the construction of dams, the increasing demands for irrigation, and evaporation, the latter being higher than the precipitation due to low rainfall amounts, high sunshine intensities and strong, dry winds, from reservoirs within the river basin [68].
from the European Space Agency (ESA) and in situ measurements.In Section 5, the method and results are discussed.In Section 6, the conclusions are presented.

Study Area
Our study area encompassed the reservoirs and rivers in the Ebro River basin (Figure 1).The Ebro River, which has a length of approximately 910 km and a drainage basin with an area of approximately 85,362 km 2 , is one of the most important rivers on the Iberian Peninsula [68].The river flow is irregular throughout the year, with low levels at the end of summer and high levels during the spring due to melt runoff in the Pyrenees, leading to a danger of flooding.The annual variation of the river flow is more than 800 m 3 s −1 at Ribarroja dam, which is the last dam controlling the water flow entering the Mediterranean Sea [69].The Ebro River is of great importance for agriculture in summer, during which drought often occurs due to the continental Mediterranean climate.Nevertheless, the mean annual flow has decreased by approximately 29% during the 20th century due to many causes: the construction of dams, the increasing demands for irrigation, and evaporation, the latter being higher than the precipitation due to low rainfall amounts, high sunshine intensities and strong, dry winds, from reservoirs within the river basin [68].The Pyrenees mountain range lies to the north of the Ebro River basin.Consequently, it is challenging to retrieve the water surface heights of the studied area as the water bodies are relatively small and are greatly influenced by the mountainous terrain.The selected water bodies, for which we have gauging stations for validation, are shown in Figure 1.The coordinates, the widths of the water bodies covered with satellite tracks, the approximate distances between the gauging station and satellite tracks, and the average slope within 5 km along the satellite track are listed in Table 2.The water bodies that can be retracked properly, along with their DEMs, are shown in Figure 2 along with the water masks, Sentinel-3 tracks and the location of gauging stations.Only tracks overpassing the upper stream of the river (before the dam) were considered, ensuring the results could be validated against the in situ measurements.The Pyrenees mountain range lies to the north of the Ebro River basin.Consequently, it is challenging to retrieve the water surface heights of the studied area as the water bodies are relatively small and are greatly influenced by the mountainous terrain.The selected water bodies, for which we have gauging stations for validation, are shown in Figure 1.The coordinates, the widths of the water bodies covered with satellite tracks, the approximate distances between the gauging station and satellite tracks, and the average slope within 5 km along the satellite track are listed in Table 2.The water bodies that can be retracked properly, along with their DEMs, are shown in Figure 2 along with the water masks, Sentinel-3 tracks and the location of gauging stations.Only tracks overpassing the upper stream of the river (before the dam) were considered, ensuring the results could be validated against the in situ measurements.

Sentinel-3
Sentinel-3 is an ocean and land mission based on a constellation of two satellites (Sentinel-3A and Sentinel-3B).Sentinel-3A was launched on 16 February 2016, with data available beginning in June 2016.It was followed by Sentinel-3B, which was launched on 25 April 2018, with data available beginning in December 2018.In this study, two years of Sentinel-3A data were used from June 2016 to June 2018.The SRAL instrument is the main topographic sensor used to provide water level

Sentinel-3
Sentinel-3 is an ocean and land mission based on a constellation of two satellites (Sentinel-3A and Sentinel-3B).Sentinel-3A was launched on 16 February 2016, with data available beginning in June 2016.It was followed by Sentinel-3B, which was launched on 25 April 2018, with data available beginning in December 2018.In this study, two years of Sentinel-3A data were used from June 2016 to June 2018.The SRAL instrument is the main topographic sensor used to provide water level measurements, and, hence, it was used in our research.The detailed parameters of the Sentinel-3 SRAL can be found in [70].To acquire altimeter measurements, the Sentinel-3 SRAL transmits pulses at a Ku-band frequency, which is complemented by a C-band frequency to correct range delay errors due to the varying density of electrons in the ionosphere [70].Sentinel-3 has two operational modes: SAR mode and LRM.As the SAR mode is available globally, we could retrieve inland water levels over any area tracked by Sentinel-3.The SRAL tracks the surface in two different tracking modes, namely, closed loop and open loop tracking.For closed loop tracking, the altimeter range window is autonomously positioned based on an on-board near real-time (NRT) analysis of the previous SRAL waveform; in contrast, for open loop tracking, the altimeter range window is positioned using a priori knowledge of the surface height stored in the instrument in a one-dimensional along-track DEM.The tracking modes for the studied water bodies over the Ebro River basin are listed in Table 2.
Three levels of processed altimeter data are available: Level-0, Level-1 and Level-2 products.In our study, Level-1 non-time critical (NTC) 20 Hz data were used for the water level retrieval by the three retrackers (i.e., the threshold, OCOG and two-step physical-based retrackers), and Sentinel-3 Level-2 ocean retracked data from ESA were used for comparison.
The main objectives of the Level-2 processing of SAR mode data are to provide elementary retracked altimeter estimates of the oceans, coastal zones, ice sheets and sea ice elevations [71].Different retracking algorithms are more suited to specific surfaces, that is, for ocean retracking, OCOG retracking, ice sheet retracking, ice retracking and sea ice retracking.Unfortunately, ice-related retracker results are not available for inland water bodies; hence, we used Sentinel-3 Level-2 ocean retracker results for comparison.Ocean retrackers try to fit theoretically modeled multi-look Level-1B (L1B) waveforms to real L1B SAR waveforms, thereby providing estimates of the epoch, composite sigma, amplitude and the mispointing angle [71].

In Situ Data
In situ data for the Ebro River basin are available in the Automatic Hydrological Information System, in Spanish known as the "Sistema Automático de Información Hidrológica", or the SAIH Ebro data hub [72].SAIH Ebro is an online system providing hourly and daily hydrological information, including river gauge data, reservoir levels, rainfall amounts and temperatures, over the Ebro River basin.Together with the Sentinel-3 passes, we studied eight water bodies in our research, as listed in Table 2.The in situ data were collected from June 2016 to June 2018 on an hourly basis.

Digital Elevation Model
A DEM was used as an ancillary dataset for the selection of the waveform portion.The Shuttle Radar Topography Mission (SRTM) is an international research effort that obtains DEM data on a near-global scale from 56 • S to 60 • N [73].The SRTM provides global data at two resolutions: 1 arc-second (~30 m) and 3 arc-seconds (~90 m).In our study, 1 arc-second global elevation data, which offer a worldwide coverage of void-filled data at a resolution of approximately 30 m, were used.The DEMs of water bodies are shown in Figure 2. The vertical accuracy of SRTM DEM is about ±16 m (absolute) and ±6 m (relative) [73][74][75][76][77].However, vertical accuracy of the data decreases with the increase in slope and elevation due to presence of large outliers and voids [75].

Geophysical Corrections
The space-borne radar altimeter is an essential tool for monitoring the oceans, but it can also be used for inland water surfaces, including lakes and rivers.The principle of altimetry can be found in [78] with application of geophysical corrections.In this study, corrections for the wet troposphere, dry troposphere, ionosphere, solid earth tide, geocentric pole tide and ocean loading Remote Sens. 2019, 11, 718 9 of 25 tide were considered, as discussed by Fernandes et al. [79], and the geoid correction was applied.As the objective of the study was to directly exploit official Sentinel-3 L1B data and compare against Level-2 official products, the use of the geophysical corrections available in the Level-2 product were considered.The detailed corrections and their ranges are listed in Table 3 [80].The Geoid model used is the Earth Gravitational Model 2008 (EGM2008) [81].To improve the altimeter range accuracy, which is related to the water level measurement accuracy, the waveform needs to be retracked precisely to determine the accurate tracking point located on the leading edge [87].We tested three different retrackers, which are listed in Section 3.2, to find the tracking points precisely from the land-contaminated waveforms.

Threshold Retracker
The threshold retracker is a simple retracker based on an estimation of the epoch (the leading edge position) as a percentage of the maximum peak [46,88].In principle, it works well when the maximum peak originates from a nadir water body.The epoch is calculated as the first sample or range bin that is above a percentage of the waveform peak.To provide a higher precision on the estimated epoch, a linear interpolation between the range bins adjacent to the threshold crossing is used as suggested in [46]: where n th is the first range bin of the waveform whose power, y(n th ) is right above the threshold (A peak • η th ), A peak is the amplitude of the peak of the waveform, and η th is defined as a threshold percentage (50% in our case).The term A peak • η th − y(n th − 1) /(y(n th ) − y(n th − 1)) in Equation ( 1) corresponds to the linear proportionality coefficient.

Offset Center of Gravity (OCOG) Retracker
The OCOG retracker [46,48] was designed to calculate the center of gravity of the reflected waveform based on the power levels in the bins.It is an empirical retracker that implements a combination of the center of gravity (COG) and a conventional offset to estimate the related epoch (the leading edge position).Three main parameters are estimated for the OCOG retracker: the COG, the window size (W) and the amplitude (A).
Frappart's definition [24] was considered in our analysis using squares of the power samples: where y(n) is the n th power sample of the input waveform, and n 1 and n 2 are the initial and end range bins or positions, respectively, in the waveform used for the OCOG parameter estimation.
The estimated epoch can then be implemented using the conventional OCOG definition of the epoch: The in-house isardSAT SAR ocean retracker in [65] is integrated in a two-step fitting procedure to operate over inland waters, providing robust surface height estimation with a minimal modification to the SAR ocean retracker model.The waveforms reflected from an inland water body might exhibit ocean-like shapes (winds blowing over lakes) or, most commonly, a peakier shape (smoother surfaces) for small water bodies.To adapt to the latter ones, the MSS, describing the sea surface roughness [66], was used instead of the SWH as a fitting parameter; this approach is initially considered in [59], when exploiting the SAMOSA model over Tibetan plateau lakes.
The SAR ocean retracker in [65] considers a scattering model of rough surface with Gaussian statistics, and so the normalized backscattering coefficient (surface radiation pattern) as function of the incidence angle θ inc can be modeled as [67]: where the MSS represents the variance of the surface slopes.From the relationship in Equation ( 4), it can determined that, when the variance of the surface slopes (MSS) is large, a broad surface radiation pattern is obtained (rough surfaces), while, when MSS tends to small values, the pattern is more directed as expected for specular returns (smooth surfaces).Therefore, the MSS can be accordingly adjusted as a fitting parameter to properly fit specular returns.With the objective of exploiting the SAR ocean retracker [65] with minimal changes over inland waters, the same retracker model can be used, but the SWH is fixed (to a value close to zero), while the MSS or surface roughness is fitted.A similar approach with the SAMOSA model is exploited in [89] to fit waveform returns from leads.
To ensure that the different typology of waveforms (ocean-and specular-like) can be fitted, a two-step retracking is implemented.In a first run, the normal operation of the physical-based retracker in [65] was considered, i.e., the SWH was one of the fitting parameters (keeping MSS constant).Then, a second run of the fitting was implemented but this time the SWH was fixed (to a very small value around 1e-5) and the MSS was fitted so that the radiation pattern of the rough surface was adjusted.Then, the estimate of the surface height was extracted from the fit with the best correlation coefficient between the measured waveform and the modeled waveform.In this way, it was ensured that peaky-like waveforms could be properly fitted using the MSS as a fitting parameter, since the normal SAR ocean retracker would fail to adjust the model.An example of the operation of the two-step retracker is shown in Figure 3 over a peaky waveform.
The in-house isardSAT SAR ocean retracker integrated in the two-step procedure [65] is an implementation based on the SAR ocean model proposed by Chris Ray et al. [53].It is a physically-based model of the altimetric power waveform over the ocean received by the altimetric sensor after delay-Doppler (or SAR) processing.Details regarding the definition and implementation of this retracker can be found in [65].
In SAR or delay-Doppler altimetry, different Doppler beams, also known as looks (l), are synthesized to focus on a given surface on the ground [65].Thus, averaging of all these Doppler power beams would lead to the final multilook power waveform to be retracked.Therefore, for this physical-based retracker, each single-look power waveform was modeled as follows: where k refers to the range index vector within the received waveform window, l is the Doppler index or beam pointed to the specific surface, P u is the amplitude, k 0 is the epoch (related to the surface height), SWH is the significant wave height (referred also as H s ), and MSS is the mean square slope related to the surface roughness.In the first iteration of the two-step retracker, P u , k 0 , SW H are the fitting parameters, while in the second iteration the parameters being adjusted are P u , k 0 , MSS.The in-house isardSAT SAR ocean retracker integrated in the two-step procedure [65] is an implementation based on the SAR ocean model proposed by Chris Ray et al. [53].It is a physically-based model of the altimetric power waveform over the ocean received by the altimetric sensor after delay-Doppler (or SAR) processing.Details regarding the definition and implementation of this retracker can be found in [65].
In SAR or delay-Doppler altimetry, different Doppler beams, also known as looks (), are synthesized to focus on a given surface on the ground [65].Thus, averaging of all these Doppler power beams would lead to the final multilook power waveform to be retracked.Therefore, for this physical-based retracker, each single-look power waveform was modeled as follows: where  refers to the range index vector within the received waveform window,  is the Doppler index or beam pointed to the specific surface,  is the amplitude,  is the epoch (related to the surface height), SWH is the significant wave height (referred also as  ), and MSS is the mean square slope related to the surface roughness.In the first iteration of the two-step retracker,  ,  ,  are the fitting parameters, while in the second iteration the parameters being adjusted are  ,  , .
The terms  , and  , encompass the information related to the antenna pattern, antenna mis-pointing, and surface radiation patterns.They are obtained, respectively, as the constant and first order terms of a Taylor approximation of the antenna and surface radiation pattern's product.The closed expression for these two terms can be found in [65].A Gaussian-like antenna and surface radiation (as indicated in Equation ( 5)) patterns were assumed.
In Equation ( 5),  refers to the Doppler-dependent () dilation term, and takes into account The terms B k,l and T k,l encompass the information related to the antenna pattern, antenna mis-pointing, and surface radiation patterns.They are obtained, respectively, as the constant and first order terms of a Taylor approximation of the antenna and surface radiation pattern's product.The closed expression for these two terms can be found in [65].A Gaussian-like antenna and surface radiation (as indicated in Equation ( 5)) patterns were assumed.
In Equation ( 5), g l refers to the Doppler-dependent (l) dilation term, and takes into account the instrument configuration and the significant wave height: where σ ac and σ al refer to the widths of the Gaussian functions that approximate the point target response in the across-track (or range) and along-track (azimuth) dimensions, respectively, and σ z corresponds to the standard deviation of the sea surface height probability density function (PDF) and is related to the SWH as σ z = H s 4 .The definition of the additional terms in Equation ( 6) can be found in Table 4, where the main parameters of the SAR ocean model are summarized.
Finally, f 0 and f 1 represent the range-dependent functions modulated by the Doppler-dependent dilation term and can be obtained from the general expression (for n = 0 and 1): Standard deviation of the height PDF 1 For a chirp pulse the bandwidth can be expressed in terms of the chirp rate (k r ) and pulse duration (T p ) as follows: BW = K r • T p .

Waveform Exclusion
For very small water bodies, the waveform may be highly contaminated by the land.Thus, we added additional conditions to exclude the waveforms with a poor quality.
First, if the epoch of the reference SRTM DEM is outside the waveform window, which has a measuring range of approximately 60 m, it is likely that the on-board tracker window was unable to focus on the reflecting surface.Second, within the vicinity of land, reflections from both water and land will contribute to the received waveform, resulting in more than one peak.If the number of outstanding peaks is larger than a threshold, the waveform may be contaminated substantially by land or contains no water information at all, making it not suitable for water level retrieval.In our case, the threshold was set to 5 based on the analysis of waveforms.The number of outstanding peaks is normally less five, even with land contaminations.Finally, sigma_0, which is the backscatter coefficient for water from nadir, should be much larger than the backscatter coefficient for land, where the value is normally less than 50 dB based on the analysis in our study area.With these conditions (summarized below), we could filter out highly contaminated waveforms that were unsuitable for detecting the levels of water bodies:

•
The epoch of the reference SRTM DEM must be within the waveform window After waveform exclusion, only the remaining waveforms (including both slightly contaminated and clean waveforms) were considered in our study to maximize the usage of the available measurements in the water body of interest.For a small water body, the quality of the waveform largely depends on the satellite tracks which drift daily, and the surface area changes with seasons, resulting in no clean waveform from time to time.

Selection of the Waveform Portion
Over inland water bodies, the on-board tracker may not properly locate the retrieval window to acquire the signal from the nadir corresponding to the water body.In addition, specific across-track surface contamination will be present in the waveforms, making it difficult to properly track the desired portion of the waveform.In our study, we used DEM information (SRTM DEM at a resolution of 30 m) to locate the waveform portion that comes from nadir within the receiving window, and then retracked the waveform using the selected portion only.This pre-processing approach, which prepares the portion of waveforms, can be applied and integrated with any kind of retracker.
The waveform portion selection method isolates the nadir return within the receiving window as follows: • From the geo-located surface of interest within the water body and beneath the satellite track, the associated height was interpolated from the DEM information and referred to the geodetic ellipsoid, H DEM .

•
This height was subtracted from the satellite height (H orb ) at the geo-located surface to obtain a rough estimation of the range (or equivalently window delay) from where the nadir returns were expected.

•
This range was linked to a specific bin within the received waveform, by comparing it with the vector of ranges associated to each bin in the receiving window, which could be obtained from the measured range by the radar and the range sampling.

•
The peak location closest to the previous range bin position was taken, and the portion of the waveform was selected around this peak considering the valley positions to the left and right of the peak plus some guard samples on top: A built-in Matlab function (findpeaks) was used to compute the prominent or outstanding peaks within the waveform.Prominent peaks are those peaks that drop more than a given threshold value on either side of the peak before the signal attains a higher value.The associated valley locations can be extracted using this built-in function, but taking as input the maximum of the whole waveform minus the waveform itself.
The water level for each footprint was retrieved based on the selected portion of the waveform applying the three retrackers.The time series of the water levels were calculated using a strict water mask polygon as shown in Figure 2, within which the water levels of the altimeter footprints were considered, selected and averaged for each track.The water mask polygon was created manually according to the shoreline displayed in the Google Earth image.
The workflow of the water level retrieval using the Sentinel-3 Level-1 product included four main steps: waveform exclusion, selecting the waveform portion using a DEM determining the retracking range using different retrackers, geophysical corrections and water level averaging, as shown in Figure 4. Geophysical corrections were applied to the retracked range using the Sentinel-3 Level-2 product, and the accuracy was assessed for each water body using in situ measurements from the gauging stations in the SAIH system.

Selecting the Waveform Portion Using a DEM
When a waveform is influenced by land-based contamination, the waveform will contain small peaks due to off-nadir land or land coverage that follow the large peak, which corresponds to the water reflection within the altimeter footprint.In this case, as shown in Figure 5 (two outstanding peaks in Figure 5 waveform), the portion of the selected waveform provided a more precise tracking point for the two-step physical-based retracker and OCOG retracker, especially for the OCOG retracker, whose retracking point was the front edge of the COG window.The footprint of the signal was very close to the shore, thereby introducing a large peak at the trailing edge.The DEM information was used to select the portion of the waveform closer to nadir.The closest peak and valley after nadir were selected for retracking.Figure 5a,c shows the location of nadir in a black vertical line.The portion being selected is under the SWH fitting or MSS fitting line in Figure 5a with the two-step physical-based retracker, and is marked out in red crosses as the L1B-Waveform portion in Figure 5c with OCOG retracker.For the two-step physical-based retracker, the fitting using the selected waveform portion (Figure 5a) resulted in a 99% accuracy, while the fitting using the whole waveform (Figure 5b) resulted in an accuracy of 83%, making the location of the tracking point on the leading edge slightly different.In this study, the objective of the two-step physical-based retracker was to ensure proper retrieval of the surface height over inland waters, when exploiting the SAR ocean retracker in [65] with minimal changes on its operation; the

Selecting the Waveform Portion Using a DEM
When a waveform is influenced by land-based contamination, the waveform will contain small peaks due to off-nadir land or land coverage that follow the large peak, which corresponds to the water reflection within the altimeter footprint.In this case, as shown in Figure 5 (two outstanding peaks in Figure 5 waveform), the portion of the selected waveform provided a more precise tracking point for the two-step physical-based retracker and OCOG retracker, especially for the OCOG retracker, whose retracking point was the front edge of the COG window.The footprint of the signal was very close to the shore, thereby introducing a large peak at the trailing edge.The DEM information was used to select the portion of the waveform closer to nadir.The closest peak and valley after nadir were selected for retracking.Figure 5a,c shows the location of nadir in a black vertical line.The portion being selected is under the SWH fitting or MSS fitting line in Figure 5a with the two-step physical-based retracker, and is marked out in red crosses as the L1B-Waveform portion in Figure 5c with OCOG retracker.For the two-step physical-based retracker, the fitting using the selected waveform portion (Figure 5a) resulted in a 99% accuracy, while the fitting using the whole waveform (Figure 5b) resulted in an accuracy of 83%, making the location of the tracking point on the leading edge slightly different.In this study, the objective of the two-step physical-based retracker was to ensure proper retrieval of the surface height over inland waters, when exploiting the SAR ocean retracker in [65] with minimal changes on its operation; the significance and interpretation of other geophysical parameters that Remote Sens. 2019, 11, 718 can be inferred (SWH and MSS) needs to be further investigated and it is out of the scope of the current study.
For the OCOG retracker, the COG when using the selected waveform portion was located inside the first peak.However, when using the whole waveform, it was located in between peaks, resulting in different tracking points that were the front edge of the COG window.For the threshold retracker, the tracking point did not change in this case.The threshold retracker depended only on the largest peak, which, as shown in Figure 5, was the same with or without waveform portion selection.However, in the case of the largest peak not being from the nadir water body but from the surroundings or other objectives, the threshold retracker with and without waveform portion selection led to different results.
Remote Sens. 2019, 11, x FOR PEER REVIEW 15 of 26 significance and interpretation of other geophysical parameters that can be inferred (SWH and MSS) needs to be further investigated and it is out of the scope of the current study.
For the OCOG retracker, the COG when using the selected waveform portion was located inside the first peak.However, when using the whole waveform, it was located in between peaks, resulting in different tracking points that were the front edge of the COG window.For the threshold retracker, the tracking point did not change in this case.The threshold retracker depended only on the largest peak, which, as shown in Figure 5, was the same with or without waveform portion selection.However, in the case of the largest peak not being from the nadir water body but from the surroundings or other objectives, the threshold retracker with and without waveform portion selection led to different results.The differences between the performance using the waveform portion selection method, and the method using the whole waveform, are shown in the time series results over the Ebro Reservoir in Figure 6.The Level-2 data were obtained from the Sentinel-3 Level-2 product using the ocean retracker, which does not include any waveform portion selection.The waveform portion selection method eliminated most outliers compared with the method of using the whole waveform.The differences between the performance using the waveform portion selection method, and the method using the whole waveform, are shown in the time series results over the Ebro Reservoir in Figure 6.The Level-2 data were obtained from the Sentinel-3 Level-2 product using the ocean retracker, which does not include any waveform portion selection.The waveform portion selection method eliminated most outliers compared with the method of using the whole waveform.

Time Series Validation
The time series of the water levels were calculated using a strict water mask polygon.The water levels of the altimeter footprints within the mask were considered, selected and averaged for each date.Figures 6a and 7 show the water levels for different water bodies derived using the three retrackers: the two-step physical-based retracker, OCOG retracker and threshold retracker, combined with the waveform portion selection method.The retracked results were compared with the in situ measurements and the Level-2 results using the ocean retracker.

Time Series Validation
The time series of the water levels were calculated using a strict water mask polygon.The water levels of the altimeter footprints within mask were considered, selected and averaged for each date.Figures 6a and 7 show the water levels for different water bodies derived using the three retrackers: the two-step physical-based retracker, OCOG retracker and threshold retracker, combined with the waveform portion selection method.The retracked results were compared with the in situ measurements and the Level-2 results using the ocean retracker.
Among the eight water bodies monitored in this study, the water levels over the Mequinenza Reservoir, Cavallers Reservoir and San Salvador Reservoir, all of which were tracked in an open loop tracking mode (on-board DEM dependent), could not be retracked.The difference between the on-board tracking heights derived from the received waveforms and the heights from the SRTM DEM with reference to the geodetic ellipsoid was almost 50 m over the Mequinenza Reservoir, more than 1000 m over the Cavallers Reservoir, and approximately 80 m over the San Salvador Reservoir.These findings indicate that the SRAL sensor lost its track and was unable to acquire the waveforms reflected from these water body surfaces.Among the eight water bodies monitored in this study, the water levels over the Mequinenza Reservoir, Cavallers Reservoir and San Salvador Reservoir, all of which were tracked in an open loop tracking mode (on-board DEM dependent), could not be retracked.The difference between the on-board tracking heights derived from the received waveforms and the heights from the SRTM DEM with reference to the geodetic ellipsoid was almost 50 m over the Mequinenza Reservoir, more than 1000 m over the Cavallers Reservoir, and approximately 80 m over the San Salvador Reservoir.These findings indicate that the SRAL sensor lost its track and was unable to acquire the waveforms reflected from these water body surfaces.
The RMSE and unbiased RMSE (ubRMSE) (Equation ( 8)) were calculated, as listed in Table 5.The unbiased RMSE is more reliable as it removes the bias.

𝑢𝑏𝑅𝑀𝑆𝐸 = 𝐸{[(𝑊𝐿 − 𝐸[𝑊𝐿]) − (𝑖𝑛𝑠𝑖𝑡𝑢 − 𝐸[𝑖𝑛𝑠𝑖𝑡𝑢])] }
where E[•] is the expectation operator and WL is the water level derived from satellite data.The RMSEs for different water bodies varied from 17 cm to more than 1 m.A comparison of the performances for the different retrackers over different water bodies is shown in  The RMSE and unbiased RMSE (ubRMSE) (Equation ( 8)) were calculated, as listed in Table 5.The unbiased RMSE is more reliable as it removes the bias.
where E[•] is the expectation operator and WL is the water level derived from satellite data.The RMSEs for different water bodies varied from 17 cm to more than 1 m.A comparison of the performances for the different retrackers over different water bodies is shown in

Discussion
In this study, the size of the target water bodies vary from 130 m to 4.5 km and change with seasons.The geometries of the water bodies are complex with irregular lake shores and various sizes, making the along-track width (the width as seen by the satellite flying along the water body) very different from one satellite overpass to another due to the drifting satellite ground tracks.As a result, continuous clean waveforms rarely exist in our study area.
The shore of the Ebro Reservoir is irregular, and lands exist within water, resulting in no clean waveforms available over many tracks.Figure 6b shows the water level time series derived using the whole waveform, while Figure 6a shows the water level time series derived using the waveform portion selection.Retracking of the waveform portion accordingly reduced land contamination effects on the waveform.Figure 6a shows higher consistency with in situ time series, and thus shows that most of the outliers were caused by the land contamination.The land contaminated waveforms can be originated by surface area changes or satellite track drifting, making the altimeter footprint closer to the side bank, introducing outliers in the time series.
Over the water bodies where Sentinel-3 is operated in closed loop, the water levels of all water bodies could be tracked.Over the water bodies where Sentinel-3 is operated in open loop, three water bodies were missed and only one was tracked accurately.The results over the Ebro Reservoir, Sotonera and Ribarroja Reservoirs are satisfactory.The Ebro Reservoir is located in the north near the Pyrenees mountain range with elevations exceeding 820 m, and an average slope within 5 km along track of about 4%.The RMSE was 30 cm compared with the in situ measurements.The Sotonera Reservoir is the biggest reservoir studied, and is located in a relatively flat area with a slope of about 3% and an elevation of about 410 m.The result is good as there is no influence from the mountainous terrain, and it is big enough to be well monitored by satellite.The Ribarroja Reservoir is only 400 m wide and with a slope of more than 20%, but an elevation of only around 70 m.Its small width results in only one altimeter footprint per track available within the water body.The ubRMSE of Ribarroja Reservoir shows very good performance with a value of about 16 cm.However, the results for a few tracks were excluded due to poor waveform quality, which may be influenced by the rapid changes of surrounding terrain.The RMSEs over the Itoiz Reservoir and Irabia Reservoir exceeds 1 m.One reason for this could be the rapid change in the elevations along the satellite track as both water bodies are located in the Pyrenees, making it easy for a satellite to lose its track.The width of the Itoiz Reservoir with satellite data varies from 400 m to more than 2 km depending on the different satellite tracks, which drift daily.In addition, the Irabia Reservoir is a very small water body with a width of approximately 130 m, making the retrieval of its surface height more challenging.Besides the waveform quality, the size of the water body and the surrounding terrain, another reason for the differences with in-situ data from the gauging station is the distance between the gauging station and the Sentinel-3 tracks, varying from a few meters to several kilometers.Satellite tracks are always drifting (up to 1.5 km), making the footprints always changing in their locations for each water body.With the altimetry mission originally designed for ocean, the accuracy of geophysical corrections over inland water ranges from a few centimeters to tens of centimeters, depending on the size of the water body and wind conditions [79,90,91].These effects will most definitely influence the final results, and need to be further considered quantitatively.
The results of the three retrackers, in combination with the waveform portion selection method, always give better accuracies than the results of the Sentinel-3 Level-2 ocean retracker (the results of both sets of retrackers are shown in Figure 6 and Table 5).Therefore, given the lack of an ice retracker over inland water bodies, retrieving the water levels by retracking the Level-1 waveforms and using the waveform portion selection method constitutes an alternative approach to obtain more accurate water levels.
The waveform portion selection method using DEM information greatly improved the results shown in both Figures 5 and 6.As mentioned above, the vertical accuracy of SRTM DEM is approximately ±16 m (absolute) and ±6 m (relative).In the case of a 6 m error, the difference of choosing the right peak corresponding to nadir would be about 3 m.Therefore, if two large peaks were very close to each other, the portion of the waveform chosen might not be accurate.However, for most cases, the DEM error is small compared to the waveform tracking window, which is more than 60 m, and thus optimistic results can be obtained.Besides, waveform portion selection using DEM information improved the water level results over water bodies where only a few or no clean waveforms are available, and retained more results to have a better complete time series.In the case of the water level changing enormously and the waveform being heavily contaminated at the same time, the waveform portion selection approach may lose its capability of precise retracking.
The different retrackers obtained a good accuracy with in situ measurements.The ubRMSEs showed smaller values compared with RMSEs for all retrackers over different reservoirs, as listed in Table 5.The bias between the retrieved water levels and the in situ measurements existed in most cases, which may be related to the drift of the river course (the measurements are not co-located with the Sentinel-3 track).Other error sources include the waveform quality, which depends on the surrounding terrain.These include the geometry of the lake shore and the characteristics of the terrain, all of which contribute to the shape of the waveform.Nevertheless, Sentinel-3 has been proven to work for small water bodies.
We detected some issues related to the Sentinel-3 SRAL tracking mode.The open loop tracking mode, which depends on the on-board DEM to locate the tracking window, loses its track over the Cavallers Reservoir, which is surrounded by mountains, but also in Mequinenza and San Salvador Reservoirs, which are in relatively flat areas.We provided the correct water level height information to make Sentinel-3 able to capture the correct water level in the open loop tracking mode using the service altimeter open loop tracking command for hydrology monitoring (https://www.altimetry-hydro.eu/).We optimistically expect the water bodies that are under the Sentinel-3 coverage, but with loss of track, will be monitored in the near future.Another potential solution could be changing the tracking mode to a closed loop mode.

Conclusions
The objective of this study was to develop a methodology that can estimate water levels in small water bodies (similar size to the along-track resolution) and/or complex topography environments where many waveforms are expected to be contaminated by land.The study compared the performances of the threshold, OCOG and two-step physical-based retrackers, but no significant differences were found in the results of the three retrackers.However, the main improvements found show that selecting the portion of the waveform to be fitted significantly affected the retrievals.
The novel DEM-oriented waveform portion selection method could isolate the nadir peak from land-based contamination relatively well, and greatly improved the results with a reduction of RMSE by more than 0.5 m over Ebro Reservoir.Retracking the water levels from Level-1 waveforms using the combination of a retracker and the waveform portion selection method was more robust to land contamination than using Sentinel-3 Level-2 data directly.It also resulted in a better accuracy.
The water levels were compared with in situ measurements and the Sentinel-3 Level-2 ocean retracker from ESA.The results show good agreement with the in situ measurements.The ubRMSE over the Ribarroja Reservoir with a width of about 400 m could reach 16 cm, while that over the Ebro Reservoir whose width is about 1.8 km, it could reach 28 cm.Waveform portion selection using DEM information showed its ability to retrieve water levels over small-to medium-sized water bodies.In contrast, the accuracy over water bodies located in mountainous areas still needs to be improved, but the results demonstrate the possibility of retrieving the water levels over a very small water body with a width of approximately 130 m.
Overall, the Sentinel-3 SRAL has been proven to work over inland water bodies, even those with widths as small as 130 m.Retracking using the waveform portion selection method with an SRTM DEM as a pre-processing to the two-step physical-based retracker, the OCOG retracker and threshold retracker improved the accuracy with an optimal ubRMSE of 16 cm.Further steps need to be taken to explore the possibilities for both tracking water bodies that are not currently being tracked, and for improving the accuracies for very small water bodies with Sentinel-3 altimeter data.

Figure 1 .
Figure 1.The Ebro River basin: (a) the Ebro River network in the Iberian Peninsula; and (b) the water bodies covered by Sentinel-3 satellite tracks with all available gauging stations.

Figure 1 .
Figure 1.The Ebro River basin: (a) the Ebro River network in the Iberian Peninsula; and (b) the water bodies covered by Sentinel-3 satellite tracks with all available gauging stations.

Remote 26 Figure 3 .
Figure 3. Operation of the two-step physical-based retracker over peaky-like waveform.r (1st) is the correlation coefficient for the SWH fit with the original waveform, while r (2nd) is the correlation coefficient for the MSS fit.

Figure 3 .
Figure 3. Operation of the two-step physical-based retracker over peaky-like waveform.r (1st) is the correlation coefficient for the SWH fit with the original waveform, while r (2nd) is the correlation coefficient for the MSS fit.

Figure 4 .
Figure 4. Overview of the workflow.

Figure 4 .
Figure 4. Overview of the workflow.

Figure 5 .
Figure 5.Comparison of the waveform portion selection method with the method using the whole waveform for the two-step physical-based retracker and OCOG retracker: (a) two-step physical-based retracker with the selected waveform portion; (b) two-step physical-based retracker with the whole waveform; (c) OCOG retracker with the selected waveform portion; and (d) OCOG retracker with the whole waveform.Red stars mark the selected portion of the waveform.r (1st) is the correlation coefficient for the SWH fit with the original waveform, and r (2nd) is the correlation coefficient for the MSS fit.

Figure 5 .
Figure 5.Comparison of the waveform portion selection method with the method using the whole waveform for the two-step physical-based retracker and OCOG retracker: (a) two-step physical-based retracker with the selected waveform portion; (b) two-step physical-based retracker with the whole waveform; (c) OCOG retracker with the selected waveform portion; and (d) OCOG retracker with the whole waveform.Red stars mark the selected portion of the waveform.r (1st) is the correlation coefficient for the SWH fit with the original waveform, and r (2nd) is the correlation coefficient for the MSS fit.

Figure 6 .
Figure 6.Time series water levels over the Ebro Reservoir: using the waveform portion selection method (a); and using the whole waveform (b).

Figure 6 .
Figure 6.Time series water levels over the Ebro Reservoir: using the waveform portion selection method (a); and using the whole waveform (b).

Table 4 .
[65]nition of the main parameters of the SAR ocean model[65].

Table 5 .
The best result was obtained for the Ribarroja Reservoir, which had an ubRMSE of approximately 16 cm.The results over the Ebro Reservoir, Sotonera and Ribarroja Reservoirs were also satisfactory.The RMSEs for the Itoiz Reservoir and Irabia Reservoir were relatively large with values greater than 1 m.The results using the three different retrackers do not show large differences except for the Sotonera Reservoir and Ribarroja Reservoir.The OCOG retracker gave the smallest ubRMSE for the Sotonera Reservoir with a value of approximately 38 cm.The two-step physical-based retracker and threshold retracker performed better over the Ribarroja Reservoir than the OCOG retracker, which had an ubRMSE of approximately 16 cm.All retrackers gave better results than the Level-2 ocean retracker without the waveform portion selection method.

Table 5 .
The best result was obtained for the Ribarroja Reservoir, which had an ubRMSE of approximately 16 cm.The results over the Ebro Reservoir, Sotonera and Ribarroja Reservoirs were also satisfactory.The RMSEs for the Itoiz Reservoir and Irabia Reservoir were relatively large with values greater than 1 m.The results using the three different retrackers do not show large differences except for the Sotonera Reservoir and Ribarroja Reservoir.The OCOG retracker gave the smallest ubRMSE for the Sotonera Reservoir with a value of approximately 38 cm.The two-step physical-based retracker and threshold retracker performed better over the Ribarroja Reservoir than the OCOG retracker, which had an ubRMSE of approximately 16 cm.All retrackers gave better results than the Level-2 ocean retracker without the waveform portion selection method.

Table 5 .
Comparison of the water level validation results over all monitored water bodies.