Werner Alpers
email: werner.alpers@uni-hamburg.de

Yuan Zhao

Alexis Mouche

Wai Pak

Chan

Alpers Werner

Zhao Yuan

Mouche Alexis

Chan Pak

A note on radar signatures of hydrometeors in the melting layer as inferred from Sentinel-1 SAR data acquired over the ocean

Keywords: Rain cells, C-band radar backscatter, Melting layer, Weather radar, Sentinel-1, Crosspolarization, SAR

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

SAR images acquired over the ocean, in particular over coastal areas, often show large variability of the backscattered radar power or the normalized radar cross section (NRCS). This variability can have many reasons: It can be due, among others, to variable sea surface winds, variable the air-sea interface stability, variable ocean surface currents, slick coverage, or rain.

Identifying and explaining rain signatures in C-band synthetic aperture radar (SAR) imagery acquired over the ocean is a challenging task, since rain can lead to an increase or a decrease of the NRCS [START_REF] Braun | Multifrequency scatterometer measurements on water surfaces agitated by artificial and natural rain[END_REF][START_REF] Alpers | Rain footprints on C-band synthetic aperture radar images of the ocean -Revisited[END_REF]. The reason is that several physical mechanisms contribute to the rain signature in SAR imagery over the ocean. They are: 1) scattering of the radar pulse from the sea surface, whose roughness is changed by raindrops impinging onto the sea surface, 2) scattering and attenuation from hydrometeors in the atmosphere (volume scattering and attenuation), and 3) scattering from the sea surface whose roughness is modified by rain-related winds, like downdraft or updraft. The first scattering mechanism process is very intriguing since raindrops impinging onto the sea surface can increase the NRCS due to scattering from rain-generated ring waves and from splash products, like stalks, craters, and raindrops emitted from the sea surface. On the other hand, they also can decrease the NRCS roughness due to generation of turbulence, which attenuates the short-scale waves (Bragg waves) responsible for the radar backscattering [START_REF] Bliven | An analysis of scatterometer returns from a water surface agitated by artificial: evidence that ring-waves are the main feature[END_REF][START_REF] Bliven | Rain generated ring-waves: measurements and modelling for remote sensing[END_REF][START_REF] Contreras | Surface effect of rain on microwave backscatter from the ocean: Measurements and modeling[END_REF][START_REF] Alpers | Rain footprints on C-band synthetic aperture radar images of the ocean -Revisited[END_REF][START_REF] Zhang | Rain effects on the hurricane observations over the ocean by C-band synthetic aperture radar[END_REF].

One outstanding phenomenon often observed in Sentinel-1 SAR images acquired over the ocean are areas of very high radar backscatter, often referred to as bright patches, which are observed when rain is present. The scattering mechanism causing the bright patches has been a matter of debate in the ocean remote sensing community since long time and has not been solved yet. In the past, several scattering mechanisms have been proposed to explain this phenomenon: 1) Scattering at splash products generated by raindrops impinging onto the sea surface (Atlas, 1994a(Atlas, , 1994b)), 2) radar pulse reflection from the sea surface followed by scattering at raindrops above the sea surface [START_REF] Jameson | SIR-C/X-SAR observations of rain storms[END_REF], 3) scattering at low-salinity "puddles" in the upper ocean layer generated by intense rainfall [START_REF] Wijesekera | Surface layer response to weak winds, westerly bursts, and rain squalls in the western Pacific Warm Pool[END_REF], and 4) scattering at steep slopes on the rim of the craters produced by impinging raindrops [START_REF] Braun | Untersuchungen zur Radar-Rückstreuung und Wellendämpfung beregneter Wasserober flächen[END_REF].

In this paper, we compare Sentinel-1 co-and cross-polarization SAR data with data from radar backscattering measurements carried out by ground-based and airborne radars, and we refer to the theory of radar backscattering from the melting layer as developed by d' Amico et al. (1998). Since the bright patches in the SAR images show similar characteristics as radar backscattering from the melting layer carried out with ground-based and airborne radars, we conclude that the underlying scattering mechanism is the same, i.e., volume scattering from wobbling, non-spherical, oblate hydrometeors within the melting layer.

The results obtained in this paper are of relevance for ocean surface wind retrieval using Cband SAR data, since scattering at hydrometeors in the melting layer can cause errors in wind retrieval. When inverting NRCS values in SAR images into wind speed, it is assumed that the backscattered radar power is solely due to surface scattering and receives no contribution from volume scattering. Thus, areas with simultaneously high co-and cross-polarization backscatter values due to volume scattering from the melting layer have to be flagged as areas of corrupted winds. We expect that wind retrieval in tropical storms is especially prone to be affected by volume scattering from the melting layer.

The rest of the paper is organized as follows: In Section 2, we review theories and experimental data on radar backsattering from the melting layer. . In Section 3, we present 4 representative examples of co-and cross-polarization Sentinel-1 SAR images on which pronounced radar signatures of rain cells are visible and relate them to quasi-simultaneously acquired weather radar images from the Hong Kong Observatory (HKO). In Section 4, we first present a SIR-C/X-SAR image showing the diversity of radar signatures of rain cells and then show how rain-induced surface scattering can interfere with volume scattering from the melting layer. In Section 5, we interpret the data in terms of a theory.on radar backscattering from hydrometores in the melting layer, and in Section 6, we summarize the results and draw conclusions.

Radar backscattereing from the melting layer

Basics

The melting layer is the layer in the atmosphere, where the irregularly falling ice particles undergo a phase transition from solid to liquid and where the ice particles are coated with liquid water due to meltin [START_REF] Co Szyrmer | Modeling of the melting layer. Part I: Dynamics and microphysics[END_REF]. Here, the radar backscatter or reflectivity (Z) is strongly enhanced, which is the result of complex interactions of dynamics and microphysics as described, e.g., by D' Amico et al. (1998) and [START_REF] Co Szyrmer | Modeling of the melting layer. Part I: Dynamics and microphysics[END_REF].

These authors have shown that, due to the horizontal gradients of the buoyancy in the melting layer, the hydrometeors experience tilting (canting) in the melting layer. Scattering of the radar pulse at randomly oriented (wobbling) melting particles causes the large radar signatures at both, co-and cross-polarization. These theoretical results have been validated by radar backscatter measurements carried out with ground-based radars. Although most ground-based radar backscatter mmeasurements aimed at investigating the melting layer have been carried out at co-polarization (see, e.g., [START_REF] Brandes | Freezing-level estimation with polarimetric radar[END_REF][START_REF] Boodoo | Application of dual-polarization radar meltinglayer detection algorithm[END_REF][START_REF] Kumjian | Principles and applications of dual-polarization weather radar. Part I: Description of the polarimetric radar variables[END_REF], some measurements were also carried out at cross-polarization.

As early as 1952, [START_REF] Browne | Cross-polarization of the radar melting-band[END_REF] performed cross-polarization measurement with a ground-based radar having a wavelength of 3.2 cm. They found that, at crosspolarization, the backscattered radar power from the "melting layer" was larger than the one from raindrops below and snowflakes above the freezing level. Sometimes, they could detect the melting layer only at cross-polarization, but not at co-polarization.

When dealing with detection of melting layers by radars, the following parameters are commonly used by radar meteorologists to determine geometrical and dynamical properties of the hydrometeors.: 1) horizontal reflectivity (ZHH), 2) vertical reflectivity (ZVV), 3) crosspolarization reflectivity (ZHV and ZVH), 4) differential reflectivity (ZDR), 5) linear depolarization ratio (LDR), and 6) cross-correlation coefficient (ρVH and ρHV). ZDR, ρVH, and LDR are defined as follows:

ZDR =10 log10 (ZHH/ZVV) (1) ρVH = < ZVVZHH >/ (< ZVV 2 >< <ZHH 2 >) 1/2 (2)
LDR =10 log10 (ZHV/ZHH) or 10 log10 (ZVH/ ZVV).

(3)

The subscipts attached to Z denote the polarization; the first subscript denotes the polarization of the transmitted signal, and the second one the polarization of the received signal, V = vertical polarization, H = horizontal polarization.

Raindrops have always an oblate shape, i.e., their width is larger than their height and they have similar shapes as pebble stones found on beaches. The differential reflectivity ZDR is primarily an indicator of the shape and size of the hydrometeors. The values of ZDR are positive for scattering from raindrops (typically 1.1 to 1.2, but can be as large as 5.0 for large raindrops).

For scattering from hydrometeors in the melting laye.r, they are also positive, but usually slightly larger than for raindrops, and for scattering from ice hydrometeors they are slightly positive as well as slightly negative [START_REF] Jameson | The interpretation and meteorological application of radar backscatter amplitude ratios at linear polarizations[END_REF]. Thus, sometimes the melting layer is not detectable on ZDR images.

On the other hand, the correlation coefficient (ρVH) is a measure for the similarity of the radar backscatter in HH and VV polarizations and thus a measure of how uniformly the scatterers are distributed. As shown by [START_REF] Ryzhkov | Interpretation of polarimetric radar covariance matrix for meteorological scatterers: Theoretical analysis[END_REF], ρVH is a function of the mean canting angle of the hydrometeors. Its value is normally close to 1.0 for raindrops and snowflakes, but below 0.95 for hydrometeors in the melting layer.

As an example of such measurements, we show in Fig. The height H of the melting layer can be calculated from the position of the melting layer band visible in the ZDR or the ρVH image by applying the relationship

H = R tan θ + H0, (4)
where R denotes the distance of the circular melting layer band from the position of the radar, θ the elevation angle of the radar beam, and H0 the height of the weather station above mean sea level. Applied the radar images depicted in Fig. 8 and using the values R = 90 km, θ = 2.7 o , and H0 = 500 m in Eq. 4, we obtain for the height of the melting layer the value H = 4.7 km. This lies in the expected range of melting layer heights in this region during strong convective rain events. We have inserted this melting layer height as a red line in Fig. 8, Panel e. LDR is a measure of how much larger the cross-polarization (HV or VH) scattering is compared to the a) b)

ZDR co-polarized (HH or VV) scattering. This parameter is not available from the Hong Kong weather radar. At present, operational ground-based radars for weather monitoring provide measurements only at HH and VV polarizations.

Theories of radar backscattering from the melting layer

Since the diameter of hydrometeors is small (typically < 0.3 cm) compared to the wavelength of the Sentinel-1 SAR (5.4 cm), the Rayleigh scattering theory can be applied [START_REF] Oguchi | Electromagnetic wave propagation and scattering in rain and other hydrometeors[END_REF][START_REF] D'amico | An anisotropic model of the melting layer[END_REF]. In early scattering models [START_REF] Dissanayake | Radar and attenuation properties of rain and bright band[END_REF]McEvans, 1989: Willis and[START_REF] Willis | Structure of the melting layer in mesoscale convective system stratiform precipitation[END_REF], it was assumed that the hydrometeors have a spherical shape. In this case, a horizontally (vertically) transmitted electromagnetic wave generates a purely horizontally (vertically) backscattered electromagnetic wave. Due to the symmetry of the target, all vertical (horizontal) components of the electric field oscillations cancel and ZDR is zero. The increase of backscattered radar power from the melting layer, which is then due purely to the increase of the dielectric constant.

However, the hydrometeors usually have an oblate or a flattened spheroidal shape. They become oblate due to air resistance when falling downwards. In this case, the backscattered horizontally and vertically polarized radar signals are different and consequently, ZDR is nonzero. When the long axis of the hydrometeor is aligned horizontally, then the electric field oscillations parallel to the long axis dominate and the horizontal backscattered radar signal becomes larger than the vertical one, i.e., ZDR becomes positive [START_REF] D'amico | An anisotropic model of the melting layer[END_REF]. The larger the raindrops, the more oblate they are, and the larger is ZDR. However, a crosspolarized backscattered radar signal can only occur, when the oblate hydrometeors are tilted out of plane of incidence of the transmitted electromagnetic wave (Fig. 2). The tilting of the elongated hydrometeors is particularly strong in the melting layer due to large horizontal gradients of the buoyancy encountered there. Thus, it is tilting, which gives rise to the large cross-polarization radar signature of rain in the melting layer. If there were no tiling or wobbling of the spheroidal hydrometeors, then there would be no cross-polarization signature due to scattering. (However, there might be a small cross-polarization signature due to differential wave attenuation between horizontal and vertical polarizations). Since the cross-polarization NRCS is always smaller than the co-polarization one, LDR is always negative. Typical LDR values for raindrops and ice particles are -25 to -30 dB and for melting hydrometeors -10 to -20 dB [START_REF] Houze | Cloud Dynamics[END_REF].

Ground-based and airborne radar measurements of the melting layer

In this sub-section, we present data of co-and cross-polarization measurements carried out with a ground-based radar and a nadir looking airborne radar, which confirm the theory presented in the previous section. They show that the melting layer manifests itself not only in an increase of the co-polarized reflectivity (in these examples of ZHH), but also in an increase of the cross-polarization reflectivity. The ground-based measurements were carried out by the fully-polarimetric C-band radar of the German Aerospace Center (DLR) at Oberpfaffenhofen, Germany [START_REF] D'amico | An anisotropic model of the melting layer[END_REF], and the airborne measurements by the NASA/JPL ABMAR airborne RADAR [START_REF] Jameson | SIR-C/X-SAR observations of rain storms[END_REF]. Radar meteorologists usually do not plot the crosspolarization reflectivity, but the linear depolarization ratio (LDR) defined by Eq. 3, which is the ratio of the reflectivity at cros-polarization to the one at co-polarization or in logarithmic scale, the difference between 10 log ZHV and 10 log ZHH (or 10 log ZVH and 10 log ZVV), see, e.g., D' Amico et al.,1998. Figs. 3 and 4 show that the reflectivity Z and LDR are strongly enhanced in the melting layer where LDR takes the values are -18 dB (Fig. 3) and -11 dB (Fig. 4). Note that in both measurements, the peaks in Z and LDR are not collocated, which is a commonly observed phenomenon, but its explanation is beyond the scope of this paper. The LDR of the melting layer is -11 dB. Adapted from [START_REF] Jameson | SIR-C/X-SAR observations of rain storms[END_REF].

Concurrent Sentinel-1 SAR and weather radar measurements of rain cells

In (2015) for cross-polarization.

The 6 June 2018 rain event

The SAR images of this event (Fig. 5

The 12 April 2017 rain event

The SAR images of this event (Fig. 6, Panels a and b) show a strong rain band extending in the NE-SW direction, which gives rise to pronounced peaks in the VV NRCS and VH NRCS transects. The VV NRCS and the VH NRCS values of the highest peaks are -10 dB and -22 dB, respectively, such that the ratio of VH NRCS to VV NRCS is, in logarithmic scale, -12 dB.

However, here the heights of the peaks relative to the background are lower than in the previous event caused by higher ambient winds. A noteworthy feature visible in Panel d is the pronounced dip in the VV NRCS profile to the right (east) of the peak. We interpret it as caused by the decrease of surface scattering due to damping of the short waves (Bragg waves) due to rain-generated turbulence.

The 3 September 2017 rain event

The SAR images of this event (Fig. 7 ECMWF wind ratio of VH NRCS to VV NRCS is, in logarithmic scale, -13 dB. For this event, we were able to retrieve the height of the melting layer from the HH and VV polarization data of the Hong Kong weather radar, see Fig. 1. This height is 4700 m, while the zero-degree height level measured by the radiosonde is 5200 m.

The height difference is due to the fact that melting starts at the zero-degree Celsisus height level, but reaches its maximum at a lower height.

Interplay between volume and surface scattering

As stated before, the backscattered radar signal of Sentinel-1 SAR receives contribution from volume scattering as well as from surface scattering. This can give rise to significant distortions of the theoretically expected radar signature of a rain cells due to volume scattering as evident in the Sentinel-1 SAR image depicted in Fig. 6. To further illustrate this complexity, we show in Fig. 9 simultaneously acquired multi-polarization/multi-polarization SAR images with rain cells from the SIR-C/X-SAR mission in 1994 [START_REF] Melsheimer | Investigation of multifrequency/multipolarization radar signatures of rain cells over the ocean using SIR-CX-SAR data[END_REF]. These images clearly show the strong variability of radar signatures of rain cells with radar frequency and polarization as stated in the Introduction. For L-band, the radar signal penetrates the melting layer almost unattenuated and therefore, the L-band radar signature of the rain cell has its origin almost exclusively in surface scattering. The black patch in the center of the L-band images is due to impinging raindrops, which generate turbulence in the upper ocean layer and attenuate there the L-band Bragg waves. When comparing the C-band VV and HH images with the HV image, we see that rain cells have a higher contrast in the VH image than in the VV and HH images.

The reason is that the background C-band NRCS due to wind-induced surface scattering is in the VH image much smaller than in the VVand HH images. Furthermore, we also note that for C-band, the size of the bright patch is smaller in the VV and HH images than in the HV image and that there is a small dark patch adjacent to the bright one in the HH/VV images. We interpret it as caused by the superposition of volume scattering and surface scattering: In the bright area, the increase of the NRCS induced by volume scattering at hydrometeors in the melting layer is dominant, and in the dark area, the attenuation of the Bragg waves is dominant due to raingenerated turbulence, thereby inducing the reduction of the NRCS.

Another remarkable feature visible in Fig. 9 is the small bright band located west of the large bright patch in the center, which is visible only in the CHV image, but not in the XVV image nor in the CVV and CHH images. We interpret this feature as a radar signature of a small rain band in the melting layer, which is less strong than the signature of the large rain cell located in the center of the CHV image (see Fig. 9 f). This signature is not visible in the CHH and CVV images, because its NRCS value is much smaller than the NRCS of the background due to high winds. This suggests that cross-polarization is better suited for detecting scattering from the melting layer than co-polarization, because it receives less contribution from surface scattering, and therefore the signal-to-noise ratio is higher. Note, that this small rain band is also visible in the L-band VV /HH images as a faint dark band, which has its origin in the damping of Bragg waves by the rain-generated turbulence.

Interpretation of the data

The analysis of the quasi-simultaneously acquired Sentinel-1 SAR images and weather radar data presented in Section 3 show that the bright patches visible in the Sentinel-1 SAR images are related to rain cells (Fig. 5), rain bands Figs. 6 and7), or larger rain areas (Fig. 8). Note that the weather radar images depicted in the Panel c of Figs. 5 -8 show the radar reflectivity at a height of 3 km, while the melting layer lies above this height. The height of the melting layer can be estimated from the height of the zero-degree level, which usually is located few hundred meters higher than the the center of the melting layer. In one case, in the 25 June 2019 event (Section 3.4), we were able to deternine the height of the melting layer from weather radar data of the HKO (see ZHV and ρHV images depicted in Fig. 1). Panels e of Figs. 5678show that during all four rain events hydrometeors were present in the melting layer. Panels d of Figs. 5 -8 show that the peaks in the VV NRCS and VH NRCS profiles are highly correlated. In addition, peaks in the SAR NRCS profiles are also correlated reasonably well with peaks in the reflectivity profiles obtained from CAPPI weather radar data. (Note, that the radar signatures of the rain cells in the Sentinel-1 SAR images and in the CAPPI displays refer to different heights, i.e., the melting layer height and the 3 km height, respectively). Furthermore, the height of the peaks relative to the background is always higher at VH polarization than at VV polarization. The VV NRCS peak values range from -9 dB to -10 dB, and the VH NRCS peak values from -20 dB to -24 dB.

In order to prove that the bright patches visible in Sentinel-1 SAR images are due to volume scattering from hydrometeors in the melting layer, we compare the SAR data with data from radar backscattering measurements from the melting layer carried out with ground-based and airborne radars. For this comparison, we employ the parameters LDR and LDRSAR, where LDR denotes the Depolarization Ratio, defined by Eq. 3, which is used by radar meteorologists to localize the melting layer, and LDRSAR, which we define by: LDRSAR =10log 10 (σ total VH /σ total VV ) or LDRSAR =10log 10 (σ total HV /σ total H )

(5)

Here σ total denotes the sum of the NRCS due to volume scattering, σ volume , and the NRCS due to surface scattering, σ surface :

σ total = σ surface + σ volume (6)
The subscripts denote the polarizations for transmission in analogy to the definition of LDR (Eq. 3). When calculating σ total from the backscattering data, one has to correct for the attenuation of the radar pulse when it propagates through the melting layer and through the layer above containing frozen hydrometeors ( snowflakes) and the layer below containing liquid hydrometeors (rain drops)y. However, while for X-band, attenuation is a major effect (see, e.

g., [START_REF] Danklmayer | Assessment of atmospheric propagation effects in SAR images[END_REF], it is a minor effect for C-band as long as the rain rate is not too high (< 20 mm/h) (Tounadre and Morland, 1997). E.g., when using the empirical aR b relation of Olsen et al. (1998) for calculating the attenuation, one obtains for a 5 km thick rain cell with rain rate of 30 mm/h an attenuation of 1.3 dB. In this investigation, we neglect attenuation, since it has no effect on our result on the nature of the scattering mechanism.

In all Sentinel-1 SAR images presented in Figs, 5 -8, volume scattering is the dominant scattering mechanism in the bright patches. The contribution of surface scattering to the total scattering is particularly small in the events of 6 June 2018 (Fig. 5) and of 25 June 2019 (Fig. values in radar backscattering measurements from the melting layer, must be the same, i.e., the scattering mechanism is volume scattering from wobbling, non-spherical, oblate hydrometeors within the melting layer.

Summary and conclusions

In this paper, we have provided evidence that the scattering mechanism causing the large radar backscatter values at co-and cross-polarization often observed in Sentinel-1 SAR images of the sea surface in the presence of rain ("bright patches") is scattering from wobbling, nonspherical, oblate hydrometeors within the melting layer. In the past, several other scattering mechanism have been proposed, which are all based on the modification of the sea surface roughness due to rain drops impinging onto the sea as described in the Introduction. In a recent review paper dealing with rain footprints on C-band SAR images of the ocean [START_REF] Alpers | Rain footprints on C-band synthetic aperture radar images of the ocean -Revisited[END_REF], it was stated in the Conclusion section: "The scattering mechanism causing the bright patches in C-band, co-polarized SAR images of rain cells could not be determined". In this paper, we have identified the scattering mechanism by comparing Sentinel-1 co-and cross-by radar scattering from hydrometeors in the melting layer: VV NRCS must have values around -10 dB or larger and VH NRCS values around -20 dB, such that LDRSAR has values between -15 dB and -10 dB. We expect this criterion to be applicable for incidence angles between 30 and 70 degrees and for wind speeds below 20 m/s. The identication of areas in C-band SAR affected by scattering from the melting layer is of relevance for ocean surface wind retrieval using C-band SARs, since in these areas the conventional wind retrieval algorithm is not applicable.

  1 a ZDR image (Panel a) and a ρVH image (Panel b), which correspond to the weather radar image depicted in Fig. 8. The melting layer can be identified on the ZDR image (Panel a) as a partially circular band with a ZDR value slightly larger than 1.0 (marked by a black arrow) surrounded by areas with ZDR = 1.0 representing rain in the inner section of the circular band and snow in the outer section. The melting layer can also be identified on the ρVH image (Panel b) as a partially circular band of slightly decreased ρVH values (0.92 -0.94). Note, that in this case the melting layer is best visible in the ρVH image.

Fig. 1 .

 1 Fig. 1. a) PPI image of the differential reflectivity (ZDR) and (b) the correlation coefficient (ρHV) acquired by the weather radar of the HKO on 25 June 2019 at 10:25:58 UTC (18:25:58 HKT) quasi-simultaneously with the Sentinel-1A data acquisition at 10:25:07 UTC (Fig. 8). The arrows mark the melting layer as a circular band of slightly increased ZDR values and of slightly decreased ρVH values (0.92 -0.94),

Fig. 2 .

 2 Fig. 2. Orientation of a spheroidal hydrometeor defined by the direction of its axis of symmetry, u. The angles θ and Φ denote the canting angle and the azimuth angle, respectively. Reproduced from D'Amico et al. (1998).

D

  'Amico et al. (1998) have carried out detailed modeling of cross-polarization radar backscattering from the melting layer. In their model, the hydrometeors in the melting layer consist of ice and liquid water with varying concentrations depending on their height in the melting layer. The oblate spheroidal hydrometeors are subject to tilting with respect to the horizontal. Using this model, they were able to simulate quite well the profiles shown in Fig.3.

Fig. 3 .

 3 Fig. 3. Vertical profiles (in meters) of ZHH and LDR measured on 17 August 1994 at 1433 LT by the DLR ground-based radar showing the melting layer in both parameters. The LDR of the melting layer is -18 dB. Note that the peaks of ZHH and LDR are not collocated. Reproduced from D'Amico et al. (1998).

Fig. 4 .

 4 Fig. 4. Vertical profiles of ZHH and LDR measured by the NASA/JPL ABMAR airborne radar that flew over a melting layer associated with a forming tropical cyclone in the Tropical Western Pacific in 1993.

  this section, we present four representative examples of Sentinel-1 SAR images showing areas of strong C-band radar backscattering ("bright patches") at co-and crosspolarization and compare them with quasi-simultaneously acquired weather radar data of the HKO. The SAR onboard the Sentinel-1 satellites operates at C-band (5.4 GHz) and has different exclusive acquisition modes. Here we use only SAR images captured in the Interferometric Wide (IW) swath mode at VV and VH polarizations (spatial resolution: 20 m x 22 m, swath width: 250 km, incidence angle range: 29.1 o -46.0 o ). All images shown in this section are Level 1 Ground Range Detected (GRD) products provided by ESA via the Copernicus O)pen Access Hub (https://scihub.copernicus.eu). They were acquired during ascending satellite passes with the SAR antenna pointing to the right of the satellite track. The Hong Kong (HK) weather radar is a dual-polarization C-band radar with a half-power beam width of 0.5 degrees and has two modes of operation. One is volume scan with scanning at different elevation angles (the lowest one being θ = 2.7 º) with a range of 256 km and a repeat cycle of 6 minutes. The other one is horizontal scan with a range of 512 km and a repeat cycle of 12 minutes. The 64 km and 128 km images are derived from the 256 km range data. The weather radar takes measurements at HH and VV polarizations, but not at cross-polarization (VH or HV). Among other products, it provides reflectivity measurements converted into rain rate at a height of 3 km above mean sea level in the Constant Altitude Plan Position Indicator (CAPPI) display. Furthermore, also Plan Position Indicator (PPI) maps are available in real time, from which vertical cross sections of the reflectivity can be generated on demand. The PPI maps shown Panels e are from the 256 km range data. For all four rain events presented in this section (Fig. 5 -Fig. 8), we show in Panels a and b VV-and VH-polarization images captured simultaneously by the Sentinel-1 SAR, in Panel c the reflectivity image of the HK weather radar, in Panel d the profiles of the VV NRCS, VH NRCS values and radar reflectivity values along the transects inserted in Panels a, b, and c , in Panel e the reflectivity in the range-height indicator (RHI) presentation along transects shown in Panel c together with the 3 km height line (solid black line) and the zero-degree Celsius height line obtained from radiosonde data (dashed black line), and in Panel f the wind field as provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). Inserted in Panel d are the VV NRCS and VH NRCS values calculated from the background wind adjacent to the bright patches by using the C-band Geophysical Model function CMOD5n developed by Hersbach (2015) for VV polarization and the one developed by Hwang et al,

Fig. 5 .

 5 Fig. 5. Sentinel-1 A SAR images acquired on 6 June 2018 at 10:25 UTC at VV (Panel a) and at VH (Panel b) over the South China Sea off the Hong Kong coast showing radar signatures of a rain band over the sea. The inserted red lines denote the transects along which the variation of NRCS values was measured. The thick red arrow marks the direction into which the SAR antenna is pointing. Panel c shows the CAPPI reflectivity image at a height of 3 km acquired by the HKO weather radar at 10:24 UTC (18:24 Hong Kong

Fig. 6 .

 6 Fig. 6. Same as in Fig. 5, but for 12 April 2017 at 10:25 UTC (Sentinel-1) and 10:24 UTC (weather radar).

  , Panels a and b) show a strong rain band aligned in the NW-SE direction, which gives rise to the peaks in the VV and VH NRCS profiles. The VV NRCS and the VH NRCS values of the highest peaks are -10 dB and -21 dB , respectively, such that the ratio of VH NRCS to VV NRCS is, in logarithmic scale, -11 dB. Note that the weather radar image (Fig. 7, Panel c) shows also rain to the west of this rain band, but the VV and VH SAR images do not show these strong radar signatures. Here the radar signature of rain is solely due to scattering from the sea surface whose roughness is modified by raindrops impinging onto the sea surface.

Fig. 7 .

 7 Fig. 7. Same as in Fig. 1, but for 3 September 2017 at 10:25 UTC (Sentinel-1) and 10:24 UTC (weather

Fig. 8 .

 8 Fig. 8. Same as in Fig. 1, but for 25 June 2019 at 10:25 UTC (Sentinel-1) and 10:24 UTC (weather radar). The red line inserted in Panel e denotes the height of the melting layer as inferred from the ZHV and ρHV

Fig. 9 .

 9 Fig. 9. SIR-C/X-SAR images acquired at L-, C-, and X-band simultaneously over the Gulf of Mexico at 08:11 UTC on 18 April 1994, displaying the strong dependence of the radar signature on radar frequency and polarization (reproduced from Alpers et al., 2016). Note, in particular, the difference in the rain signature patterns of the rain cell in the CVV/CHH and CVH images, which are due to superposition of volume and surface scattering in the CVV and CHH images.

Fig. 10 a

 10 Fig. 10 a shows a zoom of the CVH image depicted in Fig. 9 f with a transect line inserted

Fig. 10 .

 10 Fig. 10. a) Section of the CHV SAR image shown in Fig. 9 f; b) Ratio of the NRCS values at HH and HV polarizations (solid line), and at VV and VH polarizations (dashed line). Reproduced from Braun (2003).

  8): -20 dB versus -10 dB and -23 dB versus -10 dB. Thus, in these cases, the contribution from surface scattering can be neglected in LDRSA R such that it contains, like LDR, only contributions from volume scattering, We now compare LDRSAR values measured in then bright patches with LDR values measured in radar backscattering from the melting layer by using ground-based and airborne radars. In these cases, LDRSAR has peak values in the bright patches of -12 db and -13 dB, which compare well with LDR values characterizing the melting layer (see Section 2,3). Furthermore, the SIR-C/SAR data presented in Section 4 (Fig.10 b) show LDRSAR values of -15 dB and -16 dB for horizontally and vertically polarized transmitted signals, respectively, which also lie in the range of LDR values characterizing radar backscatter from the melting layer as measured by ground-based and airborne radars. Thus. we conclude that the scattering mechanisms causing the bright patches in the SAR images and the high LDR

Acknowledgments

This study was supported by ESA by the project: "Sentinel-1A Mission Performance Center for access of Sentinel-1 Interferometric Wide Swath data" and by the China Scholarship Council (CSC), who awarded a Ph. D. scholarship to Yuan Zhao. We are grateful to the colleagues in LOPS of IFREMER for their help in data acquisition and data processing and to Nicole Braun for providing Fig. 10. Last not least, we want to thank the three anonymous reviewers for their very valuable suggestions to improve the manuscript.

polarization SAR data with data obtained from radar backscattering measurements from the melting layer carried out with ground-based and airborne radars. We have calculated the ratio of the NRCS values at co-and cross-polarization in the bright patches in the SAR images (LDRSAR) and compared them with the ratio of the radar reflectivity at cross-and co-polarization (LDR values) measured in radar backscattering experiments from the melting layer using ground-based and airborne radars. Since both ratios have similar values and since it is known that the last one (LDR) characterizes the melting layer, we conclude that the bright patches are caused by scattering from the melting layer. In this context, we would like to mention that such scattering mechanism was already suspected by [START_REF] Katsaros | Wind Fields from SAR: Could They Improve Our Understanding of Storm Dynamics?[END_REF] who analyzed copolarized Radarsat-1 C-band data of the hurricane Danielle in 1998 and noticed small white spots within a rain cell, which they then compared them passive microwave data from the Special Sensor Microwave/Imager (SSM/I) onboard DMSP satellites. They suspected that the white spots result from "randomnly oriented ice particles in th cloud". We know of no paper in which this idea, i.e., the comparison of C-band SAR with passive microwave data, was pursued further.

When screening C-band co-polarization SAR images for radar signatures caused by scattering from the melting layer, one must be aware of the fact that they are sometimes distorted by overlapping radar signatures due to rain-induced surface scattering, as evident in the SIR-C/X-SAR CHH and CVV images (compared to the CHV image) shown in Fig. 9. This suggests that cross-polarization SAR images are better suited to detect radar signatures of the melting layer than co-polarization SAR images. Furthermore, we conclude from


Quilfen Yves
email: yquilfen@ifremer.fr

Chapron Bertrand

On denoising satellite altimeter measurements for highresolution geophysical signal analysis

Keywords: Altimeter Measurement Noise, Empirical Mode Decomposition, Mesoscale Variability

Satellite radar altimeter observations are key to advanced studies in ocean dynamics, particularly those focusing on mesoscale processes. To resolve scales below about 100 km, because altimeter measurements are often characterized by a low signal-to-noise ratio (SNR), low-pass filtering or leastsquares curve fitting is generally applied to smooth the data before analysis. Here, we present an alternative method. It is based on Empirical Mode Decomposition (EMD) developed to analyze nonstationary and non-linear processes, which adaptively projects a signal on a basis of empirical AM/FM functions called Intrinsic Modulation Functions (IMFs). Applied to a Gaussian noise signal, the EMD provides a set of IMFs with a predictable distribution of noise energy that can be exploited by waveletinspired threshold methods to provide an efficient approach to data denoising. The EMD method performs a local analysis of the SNR, does not require a priori assumptions about the underlying geophysical signal, e.g., its degree of smoothness or its compliance with a particular mathematical framework. The signal is simply assumed to be the sum of a piecewise-smooth deterministic part and a stochastic part. The proposed EMD-based denoising approach is therefore well suited for mapping non-linear features, such as strong gradients, and extreme values without significant smoothing. Using Jason-2, Cryosat-2, and Saral/AltiKa significant wave height measurements, the method provides an effective means of mapping overlooked geophysical variability of sea state at scales between about 100 km and 25 km, a range largely impacted by low SNR. Below 25 km, a spectral hump caused by inhomogeneities in the altimeter footprint essentially dominates the signal. In addition, the EMD method provides a consistent approach for longterm noise analysis and monitoring under global and local conditions. The proposed method is a step forward that will enable better exploitation of the unique set of altimeter observations that now covers more than 25 years.

algorithms applied to the analysis of altimeter signals. Efforts are being made to cross-calibrate measurements from the different satellite missions, in particular to improve their usefulness for climate studies, as in the European Space Agency's (ESA) Climate Change Initiative (CCI) project.

One of the main limitations is that altimeter measurements are often characterized by a low signal-tonoise ratio (SNR), which can have a significant impact on geophysical analysis at spatial scales below 100 km. The main sources of noise are the altimeter waveform re-tracking algorithms, instrumental noise and geophysical variability in the altimeter footprint [START_REF] Sandwell | Retracking ERS-1 altimeter waveforms for optimal gravity field recovery[END_REF][START_REF] Thibault | Relative Performance of the MLE3 and MLE4 Retracking Algorithms on Jason-2 Altimeter Waveforms[END_REF][START_REF] Dibarboure | Investigating short-wavelength correlated errors on low-resolution mode altimetry[END_REF]. In addition, because the recovered parameters are obtained from the same waveform re-tracking algorithm, they can have highly correlated errors, i.e., standard MLE4 processing produces four estimated parameters with correlated errors (SSH, Hs, sigma0, and off-nadir angle). Fine-scale ocean dynamics analysis therefore requires preliminary noise filtering, and low-pass or smoothing filters (e.g., Lanczos, running mean, or Loess filter) are frequently used. Each of these filters efficiently smooths noisy data but results in a systematic loss of small-scale ( < 100 km) geophysical information or the creation of artifacts in the geophysical variability analyzed. In addition, such filters only remove some of the high-frequency noise, which is also the result of approaches that infer a correction to eliminate correlated errors from other aspects of the waveform data (Zaron and de Carvalho, 2015;[START_REF] Quartly | Removal of covariant errors from altimetric wave height data[END_REF][START_REF] Tran | Assessing the effects of sea-state related errors on the precision of high-rate Jason-3 altimeter sea level data[END_REF], leaving a substantial amount of low-frequency noise in the data. Such low-frequency noise remaining in the signals after smoothing will still interfere with the precise measurement of geophysical peak height, width, and position [START_REF] O'haver | Pragmatic Introduction to Signal Processing[END_REF].

To overcome these difficulties, we present an alternative approach providing an adaptive noise elimination method for satellite altimeter measurements. It is based on the non-parametric Empirical Mode Decomposition (EMD) method developed to analyze non-stationary and non-linear signals [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF][START_REF] Wu | A study of the characteristics of white noise using the empirical mode decomposition method[END_REF][START_REF] Huang | A review on Hilbert-Huang transform: Method and its applications to geophysical studies[END_REF]. EMD is a scale decomposition into a limited number of amplitude and frequency modulated functions (AM/FM), among which the Gaussian noise distribution is predictable [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF]. It can therefore provide the basis for a noise elimination approach, optimized and applied to a variety of simulated deterministic signals, with results often superior to those of wavelet-based techniques (Kopsinis and McLaughin, 2009). For instance, to identify the occurrence of rogue waves, [START_REF] Hadjihosseini | Stochastic analysis of ocean wave states with and without rogue waves[END_REF] used EMD to filter high-frequency variability related to the background sea state in a 1-Hz surface elevation record taken by an ultrasonic wave gauge in the Sea of Japan. Recently, EMD analysis has been successfully applied to altimeter data to analyze wave-current interactions [START_REF] Quilfen | Storm waves focusing and steepening in the Agulhas current: satellite observations and modeling[END_REF][START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF], known to predominate at scales below 100 km [START_REF] Ardhuin | Small-scale open-ocean currents have large effects on wind-wave heights[END_REF]. These studies showed how EMD analysis is suitable for processing non-stationary and non-linear signals in low SNR conditions, enabling a better exploitation of current and past satellite altimetry missions. More specifically, it can help to consistently recover power spectra of Hs and SSH variability, as well as to provide an accurate representation of non-linear features, such as strong gradients, and extreme values [START_REF] Quilfen | Storm waves focusing and steepening in the Agulhas current: satellite observations and modeling[END_REF][START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF]. As the EMD denoising method was only briefly presented in these two previous articles, it is the purpose of the present study to present it in detail, especially since it is being used for processing of the upcoming sea state CCI and Copernicus Marine Environment Monitoring Service (CMEMS) products.

Particular attention is paid to comparing the results obtained from EMD with those obtained using a smoothing technique based on a Lanczos low-pass filter. The latter method is very popular among oceanographers, due to its simplicity and favorable characteristics. It is widely used for altimeter noise filtering, as performed in the current altimeter sea level and sea state products distributed by CMEMS [START_REF] Taburet | Sea Level TAC-DUACS Quality Information Document[END_REF][START_REF] Taburet | WAVE Thematic Assembly Centre Quality Information Document[END_REF]. A Lanczos filter is a reformulated low-pass filter that incorporates a multiplicative factor to attenuate spectral ringing and signal energy leakage in the pass band (caused by overshoot ripples known as the Gibbs' phenomenon, [START_REF] Emery | Data analysis methods in physical oceanography[END_REF]. In practice, the signal is reconstructed by adding the low-frequency Fourier components to those within a frequency window, centered on the cut-off frequency, which are weighted. A side effect of the Lanczos filter is excessive signal smoothing in the presence of strong gradients and extreme values. Designing a practical Lanczos filter is thus not easy, because it requires setting of both the cutoff wavelength and the number of weighted Fourier terms required to achieve the desired roll-off between the stop and pass bands. It always becomes a compromise. Signals filtered with a Lanczos filter are therefore strongly constrained by the inherent Fourier mathematical framework decomposition, while EMD filtered signals are based on a different assumption. The EMD decomposition assumes that the observed variable is the sum of a piecewise-smooth deterministic background part and a stochastic part.

Hereafter, for demonstration purposes and because it is implemented to process data distributed to users, this study focuses on the EMD-based method for processing Hs measurements. However, it also applies to SSH and sigma0 measurements. Section 2 presents the data used, section 3 the principles of EMD, its algorithm and properties, and section 4 describes the EMD-based denoising algorithm applied to altimeter measurements. In section 5, the results are presented for a few illustrative cases and the variances of both the resulting filtered signals and associated noise are analyzed. In addition, the uncertainty parameter associated with the denoised Hs is presented and characterized. Finally, section 6 provides a summary.

-Data

The multi-mission altimeter dataset processed at IFREMER by [START_REF] Queffeulou | Global altimeter Hs data set, version 11[END_REF] was analyzed from January 2014 to December 2016. This dataset contains inter-calibrated Hs and sigma0 measurements for all altimeter missions from the ERS-1 mission starting in 1991 and lasting until April 2017 for the Saral/AltiKa, Cryosat-2 and Jason-2 missions. The present study uses measurements from these last three altimeters. Each altimeter performed measurements at nadir along the satellite track, and standard processing delivers 1 Hz data with ground sampling that varied slightly from 6 to 7 km depending on the altimeter. The study uses quality-controlled, unfiltered and nonresampled, along-track Hs measurements. Data editing performed on the [START_REF] Queffeulou | Global altimeter Hs data set, version 11[END_REF] dataset was based on testing the quality flags of the AVISO Geophysical Data Record products, which are the primary source, and on along-track data screening to eliminate anomalous Hs values that deviate too far from adjacent values. However, this did not eliminate all data contaminated by land or sea ice transitions, rain, slicks or very light winds (sigma0 blooms). Outliers associated with such conditions are therefore likely to have remained in the analyzed dataset, causing the so-called spectral hump, which adds a spatially coherent error to the Gaussian white noise, especially at scales below about 25 km [START_REF] Dufau | Mesoscale resolution capability of altimetry: Present and future[END_REF][START_REF] Dibarboure | Investigating short-wavelength correlated errors on low-resolution mode altimetry[END_REF]. Improved data editing will be a necessary to mitigate these errors.

-Empirical Mode Decomposition principles and algorithm

EMD is a method to decompose signals into a limited number of scale-dependent components, called intrinsic modulation functions (IMFs), each modulated in amplitude and frequency. An IMF represents a simple oscillatory mode as a counterpart to a simple harmonic function derived from Fourier analysis. However, instead of imposing a constant amplitude and frequency in a simple harmonic component, an IMF can have an evolving amplitude and frequency. The decomposition is based on the local characteristic sampling scale of the data and is therefore applicable to non-linear and non-stationary processes. IMFs have well-behaved Hilbert transforms, from which instantaneous frequencies can be calculated. Thus, any event can be well localized in both the time and frequency domains [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF].

Unlike a more conventional decomposition method such as the wavelet technique, which projects the signal considered onto a basis of predefined wavelet functions, IMF expansion depends only on the analyzed signal. IMF calculation is performed using an iterative procedure called sifting. By construction, IMFs share the following properties: they are zero mean, all their maxima and minima are respectively positive and negative, they have the same number (or + / -1) of zero-crossings and local extrema, and they are modulated in both amplitude and frequency (AM/FM). EMD belongs to the class of multi-scale, orthogonal and adaptive decompositions. The EMD method ranks IMFs by construction from the first, containing the finest scales, to the last, containing the global trend.

EMD adaptively decomposes a signal x(t) into a small number L of IMFs so that :

ℎ (𝑖) (𝑡), 1⩽𝑖⩽𝐿, (1)
𝑥(𝑡) = ∑ 𝐿 𝑖 = 1 ℎ (𝑖) (𝑡) +𝑑(𝑡)
where d(t) is a residue that is a non-zero mean slowly varying function with only a few extrema.

The IMF number, L, depends on the length of the record and typically varies from 1 to 10 for the lengths analyzed in the altimeter dataset.

The IMFs are calculated successively, the first containing the shortest scales by the construction of the algorithm. Each IMF is estimated using an iterative process called sifting (illustrated in Figure 1), which is successively applied to the residual multicomponent signal with :

𝑥 (𝑖) (𝑡)
to derive the first IMF

𝑥 (𝑖) (𝑡) = 𝑥(𝑡),𝑖 = 1
to derive the IMFs of rank > 1

(2)

𝑥 (𝑖) (𝑡) = 𝑥(𝑡) -∑ 𝑖 -1 𝑗 = 1 ℎ (𝑗) (𝑡),𝑖 ≥ 2
For the estimation of an IMF of rank i , during the (n + 1) th sifting iteration, a temporary IMF ℎ (𝑖) (𝑡) estimate , called proto-IMF, gets improved according to the following steps:

ℎ (𝑖) 𝑛 (𝑡)
1) Find the local maxima and minima of . ℎ (𝑖) 𝑛 (𝑡)

2) Interpolate, using natural cubic splines, along the maxima and minima of found in the first

ℎ (𝑖) 𝑛 (𝑡)
step in order to form an upper and a lower envelope.

3) Calculate the average of the two envelopes.

4) Obtain a refined estimate of the proto-IMF by subtracting the average found in the ℎ (𝑖) 𝑛 + 1 (𝑡)

previous step.

5) Repeat the process from step 1 unless the stopping criterion has been met.

The sifting process to estimate an IMF of rank i is illustrated in Figure 1. The signal called "arbitrary input" in panel (a) can be the original signal x(t) for the estimation of the first proto-IMF of the first IMF, or the residual signal for the estimation of higher ranking IMFs, or a proto-IMF obtained during the sifting process for any IMF estimation. Panel (a) presents a smooth input signal for clarity. Panels for estimation of the rank i+1 IMF, or for estimation of the residue if the number of extrema is too small. As can be seen by comparing panels (a) and (f), the sifting process has isolated the shortest scales as an AM/FM function. upper and lower envelopes (blue and green lines) and their average (dashed line) of the proto-IMF superimposed on the new input; and (f) the IMF resulting from the sifting process.

To best eliminate the noise in Hs or SSH measurements, the stopping criterion should be met when the number of iterations of the sifting process has reached a predefined number. A value of eight iterations is used to comply with the optimization work performed by [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF] and by Kopsinis and McLaughin (2009). Other stopping criteria may be used. The reason for choosing a given number of iterations is the behavior of the energy distribution of an EMD-processed Gaussian noise according to the IMF ranking, as detailed in section 4. Figure 2 gives an example of an AltiKa Hs record and associated IMFs. The shortest scales are mapped in the first IMF, the second panel from the top in

-Data denoising from Empirical Mode Decomposition

Most denoising methods consist of three steps: projecting the signal on a basis of functions capable of usefully distributing and distinguishing the underlying signal and noise; then making a threshold selection (e.g., hard or soft for wavelet filters, frequency or wavelength cut-off for low-pass filters) in the transformed domain; and finally recovering the signal by back-projecting the modified coefficients into the initial domain. Wavelet and Fourier analyses are often used to effectively decompose signals. With Fourier low-pass filtering, the density energy contained in scales below a given threshold is eliminated. As Gaussian noise contamination can extend over the entire harmonic range, this operation likely eliminates the useful signal at scales below the prescribed cut-off scale and leaves noise above this prescribed scale. This makes it impossible to determine an optimal cut-off wavelength since the actual characteristics of the SNR vary over space and time. Low-pass filtered signals are also subject to spectral ringing (i.e., forced oscillations are generated) and Fourier analysis is not the most effective method for decomposing non-linear and non-stationary signals.

Wavelet decomposition is a more powerful technique, suitable for decomposing non-stationary signals. However, it requires the design of the basis of functions on which the signal is projected, i.e., the choice of the mother wavelet and the number of vanishing moments. This design can be optimized for well-characterized signals but is difficult for Hs or SSH measurements whose geophysical variability and associated noise can have very different signatures. This is the main reason why wavelet shrinkage methods have not yet been successfully used to eliminate noise from altimeter measurements. In this context, the adaptive EMD decomposition provides a simplified framework that overcomes the need to design the basis of functions.

Once the signal is broken down into a set of IMFs, a denoising strategy inspired by those used for wavelet techniques can be applied. Indeed, the analysis to be carried out takes advantage of 1) the well-behaved and predictable distribution of Gaussian noise energy with the IMF basis, 2) the legacy of decades of wavelet-based denoising techniques, and 3) an ensemble average approach to estimate a robust noise-free signal.

Characteristics of the Gaussian noise spectrum with an IMF basis

Reporting on numerical experiments, [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF] showed that in the case of pure fractional Gaussian noise, the first IMF possesses the characteristics of a high-pass filter while the higher order modes behave similarly to a dyadic filter bank, which is reminiscent of what is observed in wavelet decompositions. A bank of dyadic filters is obtained by decomposing the input signal into octave frequency bands in different channels, for which, as descending the frequency scale, the successive frequency bands have half the width of their predecessors. For the EMD decomposition of a white noise, the logarithm (base 2) of the number of zero-crossings (which is an approximate indication of the average frequency of each mode) is shown to decrease linearly with the increasing rank of the IMF, with a slope close to -1. [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF] further found that, for IMFs ranking > 1, the IMF's power spectra are self-similar, i.e., the power spectra of all IMFs could be reduced to a single curve when properly renormalized. The result is a well-defined behavior of the decrease of noise energy with increasing IMF ranking. The Gaussian noise variance projected onto the IMF basis can be modeled as follows:

(3)

𝑣𝑎𝑟(ℎ 𝑛 (𝑡)) ∼ 2 (𝛼 -1)𝑛
where n > 1 is the IMF rank and depends on the Hurst exponent of the fractional Gaussian noise 𝛼 (i.e., = 0.5 for an uncorrelated noise, e.g., white noise; 0.5 for an autocorrelated noise).

𝛼 𝛼

With the EMD basis, noise energy decreases rapidly with the increasing IMF rank: ~ 59, 20.5, 10.3, 5.2, 2.6% of total energy for the first five IMFs, respectively. The first four IMFs represent ~ 95% of the noise energy.

The reasoning behind the EMD denoising approach is derived from this result. Since the IMF energy resulting from the decomposition of a noise signal is predicted from Eq. ( 3), then in real cases of signals containing both information and additive Gaussian noise, a significant difference between the predicted noise energy and the energy of the IMF signal indicates the presence of useful information.

Eq. (3) then gives the expected noise energy in each IMF to determine the different thresholds below which signal fluctuations are associated with noise.

As an illustration, Figure 3 shows how the EMD filter bank distributes a white noise signal and the Jason-2 altimeter Hs signal in the Agulhas Current region. The standard deviation of noise was adjusted to fit the Hs background noise at scales < 20 km. In agreement with [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF], for a white noise, the EMD filter bank is composed of a high-pass filter, IMF1, and a dyadic filter bank for higher ranking IMFs. It is interesting to note that a similar structure is observed when EMD is applied to the Hs along-track signal. Hs IMF1 resembles the noise IMF1, confirming that it contains mainly noise. A higher power spectral density at scales > 20 km is associated with a higher noise modulation, since the noise amplitude depends on Hs and therefore varies along the altimeter track and with very large geophysical Hs gradients, which are common in the Agulhas current region. Pure noise and Hs higher ranking IMFs share the same frequency ranges and, in the logarithm representation used, the SNR increases very rapidly with the increase of IMF rank. This gives the practical rule for denoising since the signal modulation in each IMF can be compared to the noise energy expected in a noise-only IMF of the same rank. Hs IMF2 and IMF3 show comparable energy at the 70-km scale, below which IMF2 energy dominates. On average, Hs IMF2 is noise only at a 30-km scale, and the SNR increases to ~ 0.5 at a 50-km scale which means that some geophysical signals can be retrieved locally under conditions where the SNR is still larger. For the altimetry dataset, the white noise hypothesis (H = 0.5) is made, following studies showing a near-white noise spectrum (Zaron and de Carvalho, 2015;[START_REF] Xu | Global variability of the wavenumber spectrum of oceanic mesoscale turbulence[END_REF][START_REF] Sandwell | Retracking ERS-1 altimeter waveforms for optimal gravity field recovery[END_REF]. The reference noise energy level is then estimated from the first IMF, which contains mainly noise. The theoretical noise energy level of each IMF of rank n > 1 is then given by:

(4)

𝐸 𝑛 = 𝐸 1 0.719 2.01 -𝑛
The two constants may vary slightly depending on the number of sifting iterations used to calculate the IMFs. The values in the above equation are optimized and correspond to eight sifting iterations [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF]. For the estimation of , a robust estimator is the absolute median deviation 𝐸 1

(MAD) from zero, calculated as follows:

(5)

𝐸 1 = (𝑚𝑒𝑑𝑖𝑎𝑛 |𝑛 1 (𝑡)| 0.6745) 2
where n 1 is the noise contained in IMF1, evaluated from a wavelet analysis as described in the following section. The formulation based on a scaled median value is well suited for estimating noise variance, with the factor 0.6745 being a scaling for estimating the noise standard deviation for a Gaussian distribution [START_REF] Johnstone | Wavelet Threshold Estimators for Data with Correlated Noise[END_REF]. It is robust to cases where the analyzed signal, here n 1 , contains residual values associated with a small amount of geophysical information (large Hs gradients) or outliers [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

Signal denoising

The denoised signal y(t) is then reconstructed as :

(6)

𝑦(𝑡) = ∑ 𝑀 2 𝑖 = 𝑀 1 ℎ (𝑖) (𝑡) + ∑ 𝐿 𝑖 = 𝑀 2 + 1 ℎ (𝑖) (𝑡) +𝑑(𝑡)
where h(t) represents the IMFs associated with signal x(t), and the tilde over h(t) indicates the IMFs that are effectively denoised with the threshold technique. M₁ and M₂ are two parameters that give the flexibility to exclude lower-order IMFs than M₁, if they are thought to contain only noise, and to give optional denoising for higher-order IMFs than M 2 , understood as containing very little noise. Indeed, in the case of a white noise signal, the fourth component contains four times less noise energy than the second, and the first four IMFs contain 95% of the total noise energy. This strongly "colored" nature of the noise with the EMD basis has the advantage that few IMFs of a real signal are significantly contaminated by noise, with the SNR increasing rapidly with IMF rank as shown in Figure 3. For M₁, a value of 1 was chosen for the analyzed dataset, because the strong transient or nonlinear features associated with the large Hs or SSH gradients are important geophysical information to be recovered, which are likely to appear in the first IMF. This is illustrated in Figure 4 2) Perform the IMF1 noise filtering to separate the IMF1 stochastic noise, n 1 (t), from the targeted signal and possible outliers. This is done using the wavelet denoising approach described in [START_REF] Huang | Deterministic/stochastic wavelet decomposition for recovery of signal from noisy data[END_REF]. Indeed, IMF1 is a purely oscillating function with a limited bandwidth, dominated by noise, and sharing the same properties and almost the same frequency range everywhere. It makes an adequate wavelet design much easier and robust 3) Perform a reconstruction of signal xs(t) by adding the residual signal, after eliminating noise and outlier in IMF1, with higher order IMFs 4) Randomly modify the positions of the noise n 1 (t) to obtain a new noise realization n k (t) and a new noisy signal x b (t) = xs(t)+n k (t) 5) Perform an EMD expansion of x b (t) 6) Carry out the denoising of IMFs ranking from M₁ to M₂ by hard thresholding, where the threshold is based on E n (Eq. 4 and 5) and reconstruct a denoised signal x k (t) using Eq. ( 6). The advantage of hard thresholding is that it preserves strong gradients better than soft thresholding (Kopsinis and McLaughin, 2009). 7) Iterate k times steps 4 to 6 in order to obtain a set of denoised signals 8) Make an ensemble average of the x k (t) to obtain a robust denoised signal and an uncertainty estimate calculated as the standard deviation of x k (t)

Step 2 of the proposed algorithm is IMF1 denoising with a double objective: first, to evaluate the zero-mean stochastic noise that will be used to generate the ensemble of noisy signals by random permutation of the noise samples; second, to identify the high SNR signal in IMF1 that will be kept in the reconstructed signal. High flexibility is provided by the use of [START_REF] Huang | Deterministic/stochastic wavelet decomposition for recovery of signal from noisy data[END_REF] wavelet shrinkage method. These authors developed a Bayesian approach to denoise various signals, assuming that it can be the sum of the underlying targeted signal composed of a piecewise-smooth deterministic part and a stochastic part, plus a noise as a zero average stochastic part. For the analyzed altimeter dataset, the signal may also contain outliers, which can be part of IMF1. Outliers can be caused in particular by rain and sigma0 blooms, resulting in corrupted altimeter waveforms and abnormal Hs values, which are difficult to edit. The impact of outliers on wavelet analysis can often result in high values of the coefficients at first level (the finest details) that will be targeted as a useful signal by the wavelet shrinkage algorithm. The best way to eliminate these outliers from both the signal and noise estimations, is to analyze the IMF1 with the wavelet bank in two different steps. First, for the estimation of the useful signal in IMF1, the first wavelet scale containing the finest details is systematically rejected in the reconstruction of the signal. It acts as a low-pass filter with a ~ 25 km cut-off wavelength, the wavelength below which the spectral hump dominates [START_REF] Dibarboure | Investigating short-wavelength correlated errors on low-resolution mode altimetry[END_REF]. However, outliers can also corrupt the higher-level wavelet scales, so better editing of the data source would still be necessary. As shown in Figure 4, the strongest geophysical Hs gradients are projected in the largest scale levels of the wavelet decomposition and are therefore not impacted by this low-pass filtering. Second, for the IMF1 noise estimation, all scale levels are used to separate noise, associated with wavelet coefficients below the prescribed threshold, from useful information and outliers, associated with wavelet coefficients above the prescribed threshold. For the wavelet basis design, Symmlet-4 is chosen as the mother wavelet.

Step 2 of the algorithm then provides an estimate of the IMF1 noise, E 1 , in order to calculate the noise energy E n in the other IMFs from Eq. ( 4), and the associated thresholds T n to be used in step 6). T n are derived from E n following the rules that are generally adopted in wavelet denoising techniques. [START_REF] Mallat | A wavelet tour of signal processing[END_REF] showed that the maximum amplitude of Gaussian noise with variance has a very high 𝜎 2

probability of being just below , where N is the number of samples. Indeed, a greater 𝑇 = 𝜎 2𝑙𝑜𝑔𝑁 amplitude of noise is likely to occur with increasing sample size. Following this threshold formulation, longer altimeter records also have a higher probability of including higher Hs values for which the associated noise is higher, which requires a higher threshold value. [START_REF] Mallat | A wavelet tour of signal processing[END_REF] has shown that this formulation is too restrictive and that a lower threshold is a better choice for reducing the risk of filtering the underlying signal. The threshold formulation is then refined by introducing a factor A, giving for an IMF of rank n and sample size N:

(7) 𝑇 𝑛 = 𝐴 𝐸 𝑛 * 2𝑙𝑜𝑔𝑁 with A < 1 being a fine adjustment factor, and the logN factor is constrained to vary within a limited range to account for a possible change in noise amplitude statistics for longer data records, as explained above. A standard value of 0.7 was determined for A in the current Hs denoising algorithm. This parameter can be refined if necessary to account for the differences in the different altimeter missions. Besides, the sensitivity of the results to a prescribed threshold is in practice reduced by the ensemble average performed to estimate a denoised robust signal.

A detailed description of the interval hard thresholding operation can be found in Kopsinis and McLaughin (2009) and subsequent publications dealing with the wavelet thresholding technique. Its adaptation to IMF denoising is illustrated in Figure 5. Black dotted oscillations between the thresholds indicated by the green lines are canceled, and the red dotted line represents the recovered signal.

-Results

Examples

Figures 6 to 7 show the results for the two altimeter records presented in the previous sections.

On the upper panels, raw data (blue solid line), EMD-denoised data (black solid line) and data filtered using the Lanczos filter with a cut-off wavelength of 60 km (red dashed line) are displayed. For the long record in Figure 6, the details do not appear and the two filtered signals seem to be in good agreement. The bottom panel shows the noise that has actually been eliminated. The estimated noise amplitude varies along the track, showing an expected dependence on Hs, with values more than twice as high for Hs > 5 m as for Hs ~ 2 m. The example in Figure 7 is a particularly good illustration of the capability of EMD-based denoising to map very large gradients and avoid smoothing the maximum values. These large Hs values and associated gradients are the result of energetic wave / current interactions [START_REF] Quilfen | Storm waves focusing and steepening in the Agulhas current: satellite observations and modeling[END_REF][START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF]. Since the algorithm is based on measurements extrema detection and removal of the estimated noise, rather than on smoothing, the denoised maximum Hs reaches a peak close to 9 m, about 20 cm less than the raw data. This contrasts with Lanczos filtering at a 60 km cutoff (value used in CMEMS operational products), which significantly modifies gradients and peak values, regardless of how smooth the roll-off is. As tested, a 30 km cut-off for Lanczos filter will effectively restore the observed Hs gradient and maximum, but at a cost of a large number of artifact oscillations produced elsewhere. This expected behavior of low-pass filtering [START_REF] Emery | Data analysis methods in physical oceanography[END_REF] is illustrated further in Figures 8 and9.

Figure 8 shows a Cryosat-2 pass across the island of Fiji, with an abrupt change in sea state on the northern side of the island, between the lagoon area and the high seas. The reason for such a marked Hs gradient is that the lagoon is protected by the Yasawa coral reefs. The data point in this area is missing from the original data as it was corrupted by the ground, but the altimeter track segments on both sides were processed as a continuous segment for illustration purposes. Indeed, high resolution altimetry can provide data very close to the coast and situations involving very rapid changes in fetch conditions will often be encountered. The lower right panel of Figure 8 shows the raw and denoised Hs when using the EMD and Lanczos processing approaches. Although the Lanczos filter is specifically designed to reduce the Gibbs effect, at the cost of gradient smoothing, residual artifact oscillations are still produced on both sides of the gradient and their amplitude depends on the number of coefficients (size of the convolution kernel) used to control the desired roll-off near the cut-off wavelength. The solid green line shows the results for a number of coefficients locally adjusted to best reduce the Gibbs oscillations, and the solid red line shows the results for a number of coefficients ensuring a sharp cut-off of the filter. This choice has a significant impact, both locally and globally, as discussed in the following section. EMD denoising provides a more realistic estimate of the Hs profile. Indeed, as shown in the upper right panel of Figure 8, the IMF1 wavelet analysis effectively separates noise from the frontal signature that occurs near 16.7° of latitude. To avoid artifacts and data loss near the edges of the data record for the Lanczos processing, the same mirror padding was applied for both EMD and Lanczos processing. In addition, the data filtered by the Lanczos filter contain a significant amount of statistically insignificant modulations, constrained by the Fourier analysis; this is illustrated in Figure 9, which shows a Cryosat-2 pass across the Galapagos islands. This is an area where the SAR mode is systematically put into operation, resulting in a 1-Hz Pseudo Low Resolution Mode (P-LRM) noise increase as shown in the right panel of the figure. In this low SNR context, small mesoscale variability is embedded in noise with similar or greater amplitude. In such conditions, Hs derived from EMD varies more smoothly as a result of the analyzed low SNR, while Hs derived from the Lanczos filter shows forced oscillations, for example near 5°S latitude. The EMD result is constrained by the local SNR analyzed while the Lanczos result is strongly constrained by the Fourier analysis framework.

-Analysis of the denoised signal and noise signatures

Figure 10 shows the Hs wavenumber spectra in the Agulhas current region (10 -35 E ; 45 -33 S) for the period 2014-2016. The raw measurements spectra flatten at scales < 100 km where the measurements are more significantly contaminated by noise. The differences in the average noise levels of the different altimeters are as expected, with Saral data being less noisy than Cryosat-2 data.

Agulhas is one of the regions where the Cryosat-2 Synthetic Aperture Radar mode is systematically used, resulting in noisier data for the 1-Hz low rate mode. In other regions where SAR mode is not used, Cryosat-2 noise should be lower than Jason-2 noise. The spectra of denoised measurements are almost identical for all three altimeters, which gives a good indication that the EMD processing worked satisfactorily. The PSD obtained from the Lanczos filtering is shown for a sharp (solid green line) and smooth (dashed green line) roll-off at the cut-off wavelength. As expected, the results are very different, which raises the question of the choice of this control parameter. Surprisingly, this has not been addressed in studies using the Lanczos filter for altimeter data filtering. For scales between 300 and 30 km, the PSD slope calculated for data filtered by EMD is not linear. However, it compares favorably with a k -3 law, result in agreement with the study by [START_REF] Ardhuin | Small-scale open-ocean currents have large effects on wind-wave heights[END_REF] The total noise spectrum exhibits the same behavior for all three altimeters, with a "plateau" for scales below ~ 25 km, at the same PSD level as the raw Hs signals. This part of the signal is completely suppressed as explained above, as it essentially corresponds to the sum of the instrumental Gaussian white noise and the hump artifact [START_REF] Dibarboure | Investigating short-wavelength correlated errors on low-resolution mode altimetry[END_REF]. The PSD of total noise gradually diminishes with the increase in scale, indicating that significant low-frequency noise energy is removed from IMFs of order > 1, which can be evaluated from the difference between the total noise and the IMF1 noise PSDs (Figure 10, left and right panels).

The IMF1 noise PSD, which takes into account high frequency Gaussian noise at scales below about 25 km, is slightly lower than the raw Hs PSD. This difference should correspond to outlier signatures and is small because the sigma0 bloom occurrence is weak in the Agulhas region [START_REF] Tournadre | A satellite altimeter model for ocean slick detection[END_REF]. A different characterization of the measurement noise, essentially free of outliers, can then be made using the high frequency noise estimated from the processing of IMF1. Indeed, EMD analysis is an efficient way to isolate the high frequency part of the total noise variance in IMF1, which theoretically represents about 59% of the total noise for a white noise signal. Since the noise variance left in other IMFs is simply related to IMF1 noise variance through Eq. ( 3 As mentioned earlier, outliers affect the EMD analysis. A more thorough verification of the data is therefore carried out before the noise analysis. The data are rejected when the sigma0 is greater than the 99.5 th percentile, when the IMF1 absolute value is greater than the 99.9 th percentile, and when the distance to the coast is less than 50 km. The sigma0 threshold allows the editing for sigma0s, and corresponds to 18.9, 18.7, 19.4 dB for Cryosat-2, Jason-2, and Saral/AltiKa, respectively. The standard deviation of IMF1 noise in one-meter Hs bins is shown in Figure 11. The results are consistent with simulations performed with recent retracking algorithms [START_REF] Thibault | Convergent solutions for retracking conventional and Delay Doppler altimeter echoes[END_REF], showing an increasing standard deviation of noise with the increase in Hs and an overall lower standard deviation of noise for Saral/AltiKa. The noise levels of Cryosat-2 and Jason-2 are of similar amplitude although this global analysis includes areas where the Cryosat-2 SAR mode is used, for which noise in 1-Hz low resolution mode is higher. Figure 11 (lower panel) shows the probability density function (PDF) of the IMF1 noise amplitude for the three altimeters, and that obtained from Monte-Carlo simulations of the IMF1 amplitude for Gaussian noise with zero-mean and standard deviation of 12 cm (green dots) and 7 cm (magenta dots). Two hundred realizations of a Gaussian noise with 512 samples are performed for the Monte Carlo simulation. The good agreement between the measured and simulated PDFs is a good indication that the chosen denoising approach is consistent. A best fit between the altimeters and the simulated PDFs gives an average estimate of the Hs measurement noise standard deviation at 12., 12.75, and 7. cm for Cryosat-2, Jason-2, and Saral/AltiKa, respectively. Although not strictly comparable, these values are consistent with the 9, 12, and 6 cm values obtained from the standard deviation of the high-resolution measurements used to calculate 1-Hz Hs [START_REF] Taburet | WAVE Thematic Assembly Centre Quality Information Document[END_REF], since our estimate for Cryosat-2 includes the areas with SAR operations.

The EMD method therefore also allows noise analysis and monitoring over time and under local conditions. The noise characteristics presented should be significantly improved by the application of improved re-tracking algorithms (e.g., [START_REF] Passaro | ALES: A multi-mission subwaveform retracker for coastal and open ocean altimetry[END_REF][START_REF] Thibault | Convergent solutions for retracking conventional and Delay Doppler altimeter echoes[END_REF], as planned in the ESA sea state CCI project.

-Summary

An adaptive method to efficiently eliminate and characterize altimeter data noise was defined and applied. Our method is essentially built on Empirical Mode Decomposition (EMD), which projects a signal on a basis of empirical AM/FM functions, called Intrinsic Modulation Functions (IMFs), whose properties offer a new way of analyzing and denoising along-track altimeter measurements. When applied to a Gaussian noise signal, EMD provides a set of IMFs with a predictable energy distribution. This property was therefore used to define rules for altimeter data denoising, which are based on existing wavelet filtering techniques. The chosen approach simply assumes that the Hs along-track variability is the sum of an underlying signal composed of a piecewise-smooth deterministic part and a stochastic part, allowing the use of a Bayesian analysis, and a zero-mean Gaussian noise. The stochastic hypothesis for the wave field is a framework common to many studies, as well as in traditional wave models such as WAVEWATCH-III, to take into account the random nature of the wave field.

The method performs a local SNR analysis for each processed data record. For each record, whose length varies from about 50 km to 1500 km depending on the continuity of the data, the final noisefree altimeter signal is estimated as an ensemble average of several signal realizations. A fine adjustment of the result is possible thanks to a global scaling parameter, making it possible to adapt the thresholds used to separate the noise from the signal. It allows, if necessary, to optimize the results according to the different sensors on-board the different altimetry missions. Complementary to methods based on crossovers analysis, EMD provides also a mean to effectively estimate the noise characteristics for the various altimeter sensors. The noise variance parameter derived from the EMD should therefore be useful for noise monitoring over time and under local conditions and can provide practical information in various applications such as altimeter data assimilation to numerical models.

The proposed method is free of systematic artifacts, preserves the amplitude of spatial gradients and extreme values, and eliminate the noise over the whole frequency range, which is not the case for the boxcar or low-pass filters currently applied to altimetry products distributed to users. Indeed, such filters only remove some of the high-frequency noise, which is also the result of approaches that infer a correction to eliminate correlated errors from other aspects of the waveform data (Zaron and de Carvalho, 2015;[START_REF] Quartly | Removal of covariant errors from altimetric wave height data[END_REF][START_REF] Tran | Assessing the effects of sea-state related errors on the precision of high-rate Jason-3 altimeter sea level data[END_REF], leaving a substantial amount of low-frequency noise in the data. Low-frequency noise remaining in the signals after smoothing will still interfere with the precise measurement of geophysical peaks height, width, and position (O'Haver, 2019).

The EMD-based denoising method robustness stems from the self-consistent, data-driven approach, and from the ensemble average computed as the resulting denoised signal. The method developed for processing altimeter data benefited from previous optimization work and sensitivity studies carried out by [START_REF] Flandrin | Empirical mode decomposition as a filter bank[END_REF] Kopsinis and McLaughin (2009). Since the EMD-based denoising method does not contain an implicit smoothing operator, sparse data may still appear noisy in particular when the data editing performed before denoising did not identify large errors. It can occur mainly when altimeter waveforms are corrupted by rain, surface slicks or very light winds generating sigma0 blooms, or during transitions between land, sea ice and sea. However, the use of more efficient data editing and retracking algorithms should mitigate these errors (e.g., [START_REF] Passaro | ALES: A multi-mission subwaveform retracker for coastal and open ocean altimetry[END_REF][START_REF] Thibault | Convergent solutions for retracking conventional and Delay Doppler altimeter echoes[END_REF].

As reported, signals down to scales of nearly 30 km can be recovered, provided that the signal-to-noise ratio is sufficient. This is particularly the case in the presence of very strong gradients and in main dynamic and/or eddy-rich regions. Although recent studies (Ardhuin et al., 2017, Quilfen and[START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF] reveal that wave-current interactions can trigger a significant increase in energy below the 100-km scale, very little is known and very few measurements are available for this scale of sea state variability. Until now, the sea state variability at scales below 100 km was often considered to be noise only [START_REF] Zaron | Identification and reduction of retracker-related noise in altimeter-derived sea surface height measurements[END_REF].

Analysis of Jason-2, Cryosat-2, and Saral/AltiKa significant wave height power spectral density and of the spatial distribution of variance at short scales shows that the proposed method gives consistent results. As a newly developed tool for altimetry data, the EMD-based system could certainly require further testing and optimization, in particular for better screening of outliers that can shape the energy spectrum and bias noise statistics. The proposed method is a step forward to better exploit this unique set of sea state observations, which now covers more than 25 years, and is easily applicable to other altimetric observations, such as sea level and radar cross-section.

  (b) to (e) show the steps of the sifting process for two iterations, i.e., searching for extrema in panel (b), determining the upper and lower envelopes in panel (c), and estimating a refined proto-IMF by removing the average of the upper and lower envelopes in panel (d). Panel (e) is like panel (c) but an iteration later, showing that the average of the two envelopes is approaching 0. Panel (f) shows the IMF resulting from the sifting process. This IMF signal is removed from the input signal in panel (a)

Figure 1

 1 Figure 1: Simplified diagram of the sifting process for two successive iterations, panels a) to e), and the IMF obtained at the last iteration, panel f). (a) an arbitrary input; (b) identified local maxima (blue diamonds) and minima (green circles) superimposed on the input; (c) upper and lower envelopes (blue and green lines) and their average (dashed line) superimposed on the input; (d) proto-IMF (new input calculated as the difference between the solid line and the dashed line in Figure 2c), which must be refined; (e)

Figure 2 ,

 2 Figure 2, to likely contain mostly noise information. The bottom panel shows the last IMF and the

Figure 2 :

 2 Figure 2: EMD expansion in L IMFs and residue for a segment of AltiKa Hs measurements (October 7, 2016). (upper panel) raw Hs ; (intermediate panels) L-1 IMFs ; (bottom panel) last IMF and residue. As obtained, six IMFs plus the residue describe the total signal.

Figure 3 :

 3 Figure 3: Mean Power Spectral Density (PSD) of the first four IMFs for white noise (red curves) and Jason-2 Hs along-track measurements (black curves), and mean PSD of the corresponding noisy (solid blue line) and denoised (dashed blue line) Jason-2 Hs measurements. The PSD is the average of PSDs computed over all data segments covering the years 2014 to 2016, and the Agulhas region (10ºE -35ºE; 45ºS -33ºS).

  , which shows the raw and filtered Jason-2 Hs measurements (upper panel) with very large gradients (~ 1.5 m / 10 km) triggered by wave-current interactions. Locally, close to the maximum of Hs, the IMF1 shows a large amplitude modulation well above the average noise level.

Figure 4 :

 4 Figure 4: (top) Jason-2 Hs along-track raw (blue) and filtered (black) measurements as a function of latitude; (bottom) associated IMF1 obtained from EMD expansion.

Figure 5 :

 5 Figure 5: IMF amplitude before (solid blue line) and after (red dots) thresholding. Black dots represent the canceled signal and solid green lines indicate the threshold values.

  Figure 6: AltiKa 1-Hz along track Hs measure ments (Octobe r 7, 2016). (upper panel) raw Hs (solid blue line), EMD-denoised Hs (solid black line), Lanczos-filtered Hs (dashed red line). (bottom panel) noise estimated as the difference between raw Hs and EMD-denoised Hs. The abscissa spatial sampling is about 7 km.

Figure 7 :

 7 Figure 7: Jason-2 1-Hz Hs measurements (February 29, 2016). (upper panel) raw Hs (solid blue line), EMDdenoised Hs (solid black line), Lanczos-filtered Hs with smooth (solid green line) and sharp (solid red line) roll-

Figure 8 :

 8 Figure 8: Left panel: Cryosat-2 1-Hz Hs measurements over the Fiji island (January 1, 2016). Right upper panel: IMF1 (blue line) and noise estimated from wavelet analysis of IMF1 (red line). Right bottom panel: raw Hs (solid blue line), EMD-denoised Hs (solid black line), Lanczos-filtered Hs with smooth (solid green line) and sharp (solid red line) roll-off as a function of latitude.

Figure 9 :

 9 Figure 9: Left panel: Cryosat-2 1-Hz Hs measurements over the Galapagos Islands (October 7, 2016). Right panel: raw Hs (blue solid line), EMD-denoised Hs (solid black line), Lanczos-filtered Hs with smooth (solid green line) and sharp (solid red line) roll-off as a function of latitude.

  based on high resolution numerical simulations of the WaveWatch-III model. Below a scale of about 25 km, the PSD vanishes for the three altimeters as a consequence of the cancellation of IMF1's finest wavelet scales (as discussed at the end of section 4). The results can be analyzed in more detail from the interpretation of the total noise PSD (see Figure10left panel), and of the IMF1 noise PSD (see Figure10right panel). The total noise is calculated as raw Hs minus filtered, which therefore include outliers.

  ), IMF1 noise gives a good representation of the characteristics of the total noise affecting the Hs noisy measurements. It is also the noise figure used to perform the ensemble averaging operation that ultimately gives the denoised signal and associated uncertainty. The results are presented in Figure 11.

Figure 10 :

 10 Figure 10: Power spectral density (PSD) as a function of wavenumber (x-axis, labels in km) for altimeter Hs measurements: SARAL/Altika (black), Cryosat-2 (red) and Jason-2 (blue), period 2014-2016, Agulhas region. Dashed lines show raw 1 Hz data, thick solid lines show filtered data, and thin solid lines show the signal removed by the whole process (left panel) and for the IMF1 noise (right panel). The solid (dashed) green line gives the PSD for Jason-2 data processed with a Lanczos filter of 60 km cutoff wavelength with sharp (smooth) roll-off. The solid black line gives the k -3 law between 300 km and 30 km wavelength (right).

Figure 11 :

 11 Figure 11: (top) Standard deviation (m) of IMF1 noise amplitude as a function of Hs (m) in one-meter bins; (bottom) probability density function (%) of IMF1 noise amplitude (m). SARAL/Altika (black), Cryosat-2 (red) and Jason-2 (blue).The dots map the PDF, obtained from Monte-Carlo simulations, of the IMF1 amplitude for gaussian noise with zero-mean and standard deviation of 12 cm (green dots) and 7 cm (magenta dots).

Figure 12

 12 Figure 12 shows the spatial distribution of the IMF1 variance, calculated as the sum of the square amplitudes in 1 o x1 o longitude/latitude boxes. Patterns are clearly apparent in the IMF1 variance

Figure 13 :Figure 14 :

 1314 Figure 13: Statistical behavior of the  uncertainty parameter (meters) for SARAL/Altika (black), Cryosat-2 (red) and Jason-2 (blue): (top panel) probability density function; (middle panel) Hs dependency; (bottom panel) radar cross-section, ( ) dependency.𝜎 0

-IntroductionFor more than

years, satellite altimeter measurements have supplied sea surface height (SSH), significant wave height (Hs), and radar cross-section (sigma0) measurements, building a unique dataset on these aspects of the ocean. Although this data offers potential for significant progress in various research-oriented and operational applications, its use can be hampered by technological differences between the sensors used and the different processing and re-tracking

Acknowledgments

This study was conducted within the Ocean Surface Topography Science Team (OSTST) activities. A grant was awarded to the SILLAGE project by the TOSCA board in the framework of the CNES/EUMETSAT call CNES-DSP/OT 12-2118. The authors thanks Y. Kopsinis for having made the EMD-denoising software available (visit http://cgi.di.uoa.gr/~kopsinis/software.html); the version modified to process altimeter data can also be obtained upon request. The authors are grateful to the two anonymous reviewers who helped to improve the manuscript. Altimeter data were provided by Institut Français de Recherche pour l'Exploitation de la MER (IFREMER) ftp://ftp.ifremer.fr/ifremer/cersat/products/swath/altimeters/waves/).

increased sea state variability resulting from wave-current interactions that dominate Hs variance at scales below 100 km [START_REF] Ardhuin | Small-scale open-ocean currents have large effects on wind-wave heights[END_REF][START_REF] Quilfen | Storm waves focusing and steepening in the Agulhas current: satellite observations and modeling[END_REF][START_REF] Quilfen | Ocean surface wave-current signatures from satellite altimeter measurements[END_REF] .

-Uncertainty associated with the denoised significant wave height estimate

In the previous section, the noise derived from EMD processing and associated with standard 1Hz measurements is analyzed and characterized. However, this does not characterize the uncertainties associated with the "denoised significant wave height" parameter. Such uncertainties could be assessed using conventional approaches, such as comparisons with a reference dataset or a crossover analysis, but this will give an estimate of the average uncertainty rather than an estimate associated with each measurement. The EMD ensemble average approach allows this to be done in line with the estimate of denoised Hs. Indeed, the uncertainty  can be calculated for each measurement along the altimeter track as the standard deviation of the ensemble of denoised signals used to calculate the average denoised Hs. The probability density function of  is shown in Figure 13, top panel. It reflects the results presented in the previous section: the noisier the standard 1 Hz measurements, the greater the uncertainties associated with the denoised measurements. The central panel also shows that  decreases with decreasing significant wave height, except for the lower values for which Jason-2 and Cryosat-2 values are slightly higher. This can be analyzed in more detail by computing  as a function of the radar cross-section, , which is mainly related to wind speed. At low values, corresponding 𝜎 0 𝜎 0 to strong winds and large associated Hs,  values are larger as expected.  values decrease with decreasing , reaching a minimum close to 12dB, then increase with decreasing wind speed. 𝜎 0 Uncertainties associated with sea state measurements are also expected to be greater under light wind conditions, but Jason-2 and Cryosat-2  values increase more rapidly. This is probably related to the fact that the light wind zones are climatologically associated with precipitation and surface slicks.

Under such conditions, the better Saral/AltiKa SNR combined with the larger rain contamination in the Ka band makes data editing and noise filtering more efficient. Indeed, the number of edited data, before data denoising, is significantly higher for Saral/AltiKa.  is then a practical uncertainty parameter, taking into account the different errors likely to affect the denoised Hs data and their dependence on local conditions, which should be useful in particular for assimilation in numerical wave models. To further highlight its behavior, Figure 14 shows the mean geographical distribution of  computed over a three-year time period (2014)(2015)(2016). As expected, the uncertainty is related to the altimeter measurements signal to noise ratio: it is higher for Jason-2 than for Saral/AltiKa, it is higher under high sea state conditions that predominate in southern oceans and in the North Atlantic Ocean, and is higher in the regions of increased sea state mesoscale variability as depicted in Figure 12 (right panel) for which the SNR is effectively lower.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: