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, high-quality 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 assessment. Three SMOS SSS fields are distributed: a high-resolution level-3 product (with DOI https://doi.org/10.20350/digitalCSIC/12601, Olmedo et al., 2020a) consisting of binned SSS in 9 d maps at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$; low-resolution level-3 SSS computed from the binned salinity by applying a smoothing spatial window of 50 km radius; and level-4 SSS (with DOI https://doi.org/10.20350/digitalCSIC/12600, Olmedo et al., 2020b) consisting of daily $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$ maps that are computed by multifractal fusion with sea surface temperature maps. For the validation of BEC SSS products, we have applied a battery of tests aimed 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 the effective spatial resolutions are also 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 the BEC SMOS product with other satellite SSS and reanalysis products. The main outcomes of this quality assessment are as follows. (i) The bias between BEC SMOS and Argo salinity is lower than 0.02 psu at a 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. (iii) The results from triple collocation show the BEC SMOS level-4 product as the product with the lowest estimated salinity error in most of the global ocean and the BEC SMOS high-resolution level-3 as the one with the lowest estimated salinity error in regions strongly affected by rainfall and continental freshwater discharge.

1 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., 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; Kerr et al., 2016) (and references therein).

The Barcelona Expert Center (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 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 recent years, BEC has developed a SMOS SSS internal processing chain that generates SSS maps from SMOS raw data (level 0) to levels 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/, last access: 1 March 2021).

In this work we present the new reprocessing of the BEC SMOS SSS global L3 and L4 products v2.0 for a 9-year period comprising 2011 to 2019, which 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 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 has been introduced.

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 1 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 analyze the effective spatial resolution of each product, using power density spectra (PDS) and singularity power spectra (SPS) (Hoareau et al., 2018b); and (iv) triple collocation analysis to estimate the errors of the different products (González-Gambau et al., 2020).

This paper is structured as follows. Section 2 describes how the BEC SMOS SSS product v2 is generated: Sect. 2.1 introduces the data sets used in the generation of the product and Sect. 2.2 describes the algorithm itself. Section 3 presents the quality assessment of the data: Sect. 3.1 presents the different data sets used for comparison and validation, Sect. 3.2 describes the applied methods and associated metrics, and Sect. 3.3 presents the results of the validation exercise. Section 4 shows where the BEC SMOS SSS global L3 and L4 products are available. Finally, the conclusions are summarized in Sect. 5.

2 Generation of the BEC SMOS SSS global product v2.0

2.1 Data sets used in the generation of the product

2.1.1 SMOS brightness temperature

The brightness temperatures (TBs) obtained from the SMOS MIRAS L1B v620 product provided by ESA are used as the input for the SMOS SSS retrieval. This data set is freely available at https://earth.esa.int/eogateway/catalog/smos-science-products?text=SMOS+Britghness+Temperatures (last access: 1 March 2021).

The L1B v620 product contains the Fourier coefficients of the measured brightness temperatures. Starting from this product, using ESA's Earth Observation Customer Furnished Item (EOCFI) orbit propagation libraries (ESA, 2014) and following a similar procedure to the one used in the operational SMOS level-1 processor chain (Deimos, 2014), the measured TBs are obtained in the antenna reference frame (ARF). The unique difference from the standard processor is the number of points contained per snapshot. The operational processor uses, at antenna level, a hexagonal grid of 128×128 points. The projection of this antenna grid into the ground provides a nominal resolution of about 15 km at boresight. This resolution is more than twice the theoretical SMOS finer resolution (McMullan et al., 2008). Thus, we reduce the computational cost without actually losing information by using an antenna hexagonal grid of 64×64 points for a 30 km resolution at boresight.

2.1.2 Sea surface temperature

The Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) (see Donlon et al., 2012) maps are used as a template in the generation of the Level 4 SSS v2.0 products (see Sect. 2.2.7). OSTIA sea surface temperature (SST) results from the combination of satellite data provided by the Group for High Resolution Sea Surface Temperature (GHRSST) project, combined with in situ observations. The analysis product is obtained using a variant of the optimal interpolation (OI) method described in Martin et al. (2007) at a spatial resolution of 0.05^∘ (approx. 5 km) and daily frequency. OSTIA SST data are provided in netCDF format every day and are freely available at the Copernicus Marine Environment Monitoring Service (CMEMS) service desk at the following link: https://resources.marine.copernicus.eu/?option=com_csw&task=results?option=com_csw&view=details&product_id=SST_GLO_SST_L4_NRT_OBSERVATIONS_010_001 (last access: 1 March 2021).

2.1.3 Auxiliary data used in the salinity retrieval

The auxiliary data used for the SSS retrieval are provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) (Sabater and De Rosnay, 2010) and https://smos-diss.eo.esa.int/oads/access/collection/AUX_Dynamic_Open (last access: 1 March 2021). For each satellite overpass, an ECMWF auxiliary file co-located in time and space with SMOS is provided by ESA. The following fields are used for the retrieval: sea ice cover, sea surface temperature, rain rate, wave model, 10 m wind speed, 10 m neutral equivalent wind (zonal and meridional components), significant height of wind waves, 2 m air temperature, surface pressure, and vertically integrated total water vapor (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 $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ (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 the National Oceanographic Data Center (https://www.nodc.noaa.gov/cgi-bin/OC5/woa13/woa13.pl, last access: 1 March 2021). We have also used the monthly climatology at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ provided by WOA2013 for the correction of latitudinal and seasonal biases (see Sect. 2.2.5).

2.2 Algorithm description

2.2.1 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 (2011–2019) of SMOS SSS maps. This methodology consists of retrieving a single value of SSS from each first Stokes brightness temperature (TB) measurement, which we refer to as raw SSS retrieval. These raw SSS values 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 modeled 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 proposed by Klein and Swift (1977), which relates the TB at flat sea with the SSS and SST, 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_n^{\mathrm{raw}}$ with $n=\mathrm{1},\mathrm{\dots },N$ where N is the total number of raw retrievals in 9 years) 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_n^{\mathrm{raw}}$ values that are acquired under each fixed condition $\mathit{\gamma }=(\mathit{\phi },\mathit{\lambda },d,x,\mathit{\theta })$. Therefore, the systematic errors are the same for all the retrievals in the set $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ with $n=\mathrm{1},\mathrm{\dots }{N}_{\mathit{\gamma }}$ and N_γ the number of retrievals during 9 years with that specific value of γ.

We have defined an estimator of the “typical value” or central estimator of the ensemble $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ that we will call SMOS-based climatology, s^c(γ). We want, by construction, this SMOS-based climatology to represent the sum of a multiyear mean salinity value (which is a geophysical property) and the bias associated with that particular tuple γ (which is of instrumental origin and we want to remove). Therefore, the SMOS-based climatologies can be used for correcting all the retrievals in the set $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$. We have used the same central estimator as the one proposed in Olmedo et al. (2017), that is, the mean around the mode in an interval of ±σ_γ (the standard deviation of $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$).

Then, we have computed SMOS debiased SSS anomalies ($\mathit{\left\{}{s}_n^{\prime }\right(\mathit{\gamma }\left)\mathit{\right\}}$) by subtracting the corresponding SMOS-based climatology s^c(γ) from each individual $s_n^{\mathrm{raw}}\left(\mathit{\gamma }\right)$.

Finally, the debiased salinity value for the acquisition conditions γ, s_n(γ), is computed by adding an external multiyear salinity reference to the $s_n^{\prime }\left(\mathit{\gamma }\right)$; in this case, we have used the annual salinity field provided by WOA2013.

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.

2.2.2 Estimation of SSS error

Each value of raw salinity $s_n^{\mathrm{raw}}$ can be associated with a retrieval error, which is computed according to the following equation:

$$\begin{array}{}\text{(1)}& {\mathit{ϵ}}_n={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\sqrt{({\mathit{\sigma }}_n^{\mathrm{H}}{)}^{\mathrm{2}}+({\mathit{\sigma }}_n^{\mathrm{V}}{)}^{\mathrm{2}}}}{\frac{\partial I_n}{\partial s}}},\end{array}$$

where the ${\mathit{\sigma }}_n^{\mathrm{H}}$ and ${\mathit{\sigma }}_n^{\mathrm{V}}$ are the radiometric resolution for the horizontal (H) and vertical (V) polarizations of the brightness temperature, respectively (which are contained in the ESA L1B product), and $\frac{\partial I_n}{\partial s}$ is the derivative of the modeled first Stokes divided by 2 (I_n) with respect to the salinity (which can be estimated numerically).

2.2.3 Filtering criteria

Filtering out degraded measurements in the generation of the SMOS SSS maps is a key aspect. Without applying any filter, 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 $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ 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 regions.

We apply the following filtering criteria.

• Basic filtering. Any $s_n^{\mathrm{raw}}\left(\mathit{\gamma }\right)$ out of the range of [0, 50] is not considered part of the corresponding set of valid $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$.

• Discarding some full sets of $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$. For a given value of γ, we consider a particular set of $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ valid only 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.

  These filtering criteria are the same as the ones introduced in Olmedo et al. (2017). The only difference is that now the criterion corresponding to the kurtosis is more relaxed. In Olmedo et al. (2017) the set $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ was considered not valid and thus discarded when the kurtosis of the distribution was larger than 4. Now we filter only platykurtic distributions but not leptokurtic ones. Regarding the impact of the filtering criterion corresponding to the skewness, this is the same as the one proposed in Olmedo et al. (2017). This criterion aims at discarding ocean regions affected by RFI contamination. Although some geophysical events tend to be not symmetric and fresh, such as continental discharge and ice melting, and this leads to negative skewed salinity distributions, the typical skewness in these cases is around −0.5. The skewness values lower than −1 typically correspond to distributions that are affected by non-geophysical phenomena. However, we continue revisiting this criterion, and in the next version of the product we will probably analyze the impact of not including this criterion of the skewness.

• Outlier criteria. We discard specific salinity retrievals $s_n^{\mathrm{raw}}\left(\mathit{\gamma }\right)$ when the corresponding SMOS debiased salinity anomaly ($s_n^{\prime }\left(\mathit{\gamma }\right)$) is larger than σ_γ. Since we want to keep the geophysical variability, we include in the criterion a threshold defined by 5σ_φ,λ, with ${\mathit{\sigma }}_{\mathit{\phi },\mathit{\lambda }}^{\mathrm{2}}$ being the geophysical variance of the salinity expected at the grid point (φ,λ). This is new with respect to the criterion proposed in Olmedo et al. (2017). We discard the salinity retrievals that satisfy

  $$\begin{array}{}\text{(2)}& \left|s_n^{\prime }\right(\mathit{\gamma }\left)\right|=\left|s_n^{\mathrm{raw}}\right(\mathit{\gamma })-s^{\mathrm{c}}(\mathit{\gamma }\left)\right|>\sqrt{{\mathit{\sigma }}_{\mathit{\gamma }}^{\mathrm{2}}+\mathrm{25}{\mathit{\sigma }}_{\mathit{\phi },\mathit{\lambda }}^{\mathrm{2}}}.\end{array}$$

  In order to estimate ${\mathit{\sigma }}_{\mathit{\phi },\mathit{\lambda }}^{\mathrm{2}}$ we use SMOS 9 d salinity maps that are computed from a more relaxed filtering criteria, which is

  $$\begin{array}{}\text{(3)}& \left|s_n^{\mathrm{raw}}\right(\mathit{\gamma })-s_{\mathit{\gamma }}^{\mathrm{c}}|<\mathrm{2}{\mathit{\sigma }}_{\mathit{\gamma }}.\end{array}$$

  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). Generally, in the open ocean, σ_φ,λ is small, so Eq. (2) is dominated by σ_γ. Therefore, in the open ocean, the new and the previous filtering criteria have similar performances. However, in those regions with strong salinity dynamics, such as coastal regions, σ_φ,λ is not small, and its contribution in Eq. (2) becomes dominant. Therefore, in those regions with strong salinity dynamics, the new filtering criterion is more relaxed and thus allows better capture of the salinity variability.

• Temporal and geophysical consistency. We temporally and spatially collocate all the debiased retrievals s_n(γ) in 9 d maps with the fixed grid at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. The resulting collocated set of SSS is denoted as $\mathit{\left\{}s\right(t_T,\mathit{\phi },\mathit{\lambda }\left)\mathit{\right\}}$, with t_T being all the acquisition times in the 9 d 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 d period, which means combining all satellite overpass directions (d), across-track distances (x), and incidence angles (θ) that happen at all the times t_T in that period. Then, we consider as valid salinity measurements only those satisfying

  $$\begin{array}{}\text{(4)}& \left|s\right(t_T,\mathit{\phi },\mathit{\lambda })-\overline{s}(T,\mathit{\phi },\mathit{\lambda }\left)\right|<\mathit{\sigma }(T,\mathit{\phi },\mathit{\lambda }),\end{array}$$

  with $\overline{s}(T,\mathit{\phi },\mathit{\lambda })$ and $\mathit{\sigma }(T,\mathit{\phi },\mathit{\lambda })$ being the mean and standard deviation of the set $\mathit{\left\{}s\right(t_T,\mathit{\phi },\mathit{\lambda }\left)\mathit{\right\}}$, that is, all values of SSS at longitude φ, latitude λ, and the 9 d period centered around T. Finally, we average the valid salinity values of $\mathit{\left\{}s\right(t_T,\mathit{\phi },\mathit{\lambda }\left)\mathit{\right\}}$ in that period to obtain the binned 9 d map at $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$, $s_{\mathrm{0}}(T,\mathit{\phi },\mathit{\lambda })$. This criterion was also applied in the previous version of the product. Since, at this step, the salinity retrievals are already debiased and they are temporally and spatially collocated, the criterion of one-sigma applied here is expected to reduce the noise of the level-3 salinity maps only.

2.2.4 Mitigation of temporal biases

SMOS measurements are affected by biases that depend on time (see Martín-Neira et al., 2016). The methodology described in Sect. 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. We use the constant annual reference WOA13 to assess this assumption. The top plot in Fig. 1 shows the temporal evolution of the mean difference between the salinity field provided by WOA13 and the collocated uppermost salinity measurements provided by Argo floats. The results show that this hypothesis is true up to hundreds of practical salinity units. The bottom plot in Fig. 1 shows the difference between the spatially averaged salinity value of $s_{\mathrm{0}}(T,\mathit{\phi },\mathit{\lambda })$ and the spatially averaged salinity value of the annual reference (WOA13). We assume that for a given time T this difference has to be zero. Therefore, we correct each map with this difference, imposing the spatially averaged salinity value of every global map to be equal to the annual reference. We notate the temporally corrected binned salinity fields as $s_{\mathrm{1}}(T,\mathit{\phi },\mathit{\lambda })$.

Figure 1 (a) Assessment of the temporal correction assumption: the graphic represents the mean average of the difference between a constant annual salinity reference (the World Ocean Atlas, 2013) and collocated Argo salinity measurements. (b) Temporal correction applied to the BEC SSS products.

2.2.5 Correction of latitudinal and 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 top panel in Fig. 2). These latitudinal and seasonal biases are known to also happen in other L-band satellite missions and are supposed to be due to the different direct influences 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), which is computed as follows.

• We compute SMOS monthly climatologies ${\overline{s}}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })$ for each month m of the year by averaging all the $s_{\mathrm{1}}(T,\mathit{\phi },\mathit{\lambda })$ values, where T belongs to the same month m of the processed 9 years. Recall that this climatology is defined on a $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ grid, as this is the grid for $s_{\mathrm{1}}(T,\mathit{\phi },\mathit{\lambda })$.

• For each month m, we subtract the WOA2013 monthly climatology, $s^{\mathrm{WOA}}(m,\mathit{\phi },\mathit{\lambda })$, from the corresponding SMOS monthly climatology:

  $$\mathrm{\Delta }{s}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })={\overline{s}}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })-s^{\mathrm{WOA}}(m,\mathit{\phi },\mathit{\lambda }).$$

• We fit $\mathrm{\Delta }{s}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })$ by a second-degree polynomial of the latitude. That is, for every month, m, and every value of latitude in the 0.25^∘ grid, φ, we compute the polynomial p(m,φ)

  $$p(m,\mathit{\phi })=a\left(m\right){\mathit{\phi }}^{\mathrm{2}}+b\left(m\right)\mathit{\phi }+c\left(m\right),$$

  which minimizes the following cost function:

  $$\sum _{\mathit{\lambda }}{\left(\mathrm{\Delta }{s}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })-p(m,\mathit{\phi })\right)}^{\mathrm{2}}.$$

• After computing the optimal polynomials p(m,φ), we correct the maps $s_{\mathrm{1}}(t,\mathit{\phi },\mathit{\lambda })$ by interpolating the polynomial p(m,φ) daily to the specific moment of the month. We denote the latitudinal–seasonal debiased SSS by $s_{\mathrm{2}}(T,\mathit{\phi },\mathit{\lambda })$.

In the bottom plots of Fig. 2 the monthly interpolating polynomials p(m,φ) are presented (in blue), as well as the mean difference $\mathrm{\Delta }{s}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })$ (in green). As observed in the plots, the approach of this correction has some limitations at high latitudes, where the sea ice dynamics also induce ice–sea contamination.

Figure 2 (a) Latitudinal and seasonal bias affecting the BEC SMOS SSS L3 maps after applying the systematic bias correction proposed in Sect. 2.2.1 and the temporal correction proposed in Sect. 2.2.4. The plot represents a Hovmöller diagram of the temporal evolution of the differences between the BEC SMOS SSS and the Argo SSS in latitudinal bins of 0.25^∘. (b) In the last four rows the monthly interpolating polynomial (in blue) is presented as well as the mean difference between the monthly SMOS climatology and the monthly WOA13 (in green) (i.e., $\mathrm{\Delta }{s}_{\mathrm{1}}(m,\mathit{\phi },\mathit{\lambda })$; see text in Sect. 2.2.1).

2.2.6 Mitigation of residual spatial biases

After applying all the above corrections, we make the last check. By construction, at each geographical location the average salinity of the full period should be equal to the multiyear reference introduced in Sect. 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 $\mathit{\left\{}{s}_n^{\mathrm{raw}}\right(\mathit{\gamma }\left)\mathit{\right\}}$ for a given γ significantly deviate 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 (Sect. 2.2.3) are significant. In order to mitigate this last bias, we remove the map corresponding to the difference between the mean average of all the SMOS SSS in the 2011–2019 period ($s_{\mathrm{2}}(T,\mathit{\phi },\mathit{\lambda })$) and WOA2013 (see Fig. 3). The resulting salinity field is our L3 high-resolution product $s^{L\mathrm{3}}(T,\mathit{\phi },\mathit{\lambda })$. In future versions of this product we will introduce a better definition of the SMOS-based climatology to avoid this last correction step.

Figure 3 Difference between the 2011–2019 average of BEC SSS maps after applying corrections described in Sect. 2.2.1, 2.2.4, and 2.2.5 and WOA2013 SSS.

2.2.7 Multifractal fusion techniques

We use the multifractal fusion techniques introduced in Umbert et al. (2014) to increase the spatial and temporal resolutions 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 (SEs). 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, 2008b), it has been observed that other scalars, such as chlorophyll concentration 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 SEs. 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:

$$\begin{array}{}\text{(5)}& \mathrm{SSS}=a\cdot \mathrm{SST}+b,\end{array}$$

where a and b are smooth functions; that is they must have small gradients, as otherwise they would introduce additional SEs. 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-weighting function as in Olmedo et al. (2016), that is, the inverse of the fourth 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=\mathrm{2.5}{}^{\circ }$ 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 $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$, are obtained.

2.2.8 BEC SMOS SSS global product v2.0

The BEC SMOS SSS L3 global product v2.0 consists of 9 d SSS maps at a regular grid of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$ generated daily. 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

• the BEC SMOS HR SSS product, where HR stands for high resolution and contains the binned salinity field $s^{L\mathrm{3}}(T,\mathit{\phi },\mathit{\lambda })$, and

• the BEC SMOS LR SSS product, where LR stands for low resolution and is a low-pass-filtered version of $s^{L\mathrm{3}}(T,\mathit{\phi },\mathit{\lambda })$, computed by applying a smoothing windows of radius of 50 km. This product is denoted by $s_{\mathrm{low}}^{L\mathrm{3}}(T,\mathit{\phi },\mathit{\lambda })$.

The BEC SMOS L4 SSS product v2.0 (hereafter BEC L4) consists of daily SSS maps at a regular grid of $\mathrm{0.05}{}^{\circ }\times \mathrm{0.05}{}^{\circ }$. This product is denoted by $s^{L\mathrm{4}}(T,\mathit{\phi },\mathit{\lambda })$.

3 Quality assessment

3.1 Data sets for validation

3.1.1 Satellite sea surface salinity

We have compared the performance of the new BEC products with that of other satellite SSS products. We have centered the validation in the year 2017 because there are not any large-scale geophysical phenomena (such as El Niño or La Niña events) and also because SSS products produced by the National Aeronautics and Space Administration (NASA) Soil Moisture Active Passive (SMAP) mission are available (this mission has been operating since early 2015; Entekhabi et al., 2010). We have also computed the statistics with Argo floats for all the products for the year 2016 to assess the robustness of the results in 1 year in which El Niño happened. The satellite SSS products used for the intercomparison are as follows.

• CATDS SMOS products. The 9 d SMOS SSS maps provided by Centre Aval de Traitement des Données SMOS (CATDS). We use the L3 debiased v4 freely available at http://catds.ifremer.fr/Products/Available-products-from-CEC-OS/Locean-v2019 (last access: 1 March 2021). This product decreases the mean bias over the open ocean with respect to previous versions, and it improves ice filtering, which leads to an improvement of SSS at high latitudes, especially in the Southern Ocean (Boutin et al., 2016, 2018; Kolodziejczyk et al., 2016). See also https://www.catds.fr/Products/Available-products-from-CEC-OS/CEC-Locean-L3-Debiased-v4 (last access: 1 March 2021) for the description of the last improvements of the product.

• JPL SMAP products. The 8 d SMAP SSS maps are provided by Jet Propulsion Laboratory (JPL). We use the level-3 version 4.2 freely available at https://podaac-opendap.jpl.nasa.gov/opendap/allData/smap/L3/JPL/V4.2/ (last access: 1 March 2021). Updates of version 4.2 with respect to previous versions include improvement in the TB calibration using an adjusted reflector emissivity, the inclusion of a SST dependence on a flat surface emissivity model, use of updated land correction tables, and inclusion of averaged ancillary ice concentration data (Fore et al., 2016).

• REMSS SMAP products. The 8 d running Remote Sensing Systems SMAP Level 3 Sea Surface Salinity Standard Mapped Image version v4 is used, which is freely available at http://www.remss.com/missions/smap (last access: 1 March 2021). In particular, we have used the field sss_smap, which is a smoothed measurement at approximately 70 km resolution. The major change in Version 4.0 from Version 3.0 is an improved land correction, which allows for SMAP salinity retrievals closer to the coast (Meissner et al., 2018).

3.1.2 In situ salinity: Argo floats

For the purpose of direct comparison of values, we have used in situ salinity data obtained by Argo profilers Argo (2000). We consider the uppermost Argo salinity between 5 and 10 m depth (hereafter Argo SSS). Argo data are collected and made freely available by the International Argo Program and the national programs that contribute to it (http://www.argo.ucsd.edu, last access: 1 March 2021, http://argo.jcommops.org, last access: 1 March 2021). The Argo Program is part of the Global Ocean Observing System.

3.1.3 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 near-real-time ARMOR3D (Nardelli, 2012; Nardelli et al., 2016; Droghei et al., 2016) corresponding to the year 2017. We use version 4 of ARMOR3D, which is freely available at 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, last access: 1 March 2021). 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 measurement obtained from ISAS-CORA (both distributed through CMEMS) and a high-pass filter of Reynolds SST L4 satellite observations. This product assimilates SMOS SSS generated and distributed by CATDS.

3.1.4 Sea surface temperature

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

3.2 Validation methods

3.2.1 Comparison with Argo

Assuming that Argo values represent a ground truth (which 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 d in the case of BEC products) used in the generation of that map. We compare the Argo SSS with the value of the SSS product corresponding to the cell where the Argo is located. Before computing the match-ups between Argo and SSS products, 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 grey list (i.e., floats which may have problems with one or more sensors) are discarded.

• 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 between −2.5 and 40 ^∘C and salinity between 2 and 41 psu close to the surface are used.

3.2.2 Singularity analysis

Singularity analysis can be used for the assessment of the geophysical consistency among different products (Umbert et al., 2014; Hoareau et al., 2018b). From an oceanographic point of view, SEs 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 Fig. 4) 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 others. This is evident in the singularity fronts derived from OSTIA SST shown in Fig. 4. As discussed in Umbert et al. (2014), the SEs of different scalars must correspond, and as OSTIA SEs have the better quality, we take them as a reference: the SEs of all the SSS products will be compared to this reference. Since the OSTIA SST product has the highest resolution ($\mathrm{0.05}\times \mathrm{0.05}{}^{\circ }$), prior to computing 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 $\mathrm{0.25}\times \mathrm{0.25}{}^{\circ }$ and 9 d. 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 SEs.

Figure 4 SE for the 15 August 2017. (a) SE of BEC LR SSS product. (b) SE of OSTIA SST.

A way to assess the correspondence of SEs 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:

$$\begin{array}{}\text{(6)}& \mathit{\rho }\left(y\right|x)\equiv {\displaystyle \frac{\mathit{\rho }(x,y)}{\mathit{\rho }\left(x\right)}},\end{array}$$

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 to 1.

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 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 in 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 variables would decrease the dispersion. Thus, besides analyzing the range of the singularity correspondence, the conditioned dispersion (i.e the standard deviation of the h_SSS for each bin of h_SST) must also be considered.

In the following, we detail the methodology employed in this work.

• Calculation of the histograms. Regarding the calculation of the histograms shown in this paper, for a given SSS product we have accumulated all the values of the pairs (h_SST, h_SSS) of all the globe and for the whole year 2017 to create the joint histogram. Notice that, as shown in Umbert et al. (2014), the histograms of SEs do not change in space or in time, and thus this accumulation can be done. The SE values used to construct the histograms have been limited to the range at [−0.4, 0.4] (which typically contains more than 80 % of all possible values), and we fix the bin size to 0.05 for both h_SST and h_SSS. The marginal histogram ρ(h_SST) is obtained by summing the values at each column, and the conditioned histogram is obtained by a simple division of all the elements of a column by the associated value of the marginal 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₀(h_SST);

  • the mean value of h_SSS at each bin of h_SST: $\overline{H}\left(h_{\mathrm{SST}}\right)$;

  • the standard deviation of h_SSS at each bin of h_SSS: σ_H(h_SST).

• Correspondence of the SE. As a function of h_SST, both H₀ and $\overline{H}$ should ideally be close to the identity, 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 both variables differently (upwelling, rain-induced freshening, or river plumes) could lead to a deviation from identity as the SEs 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, which mainly affects the larger values of SEs and for that reason the deviations show up as the value of SE increases. To quantify the quality of the correspondence between the SEs of SST and those of SSS, we have computed the linear regression coefficients of the functions H₀(h_SST) and $\overline{H}\left(h_{\mathrm{SST}}\right)$ 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 closer a_mode and a_mean are to 1, the better the geophysical correspondence between SSS and SST in the range of h_SST<h_max.

• Uncertainty analysis. σ_H(h_SST) indicates which is the uncertainty associated with 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 SEs, and unknown physical processes). We have computed the mean value of σ_H(h_SST) in the same range of values of h_SST used to compute a_mean and a_mode, namely h_SST<h_max. The goal is to have a value of mean σ(h_SSST) as small as possible compared to the marginal standard deviation of h_SST.

• Validity range. Finally, we have also analyzed which is the range of good correspondence between h_SST and h_SSS. The larger the range, the better, as this would ensure that more geophysical structures (fronts, eddies, filaments) are described by both variables.

3.2.3 Spectral analysis

Spectral analysis has been extensively used to analyze satellite observations, in situ data, and model outputs, in both the 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-law behavior, 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 us to obtain information about the effective spatial resolution of the remote sensing data. For instance, the presence of noise makes the straight curve of log PDS vs. log wavenumber 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 are oversmoothed, and then there is a systematic lack of energy at high wavenumbers; in those cases, a faster-than-linear decay is observed for wavenumbers larger 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), which correspond 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 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 in 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 the effective spatial resolution of the data, and the 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 Fig. 5).

• At each on of these regions, we compute the PDS from each SSS product and the OSTIA SST, as well as their corresponding SPS from their SE maps. Both spectra (PDS and SPS) are given as a function of wavenumber values per degree (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. Thus, to convert wavelengths to resolved spatial scales (resolution scale), the values must be divided by 2.

• For each region and product, we have computed the mean PDS and the mean SPS over the full year of 2017 to reduce the fluctuations of each individual spectrum.

• We have computed and analyzed the slope of the averaged PDS and SPS, and we have compared them with the ones resulting from OSTIA SST.

Figure 5 Spatial distribution of Rossby radii of deformation, from Hoareau et al. (2018b). Blue boxes are the regions where the power density spectra and singularity power spectra are computed: NATL (North Atlantic) [0–64^∘ N, 27–33^∘ W], Intertropical Convergence Zone (ITCZ) [4–10^∘ N, 96–160^∘ W], SPURS (region of the oceanographic campaign Salinity Processes in the Upper-ocean Regional Study) [22–28^∘ N, 28–60^∘ W], ARC (Agulhas Current Retroflection) [39–45^∘ S, 30–94^∘ E], STP (south tropical Pacific) [27–33^∘S, 106–170^∘ W], and SPAC (southern Pacific) [14^∘ N–50^∘ S, 122–128^∘ W].

3.2.4 Triple collocation

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 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 et al., 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 of data because it has a fast convergence with the sampling size. This formulation has been proven for the characterization of radiometric errors in L-band brightness temperatures (TBs). By the use of CTC, we have access to maps of errors, so we can characterize which places are less noisy, and we can also ascertain which is the best suited product depending on the location.

The triplets used in this analysis are shown in Table 1, 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 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 latter is assimilated by the former.

Table 1 SSS triplets used in the analysis performed in Sect. 3.3.4.

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.

3.3 Validation results

3.3.1 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 the SSS products (see Sect. 3.1.1) but for the year 2017 only. In Table 2, the regions of study analyzed in this section are presented.

Table 2 Ocean regions used in the statistics with respect to Argo SSS.

Table 3 Statistics of the comparison of the BEC global SSS product v2.0 with Argo.

Table 3 summarizes the statistics of the differences of BEC SMOS SSS v2.0 products minus Argo SSS for the 2011–2019 period in two different ocean regions: (i) the whole temperate ocean (labeled as GLO in Table 2) and (ii) the tropics (labeled TR in the same table). The statistics are provided on a yearly basis.

In terms of mean differences with respect to Argo (which 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 have a mean difference with respect to Argo below 0.02 and 0.06 psu, in GLO and in TRO, 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 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 Argo since 2017, which is more significant in the case of the BEC HR product. We have analyzed several possible reasons for this reduction. One reason could be a reduction of RFI contamination, which has been observed at the global scale since 2017. Another reason could be 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 Sect. 2.1) are provided at a resolution of 7.8 km, while previously they were provided at 25 km.

We have calculated the statistics of the comparison with Argo SSS in the regions defined in Table 2. Tables 4 and 5 comprise the results of the comparison for all the gridded products for the year 2017, and Fig. 9 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 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 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, the equatorial Pacific, the 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 there are several projects for developing dedicated products at the Mediterranean Sea, Baltic Sea, Arctic Ocean, and Black Sea, in which BEC is involved. We have also computed the statistics of the comparison with Argo and all the satellite products for the year 2016 (see Table 6 to assess the performance in 1 year in which an El Niño event occurred). The results are similar to the ones obtained in 2017. The BEC products present a slight increase in the standard deviation of the differences with respect to Argo in all the regions that is consistent with the general increase observed at the first 6 years of the mission and that have been mentioned above.

Table 4 Regional statistics of the differences between the SMOS and Argo SSS for the year 2017. For each product and ocean region the mean, standard deviation, root mean square of the difference, and coefficient of correlation r² are provided separated by a colon.

Table 5 Regional statistics of the differences between the SMAP and Argo SSS for the year 2017. For each product and ocean region the mean, standard deviation, root mean square of the difference, and coefficient of correlation r² are provided separated by a colon. The table also contains the results for the CMEMS product.

Table 6 Regional statistics of the differences between the satellite and Argo SSS for the year 2016. For each product and ocean region the mean, standard deviation, and root mean square of the difference are provided separated by a colon.

Figure 6 shows the spatial arrangement of the differences between each salinity product (BEC L4 not shown) and the corresponding collocated Argo SSS, averaged in $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ latitude–longitude cells during the year 2017. In the open ocean, CATDS and REMSS are the two products that provide the lowest differences. 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. Significant positive differences (≈0.2 psu) are evidenced in the North Atlantic Ocean and also in the North Pacific. Similar differences are found in the CMEMS product (especially in the northern Atlantic and the eastern northern Pacific).

Figure 6 Spatial distribution of the mean differences with respect to Argo SSS. From left to right and top to bottom: BEC HR, BEC LR, JPL, REMSS, CATDS, and CMEMS.

Figure 7 Spatial distribution of the standard deviations of the differences with respect to Argo SSS. From left to right and top to bottom: BEC HR, BEC LR, JPL, REMSS, CATDS, and CMEMS.

We have also analyzed the spatial arrangement of the standard deviations of the differences between the gridded products and Argo SSS in $\mathrm{5}{}^{\circ }\times \mathrm{5}{}^{\circ }$ cells. Figure 7 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 HR is the product with the largest standard deviation. The rest of the satellite products show similar standard deviations. The largest standard deviations are located in regions of strong salinity gradients such as the Gulf Stream and close to the mouth of the main rivers. Discrepancies appear between the satellite products in the Southern Ocean, where SMOS products show a lower standard deviation of the difference than SMAP products.

Figure 8 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 assimilates 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).

Figure 8 Hovmöller diagrams of the mean difference between salinity gridded products and the temporally and spatially collocated uppermost salinity measurement provided by Argo floats. The y axis represents the latitudes: differences are averaged in latitude bins of 0.25^∘. The x axis represents the time (in days): only the year 2017 is considered in this analysis. The gridded products used in this analysis are the following ones: (from left to right and top to bottom) BEC HR, BEC LR, JPL, REMSS, CATDS, and CMEMS.

We have also analyzed the temporal evolution of the difference with Argo statistics in Figs. 10 and 11, 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 Fig. 10) 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 all of these products evolve in the same way. In the case of the BEC products, the differences show a seasonal behavior, being negative in summer and positive in winter. In summer, fresh water masses coming from ice melting may remain at the surface because of stratification. This could produce negative differences between surface salinity (as measured by satellite SSS) and salinity a few meters deeper (Argo SSS). In winter, high winds may induce significant evaporation, and the formation of sea ice also produces excess salt, and this could explain part of the positive differences between measurements at the surface and at the 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.

Figure 9 Regional statistics of the comparison with respect to Argo salinity. The numbers in the label plots correspond to the different products that are compared with Argo SSS: 1-BEC HR, 2-BEC LR, 3-CATDS, 4-JPL, 5-REMSS, 6-CMEMS, and 7-BEC L4.

Figure 10 Mean of the differences between the gridded SSS products and Argo salinity in different ocean regions for the year 2017.

Figure 11 Standard deviation of the differences between the gridded SSS products and Argo salinity in different ocean regions for the year 2017.

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 is expected because L-band TBs are less sensitive to the SSS in cold waters (wintertime) than in warm waters (summertime), which implies that the retrievals of SSS must be noisier in winter than in summer. However in the North Pacific, all the satellite products present the inverse behavior: a reduction of the standard deviation is present at the end of the year. The reason for this decrease is still under study. In the Amazon River region, the standard deviation increases in spring and summer. The reason for this increase could be related to the seasonal behavior of the North Brazil Current and the North Brazil Current retroflection that has a seasonal behavior which is manifested in the SSS (Castellanos et al., 2019) as well as to seasonal changes in the Amazon run off.

3.3.2 Singularity analysis

In this section we take the SE of OSTIA SST as a reference to assess the effective spatial resolution of the salinity products. The OSTIA SST product is not perfect and it also has some limitations in describing the small spatial gradients of SST, which could be reflected in the results of this comparison.

Figure 12 shows the histograms of the SEs of each one of the SSS products conditioned by the SEs of OSTIA SST. The modal (H₀(h_SST), in black) and the mean ($\overline{H}\left(h_{\mathrm{SST}}\right)$, in white) lines are also represented. In Fig. 13 the statistical descriptors (conditioned mean $\overline{H}\left(h_{\mathrm{SST}}\right)$, conditioned mode H₀(h_SST), and conditioned standard deviation σ_H(h_SST)) of all the products are shown.

Figure 12 Histogram of SEs of SSS conditioned to the SEs of OSTIA SST. For each SST SE bin, the corresponding SSS SE distribution is normalized by the total number of SSS SEs. The black line corresponds to the mode of the SE SSS at each bin of SE SST (H₀ in the text). The white line corresponds to the mean SE SSS at each bin of SE SST ($\overline{H}$ in the text). In the first row from left to right, the SSS products correspond to BEC HR, BEC LR, and BEC L4. In the second row from left to right the products correspond to CATDS, JPL, and REMSS, and in the third row CMEMS products are shown.

All the conditioned histograms present three different well-defined ranges. The first part of the curve is unstructured and noisy: this is normal because there are very few points with those values, so the statistics are scarce and fluctuations are large; it is also affected by small mismatches in the positions of fronts, lack of accuracy of the sensors, and, very occasionally, different singularity-inducing effects acting on each variable. Then, we find a central part where the relation between h_SST and h_SSS is clearly linear, as is evident from the H₀ and $\overline{H}$ curves. After this, we find a third range, where the value of the conditioned histogram is saturated and the values of both H₀ and $\overline{H}$ are horizontal lines, indicating that h_SST and h_SSS are 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 is independent. We can then separate the three ranges in the curve by two tipping points, one in the negative range (which we will denote by $h_{\mathrm{SST}}^{-}$ and has a value typically around −0.3 or −0.2) and the other in the positive range (which we will denote by $h_{\mathrm{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 SEs of SSS and the SEs of SST is observed; the larger this central range (the geophysical correspondence range), the better.

All the products present no correlation between h_SSS and OSTIA h_SST from h_SST>0.1. Therefore, we fix $h_{\mathrm{SST}}^{+}=\mathrm{0.1}$ for all SSS products.

We observe that the value of the $h_{\mathrm{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. BEC L4 (top right plot in Fig. 12) presents a good correspondence between h_SSS and h_SST even at the most negative values of h_SST. This is expected because the BEC L4 product is computed from OSTIA SST by applying multifractal fusion. The method improves the quality of the salinity maps by using the spatial structures of the OSTIA SST. Therefore, the singularity exponents of the BEC L4 are expected to be close to the ones of the OSTIA SST. Therefore, in this section, BEC L4 is expected to provide the best performance. The CMEMS product is far from the identity in the most negative values of the h_SST, and only at moderated negative values (from $h_{\mathrm{SST}}\ge -\mathrm{0.2}$) do we see that the h_SST–h_SSS correspondence becomes closer to the identity.

Quality of the geophysical correspondence. We have computed a_mean and a_mode (and also 〈σ_H〉) in three intervals of h_SST: $[-\mathrm{0.4},\mathrm{0.1}]$, $[-\mathrm{0.3},\mathrm{0.1}]$, and $[-\mathrm{0.2},\mathrm{0.1}]$. Table 7 shows the values of a_mean and a_mode. We observe some differences between a_mean and a_mode. This is because $\overline{H}$ 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₀ is numerically less stable (it depends on the histogram discretization) but it is less affected by outliers. In general, despite their differences, both metrics are consistent when we intercompare them among the different SSS products. The only exception is REMSS, which presents 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 product also presents very high values of a_mean and a_mode but only in the interval of $[-\mathrm{0.2},\mathrm{0.1}]$; as mentioned above, in the CMEMS product, the relation between h_SSS and h_SST in the most negative range is far from the identity. All the satellite SSS products provide negative values of $\overline{H}$ and H₀ for $h_{\mathrm{SST}}\in [-\mathrm{0.35},-\mathrm{0.25}]$, while the CMEMS product provides positive values. This suggests that there are some SSS fronts of moderate intensity that are captured by all the SSS satellite products but not captured at all by the CMEMS product. In the case of L3 satellite products, a_mean and a_mode range in between 0.29 and 0.5 depending on the product and the interval of analysis.

Table 7 Singularity analysis metrics of the SSS products over three different SST SE regimes: $[-\mathrm{0.4},\mathrm{0.1}]$, $[-\mathrm{0.3},\mathrm{0.1}]$, and $[-\mathrm{0.2},\mathrm{0.1}]$. For each one of these regimes we compute the linear regression coefficient of the mean SSS SE as a function of SST SE (mean, $\overline{H}\left(h_{\mathrm{SST}}\right)$); the linear regression coefficient of the most probable value of SSS SE as a function of the SST SE (mode, H₀(h_SST)) and the averaged value of the standard deviation of the SSS SE as a function of SST SE in the corresponding range (〈σ_H〉).

Residual uncertainty. For each of the ranges above, we have calculated the average on that range of values of h_SST of the conditioned standard deviation, 〈σ_h〉. The best performance is for BEC L4, presenting a 〈σ_h〉 of 0.11 in all the analyzed ranges. The worse performances are for the CMEMS product that reaches values of 0.17 in the range of $[-\mathrm{0.4},\mathrm{0.1}]$. The L3 satellite products provide similar 〈σ_H〉 values that vary from 0.13 to 0.15 (depending on the product and the analyzed range).

Extent of the geophysical correspondence range. In order to determine the largest possible range of reasonable geophysical correspondence between SEs of SSS and of SST, we have defined $h_{\mathrm{SST}}^{-}$ as the most negative value of h_SST with a_mean (or a_mode) larger than 0.35. Table 7 shows $h_{\mathrm{sst}}^{+}-h_{\mathrm{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.

3.3.3 Spectral analysis

Figures 14 and 15 represent the PDS and SPS (respectively) in the different regions for all SSS products and OSTIA SST data. We have also included the information about the standard deviation of the set of individual spectra that have been used in the mean spectra represented in these figures. The vertical line represented in the plots applies to all the frequencies on the log-log scale in which the spectra are represented.

Figure 13 From left to right in the following are shown: the most probable SSS SE value as a function of the SST SEs, the mean SSS SE value as a function of the SST SEs, and the standard deviation of the SSS SEs for each SST SE.

For a better comparison, Fig. 16 presents the spectral slopes of the PDS (top) and SPS (bottom) for all the products together. These slopes are calculated in the range of 100–1000 km wavelengths. Although the effect of rain and other geophysical phenomena can lead to differences in the PDS for SSS and SST, on the spatial and temporal scales in which we have computed the slopes, these differences are expected to be negligible.

Figure 14 Power density spectra of the different SSS products.

Figure 15 Singularity spectral analysis of the different SSS products.

Figure 16 Slopes of the power density spectra (a) and the singularity specter analysis (b) of the different SSS products. The corresponding slopes of the PDS and SPS of OSTIA SST are also included as a reference.

We observe a large diversity in the shapes of the SSS PDS (Fig. 14) for the different SSS products, significantly differing for the OSTIA PDS (purple line). In contrast, the shapes of SSS SPS are closer to the shape of OSTIA SST SPS up to the 80 km wavelength (see Fig. 15). This is also observed in Fig. 16: 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 the level of noise of each remotely sensed product, the geophysical structures of the SSS data are consistent until a 100–80 km wavelength. However, the slope values in Fig. 16 reveal some differences between the products.

• BEC HR (in orange) 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 tends to flatten the curve). However, the corresponding SPS slopes always remain in the range of −2 and −2.5. This indicates that even if this product is the one with the largest high-frequency noise, BEC HR allows a consistent description of the geophysical structures up to the 100 km wavelength.

• The CMEMS product (in grey) presents the steepest PDS slope, becoming lower than −3 in the STP and SPURS regions, which means that the product is oversmoothed and that wide regions contain just plainly interpolated data. This indicates that the CMEMS product has a loss of structures at wavelengths larger than 100 km. For example, Fig. 14 shows that in the STP and SPURS regions, the slope of the CMEMS product (grey line) becomes steeper at wavelengths around 250 km. This is partially confirmed by the SPS slopes that in SPURS and NATL present values lower than −2.5, which indicates that fronts of SEs have been lost and confirms the existence of an oversmoothing at wavelengths larger than 100 km.

• 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. The SPAC region is used in the data processing of CATDS to correct systematic and temporal biases (it corresponds to the region where the ocean target transformation (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.

• 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 the description of spatial scales up to 100 km wavelengths with the lowest presence of noise and the closest geophysical consistency with OSTIA SST data. That is partially because of the fusion method used in the generation of the BEC L4, which improves the salinity by introducing the spatial correlation consistency with respect to OSTIA SST.

Below 100 km, except for BEC HR, all spectral slope values get steeper (lack of signal variability into the data). REMSS (cyan line in Fig. 14) presents a fast valley-shaped decay around wavelengths of 80 km followed by a flattening (traduced by a steeper slope in the SPS slope). This indicates that the smoothing applied to the REMSS product may remove part of the geophysical variability at those scales. Around the 60 km 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 results indicate that despite 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 17 SSS error estimation by triple collocation, from top to bottom and left to right: BEC HR, BEC LR, BEC L4, JPL, REMSS, and CATDS.

Figure 18 Uncertainty in the SSS error estimation by triple collocation from top to bottom and left to right: BEC HR, BEC LR, BEC L4, JPL, REMSS, and CATDS.

Figure 19 Spatial distribution of the products with the minimum SSS estimated error. (a) The product with the lowest error among the BEC products (1-BEC HR, 2-BEC LR, and 3-BEC L4). (b) The product with the lowest error among the L3 satellite products (1-BEC HR, 2-BEC LR, 3-CATDS, 4-JPL, and 5-REMSS). (c) The product with the lowest error among the satellite products (1-BEC HR, 2-BEC LR, 3-BEC L4, 4-CATDS, 5-JPL, and 6-REMSS). (d) The product with the lowest error among all the products analyzed in the study (including reanalysis) (1-BEC HR, 2-BEC LR, 3-BEC L4, 4-CATDS, 5-JPL, and 6-REMSS; 7-CMEMS).

3.3.4 Triple collocation

Figure 17 shows the estimated error standard deviations for the different SSS products (CMEMS not shown) for the year 2017. As a general remark, the estimated error standard deviations are larger in those regions with higher salinity dynamics, such as the Gulf of Bengal and the equatorial Atlantic which is affected by the dynamics of the Amazon River plume. In general, BEC L4 provides the lowest error standard deviations among all the satellite products, with the exception of a region close to the Antarctic ice edge, where CATDS provides the lowest one. The uncertainty on the estimation of the error standard deviation (calculated as the standard deviation across all the possible triples of the error standard deviations) is lower than 0.03 psu for all the products and almost all ocean regions as shown in Fig. 18.

We assign one number to each product to assess which is the product with the lowest estimated error standard deviation at each ocean location. We have not included the uncertainty associated with the estimation of the salinity error in this computation. This implies that, although we represent one map with a single product only at each grid point, after considering the uncertainty of each estimation, several products could provide a similar performance from a statistical point of view. In particular, the estimated error in regions of large salinity variability presents larger uncertainty than in regions of low salinity variability (see Fig. 18). Therefore, the following maps are less accurate in these regions. Figure 19 shows the four comparisons that have been performed.

• Comparison of BEC products. We assign 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 rain events (such as the tropics) or continental fresh water discharges (such as the 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 the 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, 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 provides 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 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 improvement of 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. In some regions, SMAP products provide the best SSS and in others 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 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.

The results obtained from triple collocation provide a complementary view to the comparisons with Argo floats (see Sect. 3.3.1). Although the comparison between in situ and satellite data provides very valuable information about the quality of the satellite products, this comparison has several limitations that the triple collocation method does not.

• For sampling, the in situ measurements are provided over a few samples while the satellite data are synoptic. The dynamics displayed by the in situ measurement could be strongly conditioned by its sampling. Therefore, the results from the comparison could not be completely representative of the quality of the satellite product in the considered region.

• The spatial and temporal scales of the in situ and satellite measurements are different. The in situ measurements provide punctual and instantaneous measurements while satellite measurements correspond to an integrated measure of several days and a footprint of several square kilometers.

• The in situ measurements are typically given at several meters depth while the satellite data provide the salinity at a few centimeters depth.

The above-mentioned reasons could lead to apparent contradictory results between the results corresponding with Argo comparison and the results corresponding with triple collocation. Indeed, both comparisons constitute different metrics that provide different pieces of a complex puzzle.

4 Data availability

The product is available for visualization purposes on the website of the CMEMS Lambda project (2021) (http://www.cmems-lambda.eu/mapviewer/) and through the EMODnet (2021) website (https://www.emodnet-physics.eu/map/Products/Smos/). The access to the data is provided by the Barcelona Expert Center (2007) FTP service; for more details see http://bec.icm.csic.es/bec-ftp-service/. The DOI of the level-3 product is https://doi.org/10.20350/digitalCSIC/12601 (Olmedo et al., 2020a), and the DOI of the level-4 product is https://doi.org/10.20350/digitalCSIC/12600 (Olmedo et al., 2020b).

5 Conclusions

We have presented 9 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 respect to the previous version of these products: (i) a new latitudinal–seasonal debiasing has been included; (ii) improved filtering criteria based on the salinity geophysical variability have been applied, which allows a better description of the salinity gradients without increasing the overall noise error in the maps; (iii) new interpolation schemes are proposed to allow better description of 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; and (v) an estimation of the salinity uncertainty is provided in the new products.

We have performed an extensive validation of the BEC SMOS SSS products v2.0. For doing this, we have compared the 9-year time series of the new BEC SMOS SSS with Argo uppermost salinity, and we have also compared the performance of BEC products with the other three satellite SSS products (the SMOS product produced at CATDS and two SMAP products generated by REMSS and JPL) and the reanalysis product distributed by CMEMS, but in this case restricted to the year 2017. The main conclusions of this comparisons are as follows.

• The statistics of the comparisons with Argo salinity evidence a competitive performance in comparison with the statistics of the rest of the 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 L4) are very similar (being lower than 0.02 psu at a global scale), but the standard deviation is significantly different among them, with the BEC HR being the one with the largest standard deviation (lower than 0.34 psu at a global scale) and BEC L4 the one with the lowest deviation (lower than 0.27 psu).

• In terms of effective spatial resolution and geophysical consistency, we have used two different metrics.

  • Singularity analysis. The SEs of BEC HR and BEC LR SSS products are very similar to the ones of the other satellite salinity products in terms of correlation with OSTIA SST SEs. A clear improvement is observed in the BEC L4 that presents a higher correlation with the SEs of OSTIA SST, suggesting that the geophysical consistency is the most accurate as it is the closest to OSTIA SST SEs. This is partially because of the multifractal fusion method used in the generation of the BEC L4 product.

  • 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 presents similar spectral slopes to the ones of OSTIA SST, showing consistent slopes up to the 50 km wavelength (25 km spatial resolution). At smaller scales, BEC L4 presents evidence of lack of structure and oversmoothing, so it probably does not resolve 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 the geophysical signal consistently.

• 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.