Laura Gómez-Navarro

Emmanuel Cosme

Julien Le Sommer

Nicolas Papadakis

Ananda Pascual

Development of an image de-noising method in preparation for the Surface Water and Ocean Topography satellite mission

Keywords: SWOT, De-noising, Variational regularization, western Mediterranean

In a near future, the Surface Water Ocean Topography (SWOT) mission will provide images of altimetric data at kilometric resolution. This unprecedented 2-dimensional data structure will allow the estimation of geostrophy-related quantities that are essential for studying the ocean surface dynamics and for data assimilation uses. To estimate these quantities, i.e. compute spatial derivatives of the Sea Surface Height (SSH) measurements, the small-scale noise expected to affect the SWOT data must be smoothed out while minimizing the loss of relevant, physical SSH information. This paper introduces a new technique for de-noising the future SWOT SSH images. The de-noising model is formulated as a regularized least-square problem with a Tikhonov regularization based on the first, second, and third-order derivatives of SSH. The method is implemented and compared to other, convolution-based filtering methods with boxcar and Gaussian kernels. This is performed using a large set of pseudo-SWOT data generated in the Western Mediterranean Sea, from a 1/60 • simulation and the SWOT simulator. Based on Root Mean Square Error and spectral diagnostics, our de-noising method shows a better performance than the convolution-based methods. We find the optimal parametrization to be when only the second-order SSH derivative is penalized. This de-noising reduces the spatial scale resolved by SWOT by a factor of 2, and at 10 km wavelengths the noise level is reduced by 10 4 and 10 3 for Summer and Winter respectively. This is encouraging for the processing of the future SWOT data.

Introduction

The Surface Water Ocean Topography (SWOT) [START_REF] Fu | The SWOT Mission Science document[END_REF] mission will provide an unprecedented two-dimensional view of ocean surface topography at a pixel resolution of 2 km. The launch is scheduled for 2021. SWOT's wide-swath altimeter, based upon SAR interferometry technology, will measure Sea Surface Height (SSH) over a 120-km wide swath with a 20-km gap at the nadir. The satellite will also carry a conventional nadir altimeter. SWOT will evolve on two different orbits: the first 3 months of scientific data production will be dedicated to a fast-sampling phase, where the repeat cycle will be of 1 day. Then, the satellite will be moved to its nominal orbit with a 20.86-day repeat cycle. SWOT is a multi-disciplinary hydrology and oceanography mission, and here we focus on the latter.

The main oceanographic objective of SWOT is to observe the geostrophic fine-scale circulation at the global scale [START_REF] Fu | Observing Oceanic Submesoscale Processes From Space[END_REF][START_REF] Morrow | Global observations of fine-scale ocean surface topography with the Surface Water and Ocean Topography (SWOT) Mission[END_REF]. The measurement system is designed to resolve ocean circulation patterns at scales down to 15 km, whereas the present-day constellation of conventional altimeters only resolves scales of 150-200 km and above [START_REF] Morrow | Global observations of fine-scale ocean surface topography with the Surface Water and Ocean Topography (SWOT) Mission[END_REF]. In addition to potentially unexpected Preprint submitted to Elsevier discoveries, this order-of-magnitude gain in resolution will help quantifying several oceanic processes much more accurately than today. Among those processes are vertical motions, which are key to the vertical exchanges between the ocean surface and the atmosphere, and between the ocean surface and the deep ocean [START_REF] Thomas | Ocean Modeling in an Eddying Regime[END_REF][START_REF] Taylor | Ocean fronts trigger high latitude phytoplankton blooms[END_REF][START_REF] Lévy | Large-scale impacts of submesoscale dynamics on phytoplankton: Local and remote effects[END_REF][START_REF] Mahadevan | The impact of submesoscale physics on primary productivity of plankton[END_REF][START_REF] Mcgillicuddy | Mechanisms of physical-biological-biogeochemical interaction at the oceanic mesoscale[END_REF][START_REF] Pascual | A multiplatform experiment to unravel meso-and submesoscale processes in an intense front (alborex)[END_REF]; and the dissipation of kinetic energy, which partly determines the climatic role of the global ocean [START_REF] Sasaki | Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere[END_REF][START_REF] Uchida | Seasonality of eddy kinetic energy in an eddy permitting global climate model[END_REF].

The SWOT mission objectives will be reached if we can accurately estimate gridded maps of at least the first and second-order horizontal derivatives of SSH. Altimetry describes the upper ocean dynamics through geostrophy, which involves the horizontal SSH gradients. Geostrophy is a fairly good approximation of mesoscale dynamics, i.e. at scales larger than the for first Rossby deformation radius (about 10-15 km in our region [START_REF] Escudier | Eddy properties in the Western Mediterranean Sea from satellite altimetry and a numerical simulation[END_REF]), for which Rossby numbers are typically smaller than 1. Kinetic energy dissipation is driven by the horizontal strain rates of the ocean surface flow [e.g. 13]. Complete, gridded maps of SSH derivatives are required for climate studies and short-term operational applications. One way to make gridded maps from incomplete SSH observations (including SWOT, but not only) is to assimilate those data into dynamical models. The assimilation of SWOT is expected to be challenging because of the spatially correlated noise, and promising solutions to this rely upon the joint assimilation of SSH and its derivatives [START_REF] Ruggiero | An efficient way to account for observation error correlations in the assimilation of data from the future SWOT High-Resolution altimeter mission[END_REF][START_REF] Yaremchuk | On the approximation of the inverse error covariances of high resolution satellite altimetry data[END_REF]. All these considerations compel the scientific community to strive for getting accurate estimates of SSH derivatives.

Unfortunately, SWOT data will very likely be contaminated by small-scale noise that prevents the direct computation of SSH derivatives. The noise expected to contaminate SWOT measurements gathers several components with different spatial coherences and different amplitudes. Details are provided in the SWOT mission performance and error budget document [START_REF] Esteban-Fernandez | Swot project: mission performance and error budget document[END_REF]. To be prepared to exploit the future SWOT data, the SWOT simulator for ocean science has been developed to simulate realistic realizations of SWOT uncertainties [START_REF] Gaultier | The challenge of using future SWOT data for oceanic field reconstruction[END_REF]. Some are illustrated on Figure 2. Errors due to the satellite roll, the baseline dilation, and the path delay induced by atmospheric humidity, exhibit significant spatial correlations with different characteristic patterns. The system timing error presents errors invariant across-track, but with possible small-scale variations along-track. The KaRIn (Ka-band Radar Interferometer) noise is spatially uncorrelated, with higher amplitudes at nadir and near the edges of the swath. The path-delay component also exhibits small-scale variations due to sharp changes in air humidity.

Efforts have already been undertaken to filter out SWOT small-scale noise by Gómez-Navarro et al. [START_REF] Gómez-Navarro | Swot spatial scales in the western mediterranean sea derived from pseudo-observations and an ad hoc filtering[END_REF]. The authors show that the implementation of a diffusion-based filter allows to retrieve the dynamical spectral signature down to 40-60 km scales (20-30 km in terms of dynamical pattern scales). However, the de-noising approach here is not specifically designed to retrieve SSH derivatives, and we believe there is room for improvement in the scales to be retrieved. This paper presents a method designed to remove the small-scale noise of the future SWOT data, which explicitly relies upon the regularity (bounded variations) of the first three orders of SSH derivatives. This de-noising method is rooted in image restoration techniques of the variational type [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Chang | Spatially adaptive wavelet thresholding with context modeling for image denoising[END_REF][START_REF] Desbrun | Anisotropic feature-preserving denoising of height fields and bivariate data[END_REF][START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. The range of image restoration techniques is extremely wide and diversified. Testing all existing methods is out of reach and irrelevant here. Our approach is then to acknowledge that our image is a smooth physical field with relatively smooth derivatives, and that the estimation of derivatives is an important issue. This consideration guides the design of the de-noising method presented in Section 2. The method involves a set of parameters that must be adjusted. An essential task is to identify optimal sets of parameters.

This study suggests a methodology to identify them. The experimental setup is described in Section 3. Sections 4 and 5 present the results and Section 6 summarizes the study, draws the most relevant conclusions, discusses them, and suggests possible future research paths.

Variational de-noising of SWOT images with penalization of derivatives

Formulation of the image de-noising problem

The primary purpose of image de-noising here is to allow the computation of first and secondorder SSH spatial derivatives of SWOT data as accurately as possible. The two reasons, already mentioned in the introduction, are: (i) these quantities represent geostrophic velocities and relative vorticity, respectively, whose estimation is central to the success of SWOT mission; and

(ii) these quantities can be needed to draw maximum benefits from the assimilation of SWOT data into ocean circulation models [START_REF] Ruggiero | An efficient way to account for observation error correlations in the assimilation of data from the future SWOT High-Resolution altimeter mission[END_REF][START_REF] Yaremchuk | On the approximation of the inverse error covariances of high resolution satellite altimetry data[END_REF]. We therefore propose a method that explicitly constrains these derivatives.

The proposed de-noising model is formulated as a regularized least-square problem with a Tikhonov regularization. The de-noised SWOT image h is searched for by minimizing the following cost function:

J(h) = 1 2 m • (h -h obs ) 2 + λ 1 2 ∇h 2 + λ 2 2 ∆h 2 + λ 3 2 ∇∆h 2 (1)
where represents the L 2 -norm, h obs is the original noisy image (i.e., our observation, the pseudo-SWOT data) , ∇ = (∂/∂x, ∂/∂y) is the gradient operator, and

∆ = ∂ 2 /∂x 2 + ∂ 2 /∂y 2
is the Laplacian operator. Letter m and sign • represent a mask and the entrywise matrix product, respectively. They can be ignored for the present and the next sub-section: their role is discussed in Section 2.3 below. The regularization terms impose regularity constraints on geostrophic velocity, vorticity, and vorticity gradient, respectively. Parameters λ 1 λ 2 and λ 3 must be prescribed. The search for their optimal values is reported in Section 3.3.

Resolution of the variational problem

The variational problem displayed in eq. 1 is solved using a gradient descent method [START_REF] Biemond | Iterative methods for image deblurring[END_REF]. The gradient of J is written:

∇J(h) = m • (h -h obs ) -λ 1 ∆h + λ 2 ∆∆h -λ 3 ∆∆∆h (2)
so that the solution can be reached after convergence of the following iterations:

h k+1 = h k + τ m • (h obs -h k ) + λ 1 ∆h k -λ 2 ∆∆h k + λ 3 ∆∆∆h k (3) Stability of iterations is guaranteed if τ ≤ (1 + 8λ 1 + 64λ 2 + 512λ 3 ) -1 .
In practice, it is taken equal to this value. Two improvements on the method's implementation accelerate the gradient descent: Firstly, iterations are started with a preconditioned image obtained by applying a Gaussian filter onto the original image, including inpainting as discussed in Sections 2.3 and 2.4 (note that h obs remains the original, unfiltered image). Preconditioning considerably speeds up the algorithm convergence, in particular for the inpainted regions. Secondly, iterations are actually implemented with an acceleration of the scheme 3, based on the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], detailed in Appendix B. Iterations are stopped when h k+1 -h k < 10 -9 or if k = 10 4 . Those values have been fixed after a careful search of a trade-off between accuracy and numerical efficiency.

The Laplacian operator is discretized with finite differences using the five-point stencils of the image pixels. As commonly done in image processing, the division by pixel size is ignored;

this also reduces the probability of truncation errors due to operations with terms different by too many orders of magnitude. Pixels located at the boundaries, where the stencil is incomplete, must have a Laplacian value attributed yet; otherwise, the image would become smaller at each iteration of the gradient descent. The implementation of the Laplacian operator follows [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] and is detailed in Appendix A.

Dealing with gaps in the image

An inpainting method is implemented to deal with islands, continents, and the 20-km wide gap at the SWOT nadir, which all represent obstacles to the calculation of the second deriva-

tives

Comparison with convolution-based filters

In Section 4, the image de-noising technique described above will be compared with standardtype filters, namely convolution-based filters. In our experiments we test the two commonly used boxcar and Gaussian convolution kernels, with a large range of parameters, and we shortly refer to the boxcar filter and the Gaussian filter. Their parameters are the box size (or footprint) and the standard deviation for the Gaussian kernel (hereinafter referred to as σ). Gaps in the SWOT swath (lands, islands, and nadir gap) are inpainted to facilitate filtering and to ensure the smoothness of SSH fields. Then, SSH values created in gaps are removed for the evaluation of the methods using the mask m. Inpainting is implemented as follows: (i) Image gaps are filled with zeros; (ii) both the filled image and the mask m are filtered with the same kernel; and (iv) the filtered filled image is divided entrywise by the filtered mask. Note that in an earlier study [START_REF] Gómez-Navarro | Swot spatial scales in the western mediterranean sea derived from pseudo-observations and an ad hoc filtering[END_REF], a Laplacian diffusion filter was experimented. It is not reproduced here, since it is equivalent to the Gaussian filter implemented in this study.

Experimental setup

Simulated SWOT dataset

The input of our database is a 15-month North Atlantic simulation at a resolution of 1/60 The horizontal resolution is between 0.8 and 1.6 km (depending on latitude), and the grid has 300 vertical levels. This NEMO model configuration is referred to as NAtl60 and the source files and codes are available in [START_REF] Molines | meom-configurations/NATL60-CJM165[END_REF]. The particular simulation used herein has been described

in [START_REF] Amores | Up to what extent can we characterize ocean eddies using present-day gridded altimetric products?[END_REF][START_REF] Fresnay | Reconstruction of the 3-d dynamics from surface variables in a high-resolution simulation of north atlantic[END_REF][START_REF] Ajayi | Diagnosing cross-scale kinetic energy exchanges from two submesoscale permitting ocean models[END_REF]. Lastly, the simulation time span is from mid-June 2012 to October 2013. [START_REF] Sommer | Natl60: A north atlantic ocean circulation model dataset based on nemo for preparing swot altimeter mission[END_REF] The SWOT simulator for Ocean Science (version 2.21) [START_REF] Gaultier | The challenge of using future SWOT data for oceanic field reconstruction[END_REF] is run to generate pseudo-SWOT scenes from the NAtl60 simulation. The SWOT simulator first builds the SWOT observation grid, based on the provided satellite orbit. In this study, SWOT grid resolution is fixed at 1 km. After this work started, the resolution of the basic SWOT level 2 SSH data products has been fixed to 2 km, but this small mismatch does not modify the general approach. After building the grid, the simulator reads SSH data from NAtl60 and linearly interpolates them from model to SWOT grid (rendering the variable SSH model). In a last step, it computes random realizations of observation errors and adds them to the interpolated SSH data (rendering SSH obs). Observation errors considered at the moment are KaRIn errors, roll errors, phase errors, baseline dilation errors, timing errors, and errors due to signal alternation by atmospheric humidity. Among these errors, only the KaRIn noise is expected to be spatially uncorrelated. Technically, the SWOT simulator provides simulations of the noise-free SSH observed by SWOT, and of the noisy data that SWOT will actually yield (sum of the former and the noise: SSH obs = SSH model + noises/errors). For the evaluation of image de-noising methods, it thus provides "true", noise-free images (h true ) along with the realistic SWOT data 88) passes. The selected region belongs to the fast-sampling phase crossover in the western Mediterranean Sea. This is one of the regions selected for calibration/validation (Cal/Val) [START_REF] Ovidio | Frontiers in fine-scale in situ studies: Opportunities during the swot fast sampling phase[END_REF] in which in situ measurements have been made in the frame of SWOT [START_REF] Barceló-Llull | Pre-swot cruise report. mesoscale and sub-mesoscale vertical exchanges from multi-platform experiments and supporting modeling simulations: anticipating swot launch[END_REF]. To mitigate the computational complexity of the study and avoid the presence of continents and islands, limited

Diagnostics for evaluation

The quantitative evaluation of de-noising methods is carried out computing Root Mean Square Errors (RMSE) and Mean Spectral Ratios (MSR). RMSE for a single de-noised SWOT field h is computed as the Euclidean distance to the corresponding original, noise-free field h true :

RM SE(h) = 1 N N i=1 (h i -h true i ) 2 (4)
where N is the number of pixels and i a pixel index. Single image RMSEs are then averaged out by season for the 3 seasons considered, and are computed for SSH, |∇SSH| and ∆SSH. Thus, the test of a de-noising method with a specific set of parameters results in 9 RMSE values.

To evaluate the improvement after the application of the different de-noising techniques and parameters, we also calculate the percentage of the initial RMSE left. We calculate this RMSE residual (RM SE r ) as:

RM SE r (h) = RM SE(h) RM SE(h obs ) × 100, ( 5
)
where h is the de-noised field and h obs the original noisy field (SSH obs K or SSH obs).

The spatial spectra of the de-noised SWOT SSH are compared with the spectra of the noisefree and the noisy SWOT SSH. For each pass, we calculate the cross-track averaged, along-track power spectrum. The spectra are then averaged out over each season, leading to one spectrum per season. Information on the wavenumber spectrum calculations is given in Appendix C.

Again, to evaluate the improvement after the application of the different de-noising techniques and parameters, we compare the noise-free and de-noised fields. To do so, the Mean Spectral Ratio (MSR) is computed from the power spectral densities (PSD) of SSH. For each season, MSR is computed as:

M SR = 1 N k j=1 δk j N k j=1 log 10 P SD j (h true ) P SD j (h) 2 × δk j , (6)
where N k is the number of wavelengths considered; P SD j (h true ) and P SD j (h) are the power spectral density values at wavelength j for the original, noise-free SWOT field and the denoised SWOT field, respectively. The considered wavelengths span the interval from 9 km, the approximate effective resolution of NAtl60, to 200 km, the size of images along-track. MSR is defined above so that the best score is 0.

Exploring parameters of the de-noising methods

For all de-noising methods, a wide range of parameters are tested to identify optimal parameters according to the diagnostics presented in Section 3.2. The convolution-based methods use a single parameter that can easily be compared with the image dimensions in pixels. For the boxcar kernel, the tested parameter values go from 3 to 200 km, and correspond to the size of the box in pixels (1 km in our case). For the Gaussian kernel, the tested parameters go from 0.25 to 300 and correspond to the standard deviation, in pixels (we test up to a big sigma to have a highly oversmoothed image to reach the limit of the method). On the contrary, the geometric interpretation of the parameters of the variational method is not straightforward, and a wide exploration of the parameter space must be undertaken. However, due to computation time limitations, this cannot be performed in a strictly systematic manner. The adopted procedure is detailed below.

Orders of magnitude of the cost function terms

The orders of magnitude of the terms ∇h 2 , ∆h 2 and ∇∆h 2 composing the cost function (eq. 1) are estimated to coarsely scale the parameters λ 1 , λ 2 and λ 3 . The rationale is, for one of these terms (with its weight) to have some impact on the solution, it must be of an order of magnitude not too different from the background term m • (h -h obs ) 2 . Figure 3 shows the seasonal evolution of the derivative terms, computed from the model in a 2 • × 2 • region containing the SWOT passes used in this study. The relative ratios between ∇h 2 , ∆h 2 and ∇∆h 2 are approximately 1000:10:1. Therefore, if we want to include all three terms in the cost function, the ratios between λ 1 , λ 2 and λ 3 should coarsely be 1:100:1000. Those ratios must be only considered as a guideline to start the investigation, not a strict rule. Note that the order of magnitude of the background term after minimization of the cost function is thought to be in the range 1 to 100 in the same region. This has been estimated using the noise-free field.

Finding optimal sets of parameters

First, we created an exponential series of values to be tested for the three lambdas, consistently with the previously estimated relative ratios. For λ 1 , the series is chosen as {4 n , n = 0, ..., 7}. For λ 2 and λ 3 , the series are {10 × 4 n , n = 0, ..., 7} and {100 × 4 n , n = 0, ..., 7}, respectively. With these, six scenarios of DP de-noising are investigated, including one, two, or three penalization terms in the cost function 1. Three scenarios out of the six considered include a single penalization term (mono-parametric) of order 1, 2 and 3, successively. The other scenarios are made of terms of orders 1 and 2, 2 and 3 and the last one includes the three orders. For conciseness, particularly in the next section, we refer to the variational method with the first order term only as the λ 1 -method. We similarly refer to the λ 2 -method and to the (λ 1 + λ 2 )-method when the first two penalization terms are considered, and so on.

For each scenario, a two-step procedure is implemented to identify an optimal set of parameters. In a first step, de-noising of the full set of images is performed with all possible combinations of parameters permitted by the scenario and the parameter series defined previously. RMSEs and MSR are computed for all the combinations. In a second step, refined series of parameters are created in the neighborhood of the combination of parameters that yields the minimum RMSE and MSR scores. Image de-noising is then carried out again with all possible combinations of these series.

Optimal de-noising method

In this Section, the optimal de-noising method is searched for based on the RMSE and MSR scores described in Section 3. We investigate the KaRIn-noise-only scenario, then the all noises scenario, and finally have a closer look at the method identified as optimal. As it becomes clear in what follows, the notion of optimality does not only refer to qualitative measures.

The design of a single index summarizing the performance of the method for the different RMSEs is indeed subjective. Moreover, we take into account the ease of implementation and parameterization as a criteria in the final decision. Minimum values of RSME and MSR for each season, method, and variable are reported in tables 1 and 2 for the KaRIn-only and all noises scenarios, respectively. RMSE scores are actually expressed as the percentage of the original RMSEs, i.e. those of the original, noisy data. For SSH, |∇SSH| and ∆SSH RMSE and MSR of each de-noising configuration and parameterization, the scores do not necessarily correspond to the same optimal parameter (box size, σ or λ).

RMSE and MSR scores with KaRIn noise only

For all variables (h, ∇h, and ∆h), all seasons and all de-noising methods, minimum RMSEs are smaller in Summer than in Winter (table 1). This is expected because the oceanic surface features in Winter are smaller than in Summer [START_REF] Escudier | Eddies in the western Mediterranean Sea: characterization and understanding from satellite observations and model simulations[END_REF], so their observation is more affected by the KaRIn noise. Also, smaller structures are more affected by the smoothing due to the de-noising.

For all three seasons and all three variables, RMSEs and MSRs from the convolution-based methods and from the λ 1 -method are larger than RMSEs and MSRs from all other variational methods. Also, the λ 3 -method provides MSRs significantly higher than the other variational methods. None of these methods is the optimal de-noising one in this KaRIn-only noise configuration, and are not further discussed in the following.

In terms of both RMSEs and MSRs, and among the methods still on course, no method outperforms the others systematically and distinctly. For all three variables, RMSEs are close to each other, with differences less than a very few percents. MSRs are a bit more scattered, but without any clear predominance of a specific method in all seasons. However, the λ 2 -method exhibits the lowest MSR values in Summer, and the second lowest value in Winter, close to the

λ 2 + λ 3 -method.
Finally, this analysis persuades us to further examine the λ 2 -method for the KaRIn-only scenario (see Section 4.3). This choice is supported by the RMSE and MSR analysis above, which shows that other methods do not beat it clearly, but also by the fact that it is much easier to parametrize a single-parameter method rather than a two or three-parameter method.

RMSE and MSR scores with all noises

Normalized minimum RMSEs for h and |∇h| are higher than in the KaRIn-only scenario, by factors of 6 -12 for h, and 1.5 -4 for |∇h| (table 2). This is obviously due to the spatially correlated component of the noise (see Figure 2), which is not filtered out by any of the methods used here. Other approaches must be used to remove the correlated noise in order to obtain more accurate estimates.

Contrary to h and ∇h, RMSEs for ∆h are comparable with those obtained in the KaRInonly case. They are 5% higher only. This slight increase in RMSE is the signature of the nonlinear (quadratic, more precisely) component of the correlated error, due to the baseline dilation [START_REF] Gaultier | The challenge of using future SWOT data for oceanic field reconstruction[END_REF][START_REF] Esteban-Fernandez | Swot project: mission performance and error budget document[END_REF]. The other components are constant, linear or piecewise linear, thus are removed by the second-order derivatives.

Considering only RMSEs on ∆h, except for the boxcar and the λ 1 , no method performs significantly better than the others, and RMSEs are higher in Winter than in Summer. This is similar to the KaRIn-only scenario. The Gaussian filter performs comparatively better than in the KaRIn-only scenario.

In terms of MSRs, the methods involving λ 2 perform significantly better than the others, including the λ 3 -and the Gaussian methods. These last two exhibit MSR larger than the others by factors of 1.5 to 4. In Winter, the λ 2 -method is a little less effective than the multi-parameter methods, with a MSR twice as large.

The de-noising experiments with all noises, like those with the KaRIn noise only, lead us to favor the λ 2 -method. The reasons are similar: based on RMSEs and MSRs, the method compares favorably with the others, and a single-parameter method is much easier to parametrize.

The only result speaking against this choice is the MSR in Winter. Considering the score value though, and after the examination of the wavenumber spectra (see Figure 7), this point hardly justifies the disqualification of the λ 2 -method.

A focus on the second-order variational method

This Section investigates the sensitivity of the λ 2 -de-noising to the parameter value. Figure 4 shows the RMSEs for h, |∇h|, and ∆h, and the MSR for h as functions of λ 2 . On each graph, the three seasons are shown for both KaRIn-only (solid lines) and all noises (dashed lines) scenarios, making a total of 6 curves.

Except for h and |∇h| RMSEs in the all noises scenario, all RMSE and MSR curves exhibit a clear minimum point, which indicates the existence of an optimal, or a range of close-to-optimal λ 2 values for the de-noising. Optimal values are larger in Summer than in Winter. This is very likely because small-scale dynamics are amplified in Winter [START_REF] Sasaki | Impact of oceanic-scale interactions on the seasonal modulation of ocean dynamics by the atmosphere[END_REF]. Large λ 2 values tend to oversmooth the SSH field in Winter, leading to higher residual errors. The seasonal difference in optimal λ 2 values is particularly evident with MSR, with 100 in Winter, and 350 in Summer.

RMSEs for h and ∇h in the all noises scenario are dominated by the correlated SWOT errors, which remain present after de-noising. Consistently with the analysis of the previous section, those RMSEs are much higher in the all noises than in the KaRIn-only scenario.

In Summer, it is possible to identify a range of λ 2 values that are close-to-optimal for the four scores concomitantly. The same holds in Winter for the two scores not dominated by the correlated errors (RMSEs on ∆h and MSR). On Figure 4 suggesting that the choice of a λ 2 is inevitably subject to a part of subjectivity if no additional information on the ocean surface dynamics is available. Not detailed here, the results from the other (single or multiple-parameter) configurations of variational de-noising also exhibit such overlaps of intervals, except for λ 1 . In all cases, de-noising leads to correct orders of magnitude for all fields, and particularly for |∇h| and ∆h. As expected and already shown by [START_REF] Gómez-Navarro | Swot spatial scales in the western mediterranean sea derived from pseudo-observations and an ad hoc filtering[END_REF][START_REF] Chelton | Prospects for future satellite estimation of small-scale variability of ocean surface velocity and vorticity[END_REF], the original SWOT data affected by small-scale noise does not provide any useful information about SSH derivatives. The denoising method corrects this efficiently, and makes it possible to identify the main structural characteristics of the fields.

A strong spatially-correlated noise shows strong signatures on h, moderate signatures on |∇h|, and low signatures on ∆h, except at the outer boundaries of the swath. The low signature on ∆h was already observed in the RMSEs, and is due to the specific spatial structure of this noise. Most components are linear in the across-track direction. In the along-track direction, error correlations are high (Figure 2). Therefore, the correlated noise has a low effect on the second-order derivatives. The remaining noise at the outer boundaries is due to the finite difference method used to compute the derivatives described in Appendix A.

Although the resulting fields of ∆h fall in correct orders of magnitude and capture the structure of the true fields at the scale of the swath, the kilometric-scale fronts and filaments are smoothed out by the de-noising. Solving this issue would require the development of more sophisticated de-noising techniques, or a post-processing of the present result including, for example, some ocean dynamics through data assimilation techniques. This will be a natural step forward, since the first motivation for developing a de-noising technique constraining ∆h is precisely the combined assimilation of h and its first two derivatives, as stated in the introduction. From this spectral viewpoint, the de-noised data matches the noise-free data well at all scales down to ∼15 km. In the noisy data, the noise amplitude approaches the signal amplitude at wavelengths of 50 km in Summer and 40 km in Winter, and dominates the signal at shorter wavelengths. This is efficiently corrected by the de-noising. The process seems more efficient in Winter than in Summer, probably because of higher PSDs in Winter related to more intense ocean surface processes.

Following the definition proposed by Wang et al. [START_REF] Wang | On the spatial scales to be resolved by the surface water and ocean topography ka-band radar interferometer[END_REF] for the spatial scale resolved by SWOT, the de-noising reduces this scale by a factor of 2, leading to resolved scales of between approximately 20 and 30 km. Wang et al. [START_REF] Wang | On the spatial scales to be resolved by the surface water and ocean topography ka-band radar interferometer[END_REF] define the spatial scale resolved by SWOT by the wavelength at which the SWOT noise spectrum intersects the spectrum of the true signal (SSH model here). Figure 7 indicates resolved scales of 50, 40, and 50 km in the JAS12, FMA13 and JAS13 scenarios respectively, in both just KaRIn and all noises cases. After de-noising, the resolved scales are reduced to 25, 20, 20 km in the KaRIn-only case and 30, 20, 30 in the all noises case. Even below these scales, the noise left is very low, and within the variability of SSH model (red envelope in Figure 7). At wavelengths near 10 km, the noise is reduced by 10 4

in the JAS scenarios and 10 3 in the FMA scenario.

Discussion and Conclusions

Several objectives of the SWOT mission will be met only if the small-scale noise affecting the data can be efficiently removed. Small-scale noise, in particular the spatially uncorrelated KaRIn instrument noise, prevents the computation of horizontal SSH derivatives. This limits both the direct estimation of relevant oceanic variables on the SWOT swath, and the use of SWOT data to build gridded products of altimetry.

To remove the small-scale SWOT noise, we propose a de-noising method that performs better than conventional convolution-based methods both in terms of RMSE (physical space diagnostic) and spectra. The method, which originates from image processing applications, is based on the regularization of the SWOT SSH data by the penalization of its derivatives of orders 1 to 3 in a variational, optimization framework. This approach is chosen because it is in close connection with the oceanic variables of interest, namely geostrophic velocity and vorticity. After a thorough evaluation based on a large number of simulated SWOT scenes, the variational de-noising method exhibits better performance than standard, boxcar and Gaussian filters. We find the method performs best when only the second-order derivative (λ 2 ) is considered in the cost function. Only one parameter needs to be set, which makes the parameterization of the method as simple as a convolution-based method. We find that this parameter can be set smaller or larger in function of the characteristics of our field: the higher the intensity of the signal, the lower the penalization and thus the value of the parameter (as we find in the FMA13 λ 2 values in contrast to JAS12). Also, if the noise level in our fields is higher (all noises scenario), the more we need to penalize and the larger the parameter value. In other words, the higher the signal to noise ratio (SNR) the less we need to penalize our field, and so the smaller λ 2 .

The method will require further investigations before operational applications, since we have focused our attention to one particular region (the Western Mediterranean Sea), with an ocean circulation free of tidal forcing, and a prescribed Significant Wave Height (SWH) of 2 m. The present study shows that in one single region, the range of optimal parameters changes with the season, due to seasonal changes in the ocean surface dynamics. Similar conclusions are certainly expected with respect to regional and dynamical regimes. The NATL60 simulation used here does not include tidal forcing. The behavior and efficiency of the de-noising method may be questioned in presence of tidal motions and particularly tide-generated internal waves.

Finally, the SWH prescribed in the SWOT simulator to compute the KaRIn error amplitude is prescribed to 2 m. As the SWH varies geographically and according to the atmospheric regime, KaRIn errors smaller or larger than those computed for the present study with the SWOT simulator can be expected [START_REF] Wang | On the spatial scales to be resolved by the surface water and ocean topography ka-band radar interferometer[END_REF]. The first two aspects (geographic variations of ocean dynamics and internal tides) are presently under study using data from several high-resolution simulations that include tidal forcing: the HYbrid Coordinate Ocean Model (HYCOM) [START_REF] Chassignet | US GODAE: global ocean prediction with the HYbrid Coordinate Ocean Model (HYCOM)[END_REF],

the Massachusetts Institute of Technology general circulation model (MITgcm) [START_REF] Marshall | A finite-volume, incompressible navier stokes model for studies of the ocean on parallel computers[END_REF], and the recent extended NATL60 -eNATL60-simulation, not yet published).

The method should also benefit from additional developments to reconstruct more realistic fields of relative vorticity on the SWOT swath, and could ultimately lead to the estimation of vertical velocities. The de-noising process inevitably smoothes out the very fine-scale, elongated structures usually visible in surface relative vorticity fields [10, e.g.]. Restoring these structures should be investigated, perhaps using appropriate image processing techniques [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF][START_REF] Yan | Nonlocal hierarchical dictionary learning using wavelets for image denoising[END_REF], or methods already developed in the oceanographic community such as Lagrangian advection [START_REF] Rogé | Altimetric lagrangian advection to reconstruct pacific ocean fine-scale surface tracer fields[END_REF][START_REF] Berti | Lagrangian reconstructions of temperature and velocity in a model of surface ocean turbulence[END_REF].

Dynamical models could also be used in a data assimilation framework.

To conclude, the de-noising method opens the way to several relevant applications using the SWOT data, possibly including SWOT data validation, assimilation, and SSH mapping. We mention SWOT data validation due to the in-painting capability of the variational de-noising method, i.e. the fact that the process naturally fills the 20-km gap of the SWOT swath (the gap is in-painted, and emptied again after de-noising to restore SWOT data in the original shape). In other words, the SWOT KaRIn data are interpolated on the track of the SWOT nadir altimeter. This is obviously relevant for data comparison and validation. De-noising is also interesting to pre-process the SWOT data before their assimilation in ocean circulation models. This actually was a primary motivation for the method development. Computing spatial derivatives of the SWOT data allows the implementation of data assimilation methods that account for SWOT error correlations [START_REF] Ruggiero | An efficient way to account for observation error correlations in the assimilation of data from the future SWOT High-Resolution altimeter mission[END_REF][START_REF] Yaremchuk | On the approximation of the inverse error covariances of high resolution satellite altimetry data[END_REF]. Alternatively, the relative vorticity derived from the de-noising can be directly assimilated. This option has not been explored yet to our knowledge. Finally, the de-noising can also be combined with other techniques to improve the assimilation. We particularly think about the technique recently developed by Metref et al. [START_REF] Metref | Reduction of spatially structured errors in wide-swath altimetric satellite data using data assimilation[END_REF] to significantly reduce the impact of the geometrically structured, highly correlated SWOT errors (roll, phase, timing, and baseline errors).

  (h obs ) to process and compare with the truth. A set of 543, 121×200 km 2 pseudo-SWOT scenes are generated in the western Mediterranean Sea, covering one Winter and two Summer seasons (choice limited by the model's time span). SWOT scenes are sampled from the fast-sampling phase satellite orbit, focusing on a cross-over region, i.e., where an ascending pass crosses a descending pass, therefore providing 2 passes per day. The SWOT data simulation is carried out over three 3-month periods: July to September 2012 and 2013 (JAS12 and JAS13 hereafter), representing the Summer season, and February to April 2013 (FMA13) representing the Winter season. Summer periods provide 92 (resp. 91) of ascending (resp. descending) passes; the Winter period provides 89 (resp.

  subregions of the SWOT swaths are sampled. These subregions are 121 km-wide (the width of 2 SWOT swaths plus the gap) and 200 km-long. The region, the SWOT passes and the subregions are shown in Figure1. It is worth noting that each scene is affected by a unique realization of the SWOT error.In this work, image de-noising techniques are first applied to the pseudo-SWOT scenes affected by the KaRIn noise only (SSH obs K), then to the scenes containing all errors (SSH obs). This approach allows to discriminate the effects and the performance of image de-noising in presence of the spatially correlated SWOT errors. A few realizations of the different components of the SWOT error are shown on Figure2, where we can observe how most errors exhibit strong and long-range correlations, whilst the KaRIn error does not show any correlation at all.

Figure 1 :

 1 Figure 1: SSH model outputs [m] for cycle 1 of pass 9 (left) and 22 (right) of the JAS12 dataset. In red the subregions selected.

Figure 2 :

 2 Figure 2: Examples of noises and errors [m] added by this SWOT simulator version 2.21 to our study region fast-sampling phase for JAS12 pass 9 (a) and 22 (b). Note that these simulations are performed without the 20-km gap at nadir.

Figure 3 :

 3 Figure 3: Seasonal variations of the cost function terms ∇h 2 , ∆h 2 and ∇∆h 2 , from top to bottom. Shaded areas indicate the JAS12, FMA13 and JAS13 periods from left to right, respectively. The mean and median values are printed for each period and for the whole year (upper right corners, in bold).

, horizontal error bars indicate the range of λ 2

 2 values that provide scores higher than the minimum by less than 5%. Subjectively based on this information, we propose λ 2 intervals of [300 -400] in Summer and [100 -120] in Winter. MSR results for the two Summer seasons indicate slightly different optimal values,

Figure 4 :

 4 Figure 4: Scores of RMSE and MSR of λ 2 -method from just KaRIn (solid line) and all noises (dashed line) for all 3 seasons. Horizontal error bars in the RMSE plots show the the range of λ 2 values that provide scores higher than the minimum RMSE by less than 5%.

Figure 7

 7 Figure 7 shows power spectral densities (PSD) of h. The rows distinguish the just KaRIn noise added and the all noises cases. The columns are for Summer 2012, Winter 2013, and Summer 2013. On each graph, the spectra are shown for the noise-free data (SSH model),the noisy data (SSH obs), the de-noised data (SSH obs f), the pre-de-noising noise (noise) and post-de-noising noise (noise f). The de-noised data have been obtained with the λ 2 -method with parameter values chosen in the intervals identified in Section 4, and indicated on each graph.

Figure 5 :

 5 Figure 5: Fields of pass 09, cycle 2 (a) and 6 (b) of JAS12 dataset compared to the fields filtered with λ 2 = 430. From top to bottom: SSH, gradient of SSH and laplacian of SSH. From left to right: model interpolated to SWOT grid (SSH model), SSH model + KaRIn noise (SSH obs K), SSH model + all noises (SSH obs), filtered SSH obs K (SSH obs K f) and filtered SSH obs (SSH obs f).

Figure 6 :

 6 Figure 6: Fields of pass 09, cycle 2 (a) and 6 (b) of FMA13 dataset compared to the fields filtered with λ 2 = 95. From top to bottom: SSH, gradient of SSH and laplacian of SSH. From left to right: model interpolated to SWOT grid (SSH model), SSH model + KaRIn noise (SSH obs K), SSH model + all noises (SSH obs), filtered SSH obs K (SSH obs K f) and filtered SSH obs (SSH obs f).

Figure 7 :

 7 Figure 7: Spatial spectra of the model interpolated data (SSH model) are shown in red and of the pseudo-SWOT data (SSH obs) in black. Blue lines indicate the filtered pseudo-SWOT spectra (SSH obs f) obtained with the optimal λ 2 found with the MSR score. The dashed lines are the noise spectra of SSH obs (noise) and SSH obs f (noise f). Shaded areas show values between the 5th and 95th percentiles, showing the PSD variability. Top row shows pseudo-SWOT data with just KaRIn noise added and botom row with all noises. Columns represent the different seasonal datasets from left to right: Summer 2012 (JAS12), Winter 2013 (FMA13) and Summer 2013 (JAS13).

Figure D. 8 :

 8 Figure D.8: Fields of pass 09, cycle 2 (a) and 6 (b) of JAS12, all noises dataset. From left to right: comparison between the SSH model, SSH obs, SSH obs filtered with our approach and λ 2 = 430 (SSH obs f), with the optimal boxcar (SSH obs f bc) and with the optimal Gaussian (SSH obs f ga) methods. From top to bottom: SSH, gradient of SSH and laplacian of SSH

Figure D. 9 :

 9 Figure D.9: Fields of pass 09, cycle 2 (a) and 6 (b) of FMA13, all noises dataset. From left to right: comparison between the SSH model, SSH obs, SSH obs filtered with our approach and λ 2 = 95 (SSH obs f), with the optimal boxcar (SSH obs f bc) and with the optimal Gaussian (SSH obs f ga) methods. From top to bottom: SSH, gradient of SSH and laplacian of SSH

Table 1 :

 1 Scores summary of the different de-noising methods for the just KaRIn dataset (h=SSH).

	Season	De-noising		RMSE r		Minimum MSR
			Method	SSH	|∇SSH|	∆SSH
			Boxcar	12.43	0.094	0.300	0.2010
			Gaussian	12.88	0.084	0.251	0.1111
	JAS12		1 2	12.55 08.71	0.084 0.050	0.279 0.247	0.2028 0.0143
		DP	3	09.06	0.051	0.247	0.1021
			1 + 2	08.72	0.050	0.247	0.0192
			2 + 3	08.68	0.049	0.247	0.0205
			1 + 2 + 3	08.66	0.049	0.246	0.0259
			Boxcar	15.04	0.177	0.511	0.1066
			Gaussian	13.60	0.153	0.437	0.0746
	FMA13		1 2	15.41 10.92	0.173 0.115	0.483 0.420	0.1498 0.0178
		DP	3	10.86	0.113	0.416	0.0682
			1 + 2	10.92	0.115	0.420	0.0208
			2 + 3	10.79	0.113	0.416	0.0168
			1 + 2 + 3	10.82	0.113	0.416	0.0255
			Boxcar	11.98	0.086	0.326	0.1796
			Gaussian	12.81	0.076	0.277	0.0911
	JAS13		1 2	12.78 08.96	0.083 0.053	0.309 0.274	0.2031 0.0216
		DP	3	09.11	0.053	0.273	0.1010
			1 + 2	08.97	0.053	0.274	0.0394
			2 + 3	08.84	0.052	0.272	0.0243
			1 + 2 + 3	08.84	0.052	0.272	0.0269

Table 2 :

 2 Scores summary of the different de-noising methods for the all noises dataset (h=SSH).

	Season	De-noising		RMSE r		Minimum MSR
			Method	SSH	|∇SSH|	∆SSH
			Boxcar	90.31	0.171	0.303	0.2024
			Gaussian	90.10	0.156	0.264	0.1181
	JAS12		1 2	87.60 90.57	0.159 0.174	0.281 0.261	0.1922 0.0307
		DP	3	90.22	0.156	0.265	0.1359
			1 + 2	87.61	0.158	0.262	0.0328
			2 + 3	90.22	0.156	0.261	0.0391
			1 + 2 + 3	87.40	0.156	0.261	0.0395
			Boxcar	90.88	0.250	0.511	0.1274
			Gaussian	90.76	0.221	0.435	0.0515
	FMA13		1 2	89.89 91.11	0.237 0.226	0.484 0.432	0.1415 0.0314
		DP	3	90.96	0.226	0.432	0.0868
			1 + 2	89.90	0.223	0.435	0.0160
			2 + 3	91.00	0.226	0.430	0.0177
			1 + 2 + 3	89.82	0.220	0.430	0.0203
			Boxcar	89.73	0.137	0.328	0.1792
			Gaussian	89.18	0.126	0.289	0.1152
	JAS13		1 2	84.30 90.36	0.131 0.142	0.310 0.287	0.1895 0.0254
		DP	3	89.66	0.127	0.290	0.1251
			1 + 2	84.19	0.131	0.287	0.0237
			2 + 3	89.66	0.127	0.286	0.0285
			1 + 2 + 3	83.75	0.127	0.286	0.0267

Appendix A. Calculation of Laplacian Laplacian are computed using finite differences, following the method proposed by [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. We note h the image of size N x × N y . In a first step, the two components of the gradient are computed as (i = 1, ..., N x ; j = 1, ..., N y ):

In a second step, Laplacian is computed as the divergence of the gradient. Divergence of vector a = (a x , a y ) is computed as:

where:

The scheme implemented at the boundaries preserves the image size, contrary to what a standard five-point stencil Laplacian operator would do. Preservation of image size is essential in the gradient descent iterations to end up with a final image of size similar to the initial image.

Appendix B. FISTA

To speed up the gradient descent iterations, the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] is implemented. Setting t 0 = 1 and introducing an auxiliary variable y initialized as y 0 = h 0 , the iterative algorithm of Eq. 3 becomes:

Appendix C. Calculation of spatial spectra

The spatial spectra used as one of the scores for the de-noising parameterizations are calculated as follows:

1. Apply a linear detrending;

2. Remove the spatial mean;

3. Apply a Tukey window with a 0.5 fraction of the window inside the cosine tapered region;

4. Compute the 1D spatial Fourier power spectra along-track for each SSH swath across-track dimension.

Appendix D. Qualitative figures of different methods

To better illustrate the advantage of our de-noising approach, we show in Figures D.8

and D.9 the fields provided by the boxcar and Gaussian methods, corresponding to the λ 2 experiments presented in Figures 5 and6. We only show the all noises scenario. Boxcar derivatives fields are very noisy, as it is specially visible for the laplacian fields. With the Gaussian method, the laplacian is less noisy than with our method, but the gradient is oversmoothed.

Appendix E. Softwares

• Standard image techniques: For both boxcar and Gaussian kernel python's scipy.ndimage module was used with the following specific functions:

-Boxcar filter: scipy.ndimage.generic f ilter() -Gaussian: scipy.ndimage.gaussian()

• Variational regularization method: https://github.com/LauraGomezNavarro/SWOTmodule