Nine years of SMOS Sea Surface Salinity global maps at the Barcelona Expert Center

After more than 10 years in orbit, the Soil Moisture and Ocean Salinity (SMOS) European mission is still a unique, highquality instrument for providing Soil Moisture over land and Sea Surface Salinity (SSS) over the oceans. At the Barcelona Expert Center (BEC), a new reprocessing of 9 years (2011-2019) of global SMOS SSS maps has been generated. This work presents the algorithms used in the generation of BEC global SMOS SSS product v2.0, as well as an extensive quality as5 sessment. Three SMOS SSS fields are distributed: a high-resolution level 3 product (with doi: http://dx.doi.org/10.20350/ digitalCSIC/12601 (Olmedo et al., 2020a)) consisting of a binned SSS in 9-day maps at 0.25× 0.25; a low-resolution level 3 SSS computed from the binned salinity by applying a smoothening spatial window of 50-km radius; and a level 4 SSS (with doi: http://dx.doi.org/10.20350/digitalCSIC/12600 (Olmedo et al., 2020b)) consisting of daily, 0.05×0.05o maps that are computed by multifractal fusion with Sea Surface Temperature maps. For the validation of BEC SSS products, we have applied a 10 battery of tests aiming at the assessment of quality of the products both in value and in structure. First, we have compared BEC SSS products with near-to-surface salinity measurements provided by Argo floats. Secondly, we have assessed the geophysical consistency of the products characterized by singularity analysis, and also the effective spatial resolutions are estimated by means of Power Density Spectra and Singularity Density Spectra. Finally, we have calculated full maps of SSS errors by using Correlated Triple Collocation. We have compared the performance of BEC SMOS product with other satellite SSS and reanal15 ysis products. The main outcomes of this quality assessment are: i) the bias between BEC SMOS and Argo salinity is lower than 0.02 psu at global scale, while the standard deviation of their difference is lower than 0.34 and 0.27 psu for the high and low resolution level 3 fields (respectively) and 0.24 psu for the level 4 salinity; ii) the effective spatial resolution is around 40 km for all SSS products and regions; and iii) BEC SMOS level 4 product is globally the one with the lowest salinity error, while BEC SMOS low-resolution level 3 more accurate in regions strongly affected by rainfall and continental freshwater discharges. 20 Copyright statement. TEXT 1 https://doi.org/10.5194/essd-2020-232 O pe n A cc es s Earth System Science Data D icu ssio n s Preprint. Discussion started: 2 October 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
The European Space Agency (ESA) Soil Moisture and Ocean Salinity (SMOS) satellite was launched in November 2009, carrying the first orbiting radiometer that collects regular and global observations from space of two Essential Climate Variables (ECV) according to the Global Climate Observing System: Sea Surface Salinity (SSS) and Soil Moisture (SM) (Font et al.,25 2010; Kerr et al., 2010;Mecklenburg et al., 2009). After more than 10 years in orbit, the SMOS mission has been a success in terms of both technology and science, providing SSS and SM data derived from the SMOS measurements (Reul et al., 2020; Ker, 2016) (and references therein).
The Barcelona Expert Centre (BEC) was created in 2007 to support the Spanish contribution to SMOS mission activities.
Since the beginning, BEC's goals have been to contribute to the quality assessment and the development of algorithms for the 30 retrieval of geophysical variables from SMOS data as an ESA Level 2 Ocean Salinity Expert Support Laboratory; and to the calibration and validation activities as a Level 1 Expert Support Laboratory. In the recent years, BEC has developed a SMOS SSS internal processing chain that generates SSS maps from SMOS raw data (Level 0) to Level 3 and 4 (L3 and L4 added-value SSS maps), thus allowing the integration of improvements in the different levels of the processing. The resulting products are freely distributed through a sftp service (http://bec.icm.csic.es/bec-ftp-service/). 35 In this work we present the new reprocessing of the BEC SMOS SSS global L3 and L4 products v2.0 for a 9-years period comprising 2011 to 2019, that comes with an improvement of the currently used methodology. This new reprocessing is focused on four aspects: -Improving salinity gradients: A new filtering criterion that is more geophysically consistent has been introduced.
-Improving the latitudinal and seasonal biases: An empirical correction to reduce the latitudinal and seasonal biases that 40 affected the previous version of the product (Olmedo et al., 2019b) has been applied.
-Improving the quality of the acquisitions close to the coast: The interpolation scheme and also the level 4 fusion techniques have been adapted to preserve small-scale gradients close to the coast.
-Providing an estimate of the sea surface salinity uncertainty: An explicit expression to propagate the errors from brightness temperature uncertainties to the final SSS product have been introduced. 45 To assess the performance of the BEC SMOS SSS product v2.0, the complete 9-year time series of SSS maps is first compared with the salinity measurements provided by Argo. Secondly, an extensive battery of validation methods is applied to one year (2017) of data, and the results are compared with three other satellite and one reanalysis SSS products. Those methods are: i) statistics of the differences with Argo salinity match ups; ii) singularity analysis to assess the geophysical consistency of the data (Turiel et al. (2008b)); iii) spectral analysis to analyse the effective spatial resolution of each product, temperature, rain rate, wave model, 10-meter wind speed, 10-meter neutral equivalent wind (zonal and meridional components), significant height of wind waves, 2-meter air temperature, surface pressure, and vertically integrated total water vapour (Zine et al., 2007).
We have used as multiyear salinity reference the annual climatological salinity value provided by the World Ocean Atlas 2013 (WOA2013) at 0.25 • × 0.25 • (Zweng et al., 2013). The SSS provided by WOA2013 is taken as the reference value to be added to SMOS salinity anomalies (see section 2.2.1). We use the average decadal product, which is accessible at Na-90 tional Oceanographic Data Center (https://www.nodc.noaa.gov/cgi-bin/OC5/woa13/woa13.pl). We have also used the monthly climatology at 0.25 • × 0.25 • provided by WOA2013 for the correction of latitudinal and seasonal biases (see section 2.2.5).

Retrieval of SMOS debiased Sea Surface Salinity Anomalies
The debiased non-Bayesian (DNB) retrieval approach proposed in Olmedo et al. (2017) has been used to retrieve the 9 years 95 (2011-2019) of SMOS SSS maps. This methodology consists of retrieving a single value of SSS from each First Stokes Brightness Temperature (TB) measurement, that we refer to as raw SSS retrieval. These raw SSS are then appropriately classified, filtered and combined to build global SSS maps.
For the retrieval of raw SSS, the difference between the First Stokes TB measured and modelled is optimized as a function of the salinity value. A geophysical forward model links the modeled TB to the SSS. Besides the dielectric constant model 100 proposed by Klein and Swift (1977) for TB coming out flat sea, the forward model accounts for the contribution to TB of the sea surface roughness (Guimbard et al., 2012), the reflected emission of the atmosphere, the reflection on sea surface of the galactic emission (Tenerelli et al., 2008) and the sun glint (Reul et al., 2007). All these contributions are taken into account in the salinity retrieval.
All the raw salinity retrievals over 9 years (s raw n with n = 1, . . . , N where N is the total number of raw retrievals in 9 years) 105 are classified as a function of the satellite overpass direction (d), latitude (ϕ), longitude (λ), across-track distance (x), and incidence angle (θ). The underlying hypothesis of this approach is that the systematic errors (i.e., those which are independent of time) are the same for all the s raw n that are acquired under each fixed condition γ = (ϕ, λ, d, x, θ). Therefore, the systematic errors are the same for all the retrievals in the set {s raw n (γ)} with n = 1, . . . N γ and N γ the number of retrievals during 9 years with that specific value of γ.
Finally, the debiased salinity value for the acquisition conditions γ, s n (γ), is computed by adding an external multiyear salinity reference to the s n (γ); in this case, we have used the annual salinity field provided by WOA2013. 120 The retrieval algorithm proposed above effectively removes local biases, especially those produced by the land-sea contamination and artifacts produced by permanent Radio Frequency Interference (RFI) sources.

Estimation of SSS error
Each value of raw salinity s raw n can be associated a retrieval error which is computed according to the next equation: where the σ H n and σ V n are the radiometric sensitivities for the horizontal (H) and vertical (V) polarizations of the brightness temperature, respectively (that are contained in the ESA L1B product), and ∂In ∂s is the derivative of the modelled First Stokes divided by 2 (I n ) with respect to the salinity (that can be estimated numerically).

Filtering criteria
Filtering out degraded measurements in the generation of the SMOS SSS maps is a key aspect. Without applying any filter, 130 the error may become too large for many scientific applications; on the other hand, when the filtering criteria are too strict, the coverage of maps may be dramatically decreased and part of the geophysical variability may be lost. In Olmedo et al. (2017), filtering criteria based on the statistical properties of {s raw n (γ)} were proposed and the resulting maps led to an almost complete coverage with an acceptable salinity accuracy (see Olmedo et al. (2017) for more details). We revisit these filtering criteria in order to decrease the error of the retrieved salinity and improve the description of salinity gradients in highly dynamic 135 regions.
We apply the following filtering criteria: -Basic filtering: Any s raw n (γ) out of the range of [0, 50] psu is not considered as part of the corresponding set of valid {s raw n (γ)}.
-Discarding some full sets of {s raw n (γ)}: For a given value of γ, we consider a particular set of {s raw n (γ)} valid only 140 when: -It contains more than 100 salinity retrievals; -The standard deviation of its distribution is lower than 10 psu; -The absolute value of the skewness of the distribution is lower than 1; and -The kurtosis of the distribution is greater than 2. Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License.
-Outlier criteria: We discard specific salinity retrievals s raw n (γ) when the corresponding SMOS debiased salinity anomaly (s n (γ)) satisfies: where σ 2 ϕ,λ is the geophysical variance of the salinity expected at the gridpoint (ϕ, λ). In order to estimate σ 2 ϕ,λ we use SMOS 9-day salinity maps that are computed from a more relaxed filtering criteria, which is: Notice that σ γ is always greater than σ ϕ,λ , because σ γ contains the variability corresponding to the salinity uncertainty and the salinity geophysical variability (see Olmedo et al. (2019a)).
-Temporal and geophysical consistency: We temporally and spatially collocate all the debiased retrievals s n (γ) in 9-day maps with the fixed grid at 0.25 o × 0.25 o . The resulting collocated set of SSS is denoted as {s(t T , ϕ, λ)}, being t T all 155 the acquisition times in the 9-day period which is indexed by the time T (typically, T corresponds to the central day).
In particular, for a given geographical location (ϕ, λ), we combine all the different values of SSS under all acquisition conditions at that specific location and 9-day period, what means combining all the satellite overpass directions (d), across-track distances (x), and incidence angles (θ) that happen at all the time t T in that period. Then, we consider as valid salinity measurements only those satisfying: withs(T, ϕ, λ) and σ(T, ϕ, λ) being the mean and standard deviation of the set {s(t T , ϕ, λ)}, that is, all values of SSS at longitude ϕ, latitude λ and 9-day period centered around T . Finally, we average the valid salinity values of {s(t T , ϕ, λ)} in that period to obtain the binned 9-day map at 0.25 o × 0.25 o , s 0 (T, ϕ, λ).

Mitigation of temporal biases
165 SMOS measurements are affected by biases that depend on time (see Martín-Neira et al. (2016)). The methodology described in 2.2.1 aims at removing the systematic biases affecting SMOS measurements, i.e., those biases that depend on the acquisition conditions (γ) but not on time. To address the temporal biases, we follow the approach proposed in Olmedo et al. (2017), which consists of assuming that the global average of SSS does not change with time. Figure A1 shows the difference between the spatial averaged salinity value of s 0 (T, ϕ, λ) and the spatial averaged salinity value of the annual reference (WOA13). We 170 assume that for a given time T this difference has to be zero. Therefore, we correct each map with this difference, imposing the spatial averaged salinity value of every global map to be equal to the annual reference. We notate the temporal corrected binned salinity fields as s 1 (T, ϕ, λ). Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License.

Correction of latitudinal-seasonal biases
The corrections applied so far aim at systematic biases which are time-independent or space-independent, and therefore can be corrected separately. However, after applying both corrections, residual biases depending at the same time on time and on the geographical position are still present (see Figure A2). These latitudinal-seasonal biases are known to happen also in other L-band satellite missions and are supposed to be due to the different direct influence of the Sun on the instrument along its trajectory depending on the season of the year. We have therefore applied the latitudinal-seasonal bias correction proposed in Olmedo et al. (2019b), that is computed as follows:

180
-We compute SMOS monthly climatologiess 1 (m, ϕ, λ) for each month m of the year by averaging all the s 1 (T, ϕ, λ) where T belongs to the same month m of the processed nine years. Recall that this climatology is defined on a 0.25 o × 0.25 o grid, as this is the grid for s 1 (T, ϕ, λ).
-For each month m, we subtract the WOA2013 monthly climatology, s W OA (m, ϕ, λ), from the corresponding SMOS monthly climatology: -We fit ∆s 1 (m, ϕ, λ) by a second degree polynomial of the latitude. That is, for every month, m, and every value of latitude in the 0.25 o grid, ϕ, we compute the polynomial p(m, ϕ) that minimizes the following cost function: -After computing the optimal polynomials p(m, ϕ), we correct the maps s 1 (t, ϕ, λ) by daily interpolating the polynomial p(m, ϕ) to the specific moment of the month. We denote the latitudinal-seasonal debiased SSS by s 2 (T, ϕ, λ).

Mitigation of residual spatial biases
After applying all the above corrections, we make the last check. By construction, at each geographical location the average 195 salinity of the full period should be equal to the multiyear reference introduced in section 2.2.1. We have found significant differences between both averages that may be due to an inaccurate determination of the SMOS-based climatology. This may happen when the distributions of the values {s raw n (γ)} for a given γ significantly deviates from a Gaussian distribution (especially, if it is slightly skewed), and then differences between the modes (used in the computation of s c (γ)) and the means of all s n (γ) after applying the filtering criteria (section 2.2.3) are significant. In order to mitigate this last bias, we remove the 200 map corresponding to the difference between the mean average of all the SMOS SSS in the 2011-2019 period (s 2 (T, ϕ, λ)) and WOA2013 (see Figure A3). The resulting salinity field is our L3 high-resolution product s L3 (T, ϕ, λ). In future versions of this product we will introduce a better definition of the SMOS-based climatology to avoid this last correction step.

Multifractal fusion techniques
We use the multifractal fusion techniques introduced in Umbert et al. (2014) to increase the spatial and temporal resolutions 205 of SMOS SSS L3 maps (Olmedo et al., 2016). Multifractal fusion methods are based on the hypothesis that different ocean scalars have the same singularity exponents (SE). From a mathematical point of view, the SE of a function at a given point is a measure of the local regularity of the function at that point (Turiel et al., 2008a). It has been shown that synoptic maps of different ocean scalars show the same multifractal structure due to the effect of geophysical turbulence. Starting from the SE extracted from SST maps (Turiel et al., 2005(Turiel et al., , 2008b, it has been observed that other scalars, such as chlorophyll concentration 210 maps (Nieves et al., 2007;Umbert et al., 2020) and even brightness temperature maps at given frequencies (Isern-Fontanet et al., 2007) present the same structure and even values of SE. This correspondence can be used to improve the quality of the SMOS SSS L3 maps by using as a template OSTIA SST, which is an ocean scalar measured with better spatio-temporal resolution and quality than SSS. Assuming that both variables have the same SE, it can be shown (Umbert et al., 2014) that as a first order approximation the following local relationship holds: where a and b are smooth functions, that is they must have small gradients, as otherwise they would introduce additional SE.
The estimation of the smooth functions a and b is done by means of a local weighting average (see Umbert et al. (2014) and Olmedo et al. (2016) for more details). Taking advantage of the fact that a and b do not have sharp variations over large regions, the evaluation of a and b is performed by locally-weighted linear regression. We have employed a similar local-220 weighting function as in Olmedo et al. (2016), that is, the inverse of the 4-th power of the distance to the central point. In order to better describe the small-scale features, the local weighting is limited to points at most at a distance of R = 2.5 • from the central point.
By means of this multifractal fusion method, SMOS L4 SSS maps with the same spatial and temporal resolutions as the template (OSTIA SST), i.e., daily maps at a spatial grid of 0.05 • × 0.05 • , are obtained. The product is distributed in netCDF files and it contains two different salinity fields and one estimation of the SSS uncertainty.
The two SSS fields are denoted as: -BEC SMOS HR SSS product: where HR stands for High Resolution and contains the binned salinity field s L3 (T, ϕ, λ).

230
-BEC SMOS LR SSS product: where LR stands for Low Resolution and is a low pass filtered version of s L3 (T, ϕ, λ), computed by applying a spatial radio of 50km. This product is denoted by s L3 low (T, ϕ, λ).
The BEC SMOS L4 SSS product v2.0 (hereafter BEC L4) consists of daily SSS maps at a regular grid of 0.05 o × 0.05 o .
This product is denoted by s L4 (T, ϕ, λ). Passive (SMAP) mission are available (this mission has been operating since early 2015 (Entekhabi et al., 2010)). The satellite SSS products used for the intercomparison are: -CATDS SMOS products: 9-day SMOS SSS maps provided by Centre Aval de Traitement des Données SMOS (CATDS).
This product decreases the mean bias over the open ocean with respect to previous versions and it improves ice filtering, 245 which leads to an improvement of SSS at high latitudes, especially in the Southern Ocean (Boutin et al., , 2018Kolodziejczyk et al., 2016). -REMSS SMAP products: 8-day running Remote Sensing Systems SMAP Level 3 Sea Surface Salinity Standard Mapped Image version v4 which is freely available at www.remss.com/missions/smap. In particular, we have used the field sss_smap, which is a smoothened measurement at approximately 70km resolution. The major change in Version 4.0 255 from Version 3.0 is an improved land correction, which allows for SMAP salinity retrievals closer to the coast (Meissner et al., 2018).

In situ salinity: Argo floats
For the purpose of direct comparison of values, we have used in situ salinity data obtained by Argo profilers. We consider the uppermost Argo salinity between 5 and 10 m depth (hereafter Argo SSS). Argo data is collected and made freely available by Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License.

Reanalysis Sea Surface Salinity
We are also interested in analyzing the strengths and weaknesses of the satellite products when compared with a reanalysis product. To this end and for completeness, we have used the SSS fields provided by ARMOR3D Near Real-Time (Nardelli,265 2012; Nardelli et al., 2016;Droghei et al., 2016) corresponding to the year 2017. We use version the 4 of ARMOR3D which is freely available in the Copernicus Marine Environment Monitoring Service (CMEMS) service desk (https://resources.marine. copernicus.eu/?option=com_csw&task=results?option=com_csw&view=details&product_id=MULTIOBS_GLO_PHY_REP_ 015_002). In the generation of the SSS fields provided by ARMOR3D (hereafter CMEMS SSS), a correction based on the ISAS-CORA SSS field is applied as well as a combination of a Quality Control SSS measurements obtained from ISAS-270 CORA (both distributed through CMEMS) and a high-pass filter of Reynolds SST L4 satellite observations.

Sea Surface Temperature
OSTIA SST is used as a reference to assess the spatial structure, geophysical consistency, and the effective resolution of the SSS satellite products (see section 2.1.2 for the complete description).

Comparison with Argo
Assuming that Argo values represent a ground truth (that is, we neglect representative errors that are however significant) we have used Argo SSS to assess the biases and the standard deviations of the errors of the different SSS products. To that goal, we temporally and spatially collocate the Argo SSS with the SSS maps as follows: every map is compared with the Argo SSS available during the same period (9 days in the case of BEC products) used in the generation of that map. We compare the Argo 280 SSS with the value of the SSS product corresponding to the cell where the Argo is located. Before computing the s match-ups between Arago and SSS product, we apply the following quality control over the values of Argo SSS: -The cut-off depth for Argo profiles is taken between 5 and 10 m.
-Profiles from BioArgo and those included in the greylist (i.e., floats which may have problems with one or more sensors) are discarded.

285
-We use WOA2013 as an indicator: Argo float profiles with anomalies larger than 10 • C in temperature or 5 psu in salinity when compared to WOA2013 are discarded.
-Only profiles having temperature close to surface between -2.5 and 40 o C and salinity between 2 and 41 psu are used.

Singularity analysis
Singularity analysis can be used for the assessment of the geophysical consistency among different products (Umbert et al.,290 https://doi.org/10.5194/essd-2020-232 Open Access Earth System Science Data Discussions Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License. a variable around that point (Turiel et al., 2005(Turiel et al., , 2008a(Turiel et al., , 2009. From an oceanographic point of view, SE are related to the advection term, and therefore, they are intrinsic characteristics of the flow and not specific to the chosen scalar (Nieves et al., 2007). The singularity fronts (bright white streamlines in Figure A11) clearly correspond to the general circulation features, such as the Gulf Stream, the Kuroshio Extension and the tropical instability waves in both the Pacific and Atlantic, among 295 others. This is evident in the singularity fronts derived from OSTIA SST shown in Figure A11. As discussed in Umbert et al. (2014), the SE of different scalars must be in correspondence, and as OSTIA SE have the better quality, we take them as a reference: the SE of all the SSS products will be compared to this reference. Since the OSTIA SST product has the highest resolution (0.05 × 0.05 • ), prior to compute its SE it is first regridded at the same resolution of the SSS product it is going to be compared with. In the case of L3 satellite products and CMEMS, this implies reducing the resolution to 0.25 × 0.25 • and 300 9-day. In the case of BEC L4, since it has the same spatio-temporal grid as OSTIA SST, no regridding is required prior to the calculation of SE.
A way to assess the correspondence of SE is to calculate conditioned histograms of one product SSS SE, h SSS , by the values of OSTIA SE, h SST . Let us recall the definition of the histogram of one random variable y conditioned by one random variable x, denoted by ρ(y|x) and defined by the following expression: where ρ(x, y) is the joint histogram of both variables (the standard 2D histogram) and ρ(x) is the marginal histogram of the variable x (the standard histogram of this variable). In essence, if we put the variable x in the abscissa axis and the variable y in the ordinate axis, the conditioned histogram corresponds to taking the joint histogram of x and y and normalizing it by columns so each column sums up 1.

310
The histogram of a variable conditioned by the value of another variable serves to evidence any functional dependence between both. Conditioned histograms have been used to put in evidence the correspondence of singularity exponents of two variables (Hoareau et al., 2018b), and we will use them for the same purpose here, with h SST the conditioning variable and h SSS the conditioned variable. When the modal line, defined by the maximum probability of h SSS for each value of h SST , is close to a straight line of slope 1 (the identity function), then the singularity exponents of the SSS maps are in good 315 correspondence with those of the reference (Turiel et al., 2008a;Olmedo et al., 2016). In practice, this line becomes horizontal beyond a certain threshold value due to the increase of the error in the estimation of SE for larger values of SE (Turiel et al., 2006), as h SSS and h SST become independent of each other because they are dominated by noise. Notice that when two variables are independent, the conditioned histogram is constant in x and thus the modal line becomes horizontal. In addition to becoming horizontal, the conditioned dispersion will also be larger beyond that threshold, since any dependence between the https://doi.org/10.5194/essd-2020-232 histogram.
-Statistical descriptors: We have computed three descriptors to characterize the conditioned histograms: -The most probable value of h SSS at each bin of h SST : H 0 (h SST ).
-The mean value of h SSS at each bin of h SST :H(h SST ).
-The standard deviation of h SSS at each bin of h SSS : σ H (h SST ).

335
-Correspondence of the SE: As a function of h SST , both H 0 andH should ideally be close to the identitiy, what would be the best correspondence between h SSS and h SST . Although the presence of physical processes other than flow advection, capable of creating singularities and affecting differently both variables (upwelling, rain-induced freshening or river plumes) could lead to a deviation from identity as the SE would not correspond, in general this effect is small and limited to very narrow areas, so the most frequent cause of deviation is just noise, that mainly affects the larger values of andH(h SST ) up to a given threshold value of h SST , h max . We then have compared the corresponding linear regression coefficients, a mode and a mean with 1: the closest a mode and a mean are to 1, the best the geophysical correspondence between SSS and SST in the range of h SST < h max .

345
-Uncertainty analysis: σ H (h SST ) indicates which is the uncertainty associated to the correspondence between the SE of SST and SSS, which is due to effects that we cannot control (including noise, numerical inaccuracies of the algorithm used to compute SE and unknown physical processes Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License.

Spectral analysis
Spectral analysis has been extensively used to analyze satellite observations, in situ data and models outputs, both in the 355 atmosphere and the ocean (Stammer, 1998;Reynolds and Chelton, 2010;Kolodziejczyk et al., 2015;Olmedo et al., 2016;Hoareau et al., 2018b). It is well-known that the Power Density Spectra (PDS) of a scalar submitted to turbulence should follow a power-lay behaviour, characterized by a scaling exponent, sometimes referred to as "spectral slope" as it corresponds to the slope of the log-log representation of the PDS as function of the wavenumber. The analysis of spectral slopes allows obtaining information about the effective spatial resolution of the remote sensing data. For instance, the presence of noise 360 makes the straight curve of log PDS vs log wavenumber to bend and become horizontal at high wavenumbers; this happens because noise is independent of the wavenumber but the amplitude of the signal decreases with the wavenumber, so at a given wavenumber large enough (and thus, at a given spatial scale small enough) noise is dominant: the crossover point signals the effective resolution of the data. Another situation that can appear is when the data is oversmoothened and then there is a systematic lack of energy at high wavenumbers; in those case. a faster-than-linear decay is observed for wavenumber larger 365 than the resolution threshold.
Theoretical studies predicted that temperature and salinity should have the same spectral slopes (Blumen, 1978;Charney, 1971). We have used the spectral analyses based on the PDS and on the Singularity Power Spectra (SPS), that corresponds to the PDS of SE maps (Hoareau et al., 2018b), to estimate the effective spatial resolutions of the SSS products. For this, we compare the spectral slopes of the SSS products with the ones computed from OSTIA SST. PDS slopes are expected to be in 370 between -1 and -3 depending on the dynamical regime that drives the ocean, while the expected SPS slopes for SST and SSS maps range between -2 and -2.5 (Hoareau et al., 2018b). Note that, while the PDS slope is affected or distorted by the presence of noise into the data, the SPS slope is not because the SE algorithm we use (Pont et al., 2013) is designed to filter noise. For the same reason, the use of SPS reduces the uncertainty in the determination of the spectral slopes (Hoareau et al., 2018b).
Therefore, both spectral methods are complementary from a validation point of view as PDS analysis gives information on 375 the effective spatial resolution of the data, and SPS method assesses the existing geophysical structures beyond the remaining sources of uncertainty in remote sensing products.
The spectral analysis approach that we have followed in this work is the following one: -We perform the spectral analysis over the same ocean regions proposed in Hoareau et al. (2018b) (see Figure A4).
-At each on of these regions, we compute the PDS from each SSS product and the OSTIA SST, as well as their cor-380 responding SPS from their SE maps. Both spectra (PDS and SPS) are given as a function of wavenumber values in degree −1 (latitude degrees for meridional regions and longitude degrees for zonal regions) and as wavelength values in kilometers. Recall that a wavelength contains a full oscillation from 0 to +1, then to -1 and finally back to 0, and therefore it contains two resolved points; therefore, to convert wavelengths to resolved spatial scales (resolution scale) the values must be divided by 2. -We have computed and analysed the slope of the averaged PDS and SPS, and we have compare them with the ones resulting from OSTIA SST.

Triple collocation 390
The triple collocation (TC) technique is a powerful tool to estimate the standard deviation of errors of three spatio-temporally collocated measurements of the same target. TC has been used to assess the quality of many remotely sensed variables, and in particular SSS (Hoareau et al., 2018a). The major assumptions of TC are that errors must be uncorrelated with the target variable and also that the errors of the different data sets must be uncorrelated between them. Some refined formulations have been developed in recent years for taking into account the presence of cross-correlated errors between two of the data sets, but 395 they require at least four data sets (Gruber et al., 2016;Pierdicca et al., 2017).
We have used a recently developed formulation of the triple collocation method, the Correlated Triple Collocation (CTC), for the case of three data sets that resolve similar spatial scales from which two of them present correlated errors (González-Gambau, 2020). This TC can be particularly beneficial for the error characterization of variables for which getting measurement systems with uncorrelated errors is challenging or not feasible, and it is particularly well suited to work with limited samples 400 of data because it has a fast convergence with the sampling size. This formulation has been proved for the characterization of radiometric errors in L-band brightness temperatures (TB). By the use of CTC, we have access to maps of errors, so we can characterize which places are less noisy, and also we can ascertain which is the best suited product depending on the location.
The triplets used in this analysis are shown in Table A1, indicating which of the three products is considered with uncorrelated errors with respect to the other two data sets. Errors among the different SMOS products are assumed to be correlated, as 405 well as errors among the different SMAP products since they are measured by the same sensor. Additionally, we have assumed that errors of CMEMS and SMOS CATDS are correlated, since the later is assimilated by the former.
In order to estimate the SSS error of each product, we average the estimated errors resulting from each one of the triplets where the product is considered. We also compute the standard deviation of the estimated error in the different triplets as a metric of the uncertainty in the error estimation.

Comparison with Argo SSS
In the first part of this subsection, we present the comparison of the 9-year time series of the BEC SMOS SSS products with Argo SSS; then, in the second part, we extend the comparison to the rest of SSS products (see section 3.1.1) but for the year 2017 only. In Table A2, the regions of study analyzed in this section are presented.   TR in the same table). The statistics are provided on a yearly basis.
In terms of mean differences with respect to Argo (that would account for biases in the products, if we consider Argo to be the ground truth), the three BEC products (HR, LR, and L4) provide very similar performances. The BEC SSS v2.0 products 420 have a mean difference with respect to Argo below 0.02 and 0.06 psu, in GLO and in TRO regions, respectively. Those values are rather small and may be statistically non-significant.
Regarding the standard deviation of the differences with respect to Argo salinity, among the three salinity fields of the v2 product, BEC HR is the one with the largest standard deviation and L4 the one with the lowest (BEC LR is in between the other two). This is expected since BEC HR is known to have some high spatial frequency noise, while BEC LR is a smoothed 425 version of the BEC HR salinity with a radius of 50 km and therefore a reduction in the noise is expected. The fusion technique used in the generation of the L4 also leads to a reduction of the random noise of the salinity maps, and it is even better than a simple low-pass filtering as it preserves fine-scale structures. The standard deviations of the differences between BEC products and Argo salinity range between 0.34 and 0.26 psu in the case of L3 HR; 0.27 and 0.24 psu in the case of BEC LR; and 0.24 and 0.21 psu in the case of L4. There is a significant reduction of the standard deviation of the differences between SMOS and 430 Argo since 2017, which is more significant in the case of the BEC HR product. This reduction is probably due to a change in the spatial scales of the auxiliary data provided by ESA and what are used in the retrieval. Since 2017, the ECMWF auxiliary data (see section 2.1) are provided at a resolution of 7.8 km, while previously they were provided at 25 km.
We have extended the comparison with Argo to the other SSS products, for the year 2017 only. Figure A5 shows the spatial arrangement of the differences between each salinity product (BEC L4 not shown) and the corresponding collocated Argo SSS, 435 averaged in 5 • × 5 • latitude-longitude cells during the year 2017. In the open ocean, CATDS and REMSS are the two products that provide the lowest differences. CATDS uses In Situ Analysis System (ISAS) SSS (Gaillard et al., 2016;Kolodziejczyk et al., 2018) for the generation of the debiased L3 maps (Boutin et al., 2018). ISAS SSS is an interpolated product of Argo SSS, and this could explain the low differences of CATDS with respect to Argo SSS. JPL displays the largest (positive) differences with respect to Argo SSS. Regarding the BEC products, the three products display similar differences with respect to Argo SSS. We have also analyzed the spatial arrangement of the standard deviations of the differences between the gridded products and Argo SSS in 5 • × 5 • cells. Figure A6 shows these standard deviations for the different products (BEC L4 not shown). CMEMS SSS presents the lowest standard deviation in some regions, such as the Southern Ocean and the Western North Atlantic. BEC  Figure A7 shows the temporal evolution of the mean difference of the gridded products with respect to Argo SSS in 0.25 • bands of latitude (BEC L4 not shown). The product with the lowest differences with respect to Argo SSS is CMEMS, probably because CMEMS assimilate Argo data. Among the satellite products, BEC (HR, LR and L4) and REMSS present the lowest latitudinal differences with respect to Argo. JPL shows the largest differences with respect to Argo SSS (increasing at high latitudes).
We have calculated the statistics of the comparison with Argo SSS in the regions defined in table A2. Table A4 comprises 455 the results of the comparison for all the gridded products and Figure A8 summarizes the statistics by showing the mean (blue square), standard deviation (purple bar), and root mean square difference (RMS, green circle) in all the defined regions. In the plots, the labels correspond to each one of the analyzed products as follows: 1 corresponds to BEC HR; 2 to BEC LR; 3 to CATDS; 4 to JPL; 5 to REMSS; 6 to CMEMS; and 7 to BEC L4. In general, BEC LR and BEC L4 provide very competitive statistics with respect to Argo among all the satellite products. The only products that provide lower RMS differences with 460 respect to Argo in some regions are CMEMS with lower RMS in the Arctic Ocean and North Atlantic Ocean, and REMSS with slightly lower RMS in Tropical and Equatorial regions, Equatorial Pacific, Amazon region and the Indian Ocean. The processing of SMOS data in polar regions and also in semi-enclosed seas deserves specific algorithms, and in fact ESA is funding several projects for developing dedicated products at the Mediterranean Sea, Baltic Sea, Arctic Ocean and Black Sea, in which BEC is involved. 465 We have also analyzed the temporal evolution of the difference with Argo statistics in Figures A9 and A10, for the mean and standard deviation, respectively. The temporal evolution of the mean differences between the three BEC products and Argo are very similar (see blue, green and orange lines in Figure A9) and stable, i.e., no significant oscillations are observed, specially in GLO and TRO regions where the oscillations have an amplitude lower than 0.05 psu. The largest oscillations are observed in the Arctic Ocean, where all the satellite products (except CMEMS) present annual variations larger than 0.4 psu, although not 470 all of these products evolve in the same way. In the case of the BEC products, the differences show a seasonal behaviour, being negative in summer and positive in winter. In summer, fresh water masses coming from ice melting may remain in surface because of stratification. This could produce negative differences between surface salinity (as measured by satellite SSS) and a few-meter deeper salinity (Argo SSS). In winter, high winds may induce significant evaporation and also the formation of sea ice produces excess salt, and this could explain part of the positive differences between measurements at surface and at the 475 first meters. Also notice that in the North Pacific, large differences appear at the end of the year. All the SMOS SSS products present a negative difference with respect to Argo of -0.4 psu. REMSS also presents this negative difference at the beginning but then it suddenly jumps to a positive difference.
Regarding the temporal evolution of the standard deviation, there are some significant seasonal effects. For example, in the northern hemisphere, Arctic and North Atlantic regions, the standard deviation is larger in winter than in spring-summer. This Brazil Current Retroflection that has a seasonal behaviour which is manifested in the SSS (Castellanos et al., 2019), as well as to seasonal changes in the Amazon run off. Figure A12 shows the histograms of the SE of each one of the SSS products conditioned by the SE of OSTIA SST. The modal (H 0 (h SST ), in white) and the mean (H(h SST ), in black) lines are also represented. In Figure A13 the statistical descriptors f independent. This is also expected, as noise becomes dominant as we go to the largest values of the SE, and the noise in SSS and in SST are independent. We can then separate the three ranges in the curve by two tipping points, one in the negative range 500 (that we will denote by h − SST and has a value typically around -0.3 or -0.2) and the other in the positive range (that we will denote by h SST + and has a value around 0.1). The most interesting range is the central one, which is delimited by these two tipping points, where the linear dependence between the SE of SSS and the SE of SST is observed; the larger is this central range (the geophysical correspondence range), the better.

Singularity analysis
All the products present uncorrelation between h SSS and OSTIA h SST from h SST > 0.1. Therefore, we fix h + SST = 0.1 for 505 all SSS products.
We observe that the value of the h − SST depends on the SSS product. Some of the products have good correlation from the most negative SE values while others start the correlated range in moderate negative SE values. For example BEC L4 (top right plot in Figure A12) presents a good correspondence between h SSS and h SST even at the most negative values of h SST .
However, CMEMS product is far from the identity in the most negative values of the h SST , and only at moderated negative 510 values (from h SST ≥ −0.2) we see that the h SST -h SST correspondence becomes closer to the identity.  Table A5 shows the values of a mean and a mode . We observe some differences between a mean and a mode . This is because ofH is a numerically more robust metric (it does not depend on the histogram discretization), but it is more affected by outliers. On the contrary H 0 is numerically less stable (it depends on the histogram 515 discretization) but it is less affected by outliers. In general, despite of their differences, both metrics are consistent when we intercompare them among the different SSS products. The only exception is REMSS that present a much lower a mean than a mode . BEC L4 presents the best performance in all the intervals, having a mean > 0.68 and a mode > 0.88. The CMEMS https://doi.org/10.5194/essd-2020-232  Table A5 shows h + sst − h − sst defined from a mean (11th column) and from a mode (12th column). The products that present the best performances are BEC L4, BEC LR and JPL.

Spectral analysis
Figures A14 and A15 represent the PDS and SPS (respectively) in the different regions for all SSS products and OSTIA SST data. For a better comparison, Figure A16 presents the spectral slopes of the PDS (top) and SPS (bottom) for all the products 535 together. These slopes are calculated in the range of 100-1000 km wavelengths. We observe a large diversity on the shapes of the SSS PDS ( Figure A14) for the different SSS products, significantly differing of the OSTIA PDS (purple line). In contrast, the shapes of SSS SPS are closer to the shape of OSTIA SST SPS up to the 80km wavelength (see Figure A15). This is also observed in Figure A16: the values of PDS slopes vary on a range larger than the one of SPS slopes, which are more concentrated around the theoretically expected range (between -2.0 and -2.5). These results indicate that, despite of the level of 540 noise of each remote sensed product, the geophysical structures of the SSS data are consistent until a 100-80 km wavelength.
However, the slope values in Figure A16 reveal some differences between the products: -BEC HR (in blue) presents the flattest values of the PDS slopes being higher than -1.5 in most of the regions, indicating a strong influence of noise (noise tend to flatten the curve). However, the corresponding SPS slopes remain always in the range of -2 and -2.5. This indicates that even if this product is the one with largest high-frequency noise, BEC HR allows 545 describing consistently the geophysical structures up to 100km wavelength.
-CMEMS product (in grey) presents the steepest PDS slope, becoming lower than -3 in STP and SPURS regions, what means that the product is oversmoothened and that wide regions contain just plainly interpolated data. This indicates that CMEMS have a loss of structures at wavelengths larger than 100km. For example, Figure A14 shows that in STP and SPURS regions, the slope of CMEMS (grey line) becomes steeper at wavelengths around 250km. This is partially confirmed by the SPS slopes that in SPURS and NATL present values lower than -2.5, what indicates that fronts of SE have been lost and confirms the existence of an oversmoothing at wavelengths larger than 100km.
-CATDS (in green) presents the flattest PDS slope in the SPAC region (≈ −1). This suggests that the presence of noise in this region is very large in comparison with the other products and with its performance in other regions. SPAC region is used in the data processing of CATDS to correct systematic and temporal biases (it corresponds to the region where the 555 OTT is computed and applied daily to the CATDS data (Tenerelli and Reul, 2010)). This result suggests that using this region for the calibration of the SMOS measurements could lead to some issues in the resulting product.
-In general, the L3 satellite products of BEC LR and CATDS have similar PDS slopes as well as JPL and REMSS in ITCZ, SPURS and STP regions. This should happen in all the regions, indicating that the spatial resolution of the products depends on the instrument, not on the methods used in the data processing. However, in ARC, NATL and 560 SPAC, the PDS slopes of the BEC LR and CATDS and the ones of JPL and REMSS are significantly different between them, indicating the importance of the data processing methods to achieve proper spatial resolution.
-BEC L4 (in red) presents the closest PDS to those of OSTIA SST. Its SPS slopes remain in the range of -2 and -2.5 in all the regions. This indicates that BEC L4 allows describing spatial scales up to 100km wavelengths with the lowest presence of noise and the closest geophysical consistency with OSTIA SST data.

565
Below 100 km, except for BEC HR, all spectral slopes values get steeper (lack of signal variability into the data). REMSS (cyan line in Figure A14) presents a fast valley-shaped decay around wavelengths of 80km followed by a flattening (traduced by a steeper slope in the SPS slope). This indicates that the smoothing applied onto the REMSS product may remove part of the geophysical variability at those scales. Around 60km wavelength, BEC HR, BEC LR, CATDS and JPL PDS get flattened, while this does not happen in the corresponding SPS spectra. As SPS shape is less affected by noise (Hoareau et al., 2018b), these 570 results indicate that despite of the noise, the geophysical signal present in BEC HR, BEC LR, CATDS, and JPL is consistently described even at those smaller scales, so they can be considered to be valid up to a wavelength of 80 km, which corresponds to a spatial resolution of 40 km. Figure A17 shows the estimated error standard deviations for the different SSS products (CMEMS not shown) for year 2017. We assign one number to each product to assess which is the product with the lowest estimated error standard deviation at each ocean location. Figure A19 shows the four comparisons that have been performed:

Triple collocation
-Comparison of BEC products: We as the label 1 to BEC HR, 2 to BEC LR and 3 to BEC L4. In general BEC L4 is the product with the lowest SSS error. However, BEC LR and BEC HR become more accurate in regions affected by regular 585 rain events (such as the tropics) or continental fresh water discharges (such as Gulf of Mexico), where the hypothesis assumed in the generation of BEC L4 (the gradients of SSS and SST tend to be parallel) does not hold most of the year.
-Comparison of all satellite L3 products: We assign label 1 to BEC HR, 2 to BEC LR, 3 to CATDS, 4 to JPL and 5 to REMSS. In the bulk of the ocean BEC LR provides the lowest SSS error. In some specific regions, such as the Equatorial Atlantic and the Gulf of Guinea, which are regions strongly affected by the dynamics of the Amazon and Congo plumes, 590 the BEC HR provides the best SSS error. Both SMAP products provide better SSS errors in regions affected by RFI (which is expected due to SMAP on-board RFI mitigation) such as the Chinese Sea, close to Madagascar and the Mediterranean Sea. In the Southern Ocean, CATDS is providing the best SSS error.
-Comparison of all satellite products: When we include BEC L4 in the comparison (1-BEC HR, 2-BEC LR, 3-BEC L4, 4-CATDS, 5-JPL and 6-REMSS), the smallest SSS error is given by BEC L4 in the majority of the ocean. As in the 595 previous comparison, BEC HR and BEC LR provide the SSS with the lowest error in regions affected by rainfall and continental discharges. BEC L4 allows improving the SSS estimation in some regions affected by RFI with respect to the L3 products, such as in the China Sea and close to Madagascar, for example. In some regions, SMAP products provide the best SSS and in other BEC L4 is better.
-Comparison among all the products: In this case, BEC L4 remains the product with the lowest salinity error in most of 600 the ocean regions. However, CMEMS SSS also appears to be the product with the lowest salinity error in many regions.
For example, in the ocean regions close to Europe CMEMS SSS provides the best salinity estimation.

Conclusions
We have presented nine years of the new release of SMOS SSS global products generated at the Barcelona Expert Center: the BEC SMOS SSS global L3 and L4 products v2.0. The methods used in their generation include several improvements with 605 respect to the previous version of these products: i) a new latitudinal-seasonal debiasing has been included; ii) an improved filtering criteria based on the salinity geophysical variability have been applied, which allows to better describe the salinity gradients without increasing the overall noise error in the maps; iii) new interpolation schemes are proposed to allow better describing small scale spatial features that are especially relevant in coastal regions; iv) the fusion scheme used in the generation of the L4 product has been modified to preserve small scale spatial features; v) an estimation of the salinity uncertainty is -The statistics of the comparisons with Argo salinity evidence a competitive performance in comparison with the statistics of the rest of SSS products. This includes small mean and standard deviation of the differences with respect to Argo SSS (in the global and regional statistics, latitudinal biases and stable differences in terms of temporal evolution). In this sense, the mean differences with respect to Argo SSS among the three BEC products (BEC L3 (HR and LR) and BEC 620 L4) are very similar (being lower than 0.02 psu at global scale), but the standard deviation is significantly different among them, being the BEC HR the one with the largest standard deviation (being lower than 0.34 psu at global scale) and BEC L4 the one with the lowest one (lower than 0.27 psu).
-In terms of effective spatial resolution and geophysical consistency, we have used two different metrics: -Singularity analysis: the SE of BEC HR and BEC LR SSS products are very similar to the ones of the other satellite 625 salinity products in terms of correlation with OSTIA SST SE. A clear improvement is observed in the BEC L4 that present a higher correlation with the SE of OSTIA SST, suggesting that the geophysical consistency is the most accurate as it is the closest to OSTIA SST SE.
-Spectral analysis: The effective spatial resolutions of BEC HR and BEC LR are consistent with the ones of the other satellite products which are at least a wavelength of 80 km (i.e., spatial resolution of 40 km). The BEC L4 630 presents similar spectral slopes to the ones of OSTIA SST, showing consistent slopes up to 50 km wavelength (25 km spatial resolution). At smaller scales, BEC L4 presents evidence of lack of structure and oversmoothening, so probably it is not resolving scales at its nominal resolution of 0.05 • . In the case of BEC HR, it presents the flattest PDS slopes but with SPS slope values between -2 and -2.5, which indicate that even if the presence of noise is larger, BEC HR is able to represent consistently the geophysical signal.

635
-We have also computed an estimation of salinity errors by using triple collocation. Among the BEC products, BEC L4 provides the SSS field with the lowest error, but in regions strongly affected by rainfall and continental freshwater discharge, the L3 products (BEC HR and BEC LR) are better in terms of salinity error. When we compare all satellite products, BEC L4 remains as the product with the overall minimum salinity error.

Data availability 640
The product is available on the website of the CMEMS Lambda project http://www.cmems-lambda.eu/ and in the Barcelona Expert Center ftp service: http://bec.icm.csic.es/bec-ftp-service/. This product is also available trough EMODNET website             Preprint. Discussion started: 2 October 2020 c Author(s) 2020. CC BY 4.0 License. Table A4. Regional statistics of the differences between the gridded and Argo SSS for the year 2017. For each product and ocean region the mean, standard deviation and the root mean square of the difference are provided separated by :.