A Multi-Pixel Split-Window Approach to Sea Surface Temperature Retrieval from Thermal Imagers with Relatively High Radiometric Noise: Preliminary Studies

: In the following decade(s), a set of satellite missions carrying thermal infrared (TIR) imagers with a relatively high noise equivalent differential temperature (NEdT) are expected, e.g., the high resolution TIR imagers ﬂying on the future Thermal infraRed Imaging Satellite for High-resolution Natural resource Assessment (TRISHNA), Land Surface Temperature Monitoring (LSTM) and NASA-JPL/ASI Surface Biology and Geology Thermal (SBG) missions or the secondary payload on board the ESA Earth Explorer 10 Harmony. The instruments on board these missions are expected to be characterized by an NEdT of (cid:39) 0.1 K. In order to reduce the impact of radiometric noise on the retrieved sea surface temperature (SST), this study investigates the possibility of applying a multi-pixel atmospheric correction based on the hypotheses that (i) the spatial variability scales of radiatively active atmospheric variables are, on average, larger than those of the SST and (ii) the effect of atmosphere is accounted for via the split window (SW) difference. Based on 32 Sentinel 3 SLSTR case studies selected in oceanic regions where SST features are mainly driven by meso to sub-mesoscale turbulence (e.g., corresponding to major western boundary currents), this study documents that the local spatial variability of the SW difference term on the scale of (cid:39) 3 × 3 km 2 is comparable with the noise associated with the SW difference. Similarly, the power spectra of the SW term are shown to have, for small scales, the behavior of white noise spectra. On this basis, we suggest to average the SW term and to use it for the atmospheric correction procedure to reduce the impact of radiometric noise. In principle, this methodology can be applied on proper scales that can be dynamically deﬁned for each pixel. The applicability of our ﬁndings to high-resolution TIR missions is discussed and an example of an application to ECOSTRESS data is reported.


Introduction
The retrieval of the sea surface temperature (SST) from space is a well-consolidated remote sensing activity, and derived products are used in several operational applications [1].Currently, and for the foreseeable future, global ocean operational SST retrieval is based on measurements from a set of ad hoc instruments, such as VIIRS [2], SLSTR [3] and METIMAGE [4].Generally, SST retrieval involves state-of-the-art atmospheric correction algorithms based on the split window (SW) method (see, e.g., [5,6]), that exploit observations from two or three atmospheric channels and, in some cases (e.g., SLSTR and A-ATSR), from dual view observations [7].Alternatively, SST retrieval is based on the Optimal Estimation Method (OEM) [8].
SST applications require very low uncertainty of the estimated value (<0.1 K), and the requirements for the instruments and processing algorithms to fulfil such a condition are stringent in terms of channel co-registration, radiometric characteristics and the accuracy of the inversion algorithm.
Besides the state-of-the-art instruments dedicated to SST retrieval, other instruments with weaker radiometric and geometrical performance are equipped with channels that can be used for SST estimation.The consequent lower quality of the retrieved products may be derived by the secondary payload nature of the instrument (e.g., CALIPSO-IIR [9] and HARMONY thermal infrared (TIR) payload [10][11][12]) or by the high spatial resolution (HR) of the instrument (e.g., ASTER [13], ECOSTRESS [14] and LANDSAT8-9 [15,16]).In particular, within the next few years, there are several satellite missions planned to launch HR TIR sensors that, although mostly dedicated to land applications, will acquire data over coastal regions at the very least.Obtaining a high-resolution SST over coastal regions is of interest for coastal process studies, for example, the study of the dynamics of upwelling systems [17], inshore water quality [18] and the thermal discharge of cooling water from nuclear power stations [19].
The aforementioned high resolution missions are: • The NASA-JPL/ASI Surface Biology and Geology Thermal Mission (SBG), focusing on five research and applications areas: terrestrial and aquatic ecosystems, hydrology, weather, climate and solid Earth [20].

•
The Thermal infraRed Imaging Satellite for High-resolution Natural resource Assessment (TRISHNA) mission, a cooperation between the French (CNES) and Indian (ISRO) space agencies to launch a satellite in 2025 with a 5-year lifetime to measure the visible, near infrared and thermal infrared signal of the surface atmosphere system globally approximately twice a week, with a 60 m resolution for the continents and coastal oceans and a resolution of 1000 m over deep oceans [21].

•
The Land Surface Temperature Monitoring (LSTM) is an ESA mission set to join the Copernicus Sentinel system in 2028.The satellite will have TIR observation capabilities over land and coastal regions in support of agriculture management services, and possibly a range of additional services.The LSTM mission will consist of two satellites [22].
Table 1 reports the relevant characteristics of a set of past, current and future sensors with relatively high noise equivalent differential temperature (NEdT) in terms of the spectral position of the SW channels and the relevant estimated NEdT values.As a reference, the characteristics of the Sentinel 3 (S3) SLSTR equivalent channels are reported.It should be noted that, in the best of these cases, the NEdT of the reported instruments, excluding the reference SLSTR, is about 0.1 K.
Over oceans, the SW term is proportional to the total amount of water vapour in the atmosphere.Roughly, 1 cm of total precipitable water vapour in the sensed atmospheric path corresponds to 0.5 K of SW difference [23].For this reason, in conditions of a very dry atmosphere, such as in polar regions, the SW may be dominated by noise, and a single channel correction may be more effective [24].Besides the extreme case of polar regions, in general, for a given scene (i.e., a geographical region in a given day), the dynamical range of the SW difference is, on average, expected to be around 1 K. Within this frame, it is clear that a NEdT of the order of 0.1 K or higher for the two channels involved in atmospheric correction may induce artefacts in the retrieved SST mostly in the form of a noisy field.As an example, assuming uncorrelated noise, the random uncertainty associated with the SW term can be up to about 0.42 K if both channels have an NEdT of 0.3 K, which is about 40% of the expected range for a given scene.
Since for most of the aforementioned missions the objective is not absolute SST retrieval but rather the identification of structures in the SST field (e.g., gradients and patterns), a noisy SST field results in difficulties in identifying the actual SST patterns.Rather than reducing the noise by smoothing the SST field or using noise robust operators to identify spatial patterns, losing in both cases spatial information, we suggest to explore the possibility of reducing the noise in the atmospheric correction term.This is based on the following hypotheses:

•
The atmospheric correction is based on a functional form that uses the SW difference; • The scales of variability of the radiatively active atmospheric components are larger than the spatial sampling distance of the sensor.
A similar approach has been already presented and discussed for SST retrieval by several authors and applied to different instruments.
To reduce the noise in the AVHRR-derived SST based on an SW algorithm, it has been suggested to rely on the average or the median of the SW difference computed over a 3 × 3 [25], 5 × 5 [26] or up to 11 × 11 [27] pixel area.
A similar solution is suggested for SST retrieval from SEVIRI ( [28]), discussing the effect of using the average of the SW difference over a 5 × 5 pixel box.In the same study, a strategy to introduce the concept of spatial smoothing also in case of retrieval based on optimal estimation was presented and discussed.
A different approach was suggested in [29,30] for ATSR.In particular, Harris and Saunders [30] proposed to smooth the single pixel ATSR-derived SST (SST k,m ) by considering an area around the pixel in the following way: where BT11 i,j is the 11 µm brightness temperature (BT) of pixel i, j and G i,j is a clouddiscriminator operator and adopts a value of 0 when the pixel i, j is cloudy and 1 if it is cloud free.The smoothing is estimated as the average of the difference in the single pixel ATSR-derived SST (SST i,j ) and the corresponding 11 µm BT (BT11 i,j ).The authors applied smoothing up to a 29 × 29 pixel area (of 1 km size), suggesting that in the above 19 × 19 pixel area there are no significant improvements in noise reduction.The approach described in Equation ( 1) was also applied to the SLSTR SST-derived L2P distribution, applying a 3 × 3 pixel kernel (Sea Surface Temperature (SLSTR) Algorithm Theoretical Basis Document, SLSTR-ATBD L2SST-v2.5,October 2012, https://www.eumetsat.int/media/40007,accessed on 29 January 2023).
Finally, the assumption of a relatively low spatial variability of atmospheric properties with respect to surface properties is commonly adopted for similar scales by wellestablished inversion algorithms for other geophysical variables, for example, for aerosol property retrieval by the Generalized Retrieval of Atmosphere and Surface Properties (GRASP) [31].
This study, independent from the SST retrieval algorithm used in the sensor, focuses on determining under which conditions and scales it is possible to apply an atmospheric correction based on a spatial average of the SW difference to reduce the noise on the final SST field.This will be done by means of real observations from SLSTR instruments on board the Sentinel 3A and 3B satellites.
In order to document the limits of applicability of the proposed approach, a set of SLSTR cases over oceanic regions characterized by the presence of high SST spatial variability was selected.In principle, in the presence of different water masses, as, for example, in the Gulf Stream area, it is likely that an enhanced spatial variability of the lower troposphere results from different air-sea exchanges.
Section 2.1 contains a description of the SLSTR dataset used in this study.In Section 2.2, we present a short review of SST retrieval based on the SW difference with some references relative to the missions of interest.We go on to illustrate the evidences/findings of the proposed methodology in Sections 3 and 4. Finally, discussions of the results and the main conclusions are provided in Section 5.

SLSTR Study Cases Data-Set
For our investigation, we relied on the S3 (A and B) Non Time Critical products distributed as L2P SST products by EUMETSAT.These data are compliant with the specifications from the Group for High Resolution Sea Surface Temperature (GHRSST, https://www.ghrsst.org/).They include: (i) L1 single channel Top Of the Atmosphere (TOA) BTs and associated NEdT; (ii) skin SST L2 data and (iii) specific data flags to identify the type and quality of the observations (detailed below).The BTs are provided as observations obtained from the three S3 IR channels S7, S8 and S9, centered at 3.74, 10.85 and 12.02 µm, respectively.
The l2p_flags variable contains bitmasks that enable the identification and removal of land, ice, lake and river pixels from the L2 data.The quality_level variable which ranges from 0 (no data) to 5 (best quality) was used to identify the best quality (quality_level = 5) pixels.We thus selected BT pixels based on these criteria in order to analyze only high quality, non cloud-contaminated BT observations.These products maintain the original observation geometry of Level 1, i.e., a 1 km Spatial Sampling Distance (SSD) at nadir.

Selection of Study Cases
We focus on oceanic regions where SST features are mainly driven by the meso to submesoscale horizontal/vertical dynamics, thus selecting areas characterized by enhanced eddy activity (i.e., high eddy kinetic energy) [38,39] or upwelling systems [40,41].This led to the identification of the following areas Of interest (AOI) (Figure 1): Such AOIs are often characterized by the presence of fronts due to meandering of the marine surface currents as well as the presence of eddies and filaments [1].Such a scenario allows us to test the methodology proposed in this study in areas where frontal structures may trigger air-sea interactions, making the SW correction term non-trivial [42].A visual inspection of each study case was required in order to assure the presence of at least 50% contiguous cloud-free pixels over zones with recognizable frontal/eddy activity.
For each AOI, we selected four study cases in order to account for the different phases of the annual cycle.Table 2 provides detailed information on the set of study cases used in our work.
Definition of the areas of interest used to select the SLSTR study cases.

Split Window Method
Over oceans, the TOA signal measured from a satellite sensor in the spectral region of interest is, in absence of clouds, dominated by the upwelling radiance emitted by the sea surface, with a contribution due to absorption/emission of the atmosphere.In the majority of cases, the effect of the atmosphere attenuates the radiation emitted by the surface.The SW method is based on different atmospheric absorptions, and therefore on the different measured BTs, in two adequately selected spectral intervals (practically spectral channels) in the IR windows.Typically, the wavebands used for the two channels are in the range of 10.3-11.3µm and 11.5-12.5 µm, respectively [43].The true SST can be estimated using this difference in an atmospheric correction algorithm, the simplest form of which is: where BT 11 and BT 12 are the TOA BTs in the ranges 10.3-11.3µm (less contaminated by atmospheric disturbances) and 11.5-12.5 µm , respectively, and a i are regression coefficients to be determined.More accurate algorithms take into account other terms such as observation geometry, first guess SST values and/or water-vapour-dependent terms (e.g., [44][45][46]).SW-like algorithms developed for SST-dedicated missions (e.g., AVHRR, SLSTR and A-ATSR) have also been used for SST retrieval with high resolution instruments (e.g., ASTER [47], LANDSAT-8 [48] and ECOSTRESS [49]).Adopting these functional forms, the effective reduction in the retrieval noise by smoothing of the SW term depends upon the specific form of the retrieval equation [50] and the relative magnitudes of the noise levels in the sensor wavelength bands.For example, in the MultiChannel SST algorithm: where S θ = (sec(θ v ) − 1), in which θ v is the view zenith angle to account for the observation geometry.If the noise in the two SW channels is uncorrelated, then the noise of the SW term will be larger than the noise of either channel.If the regression coefficient a 2 is >1 and if a 3 S θ > 1, then the SW difference term will be amplified by the retrieval equation and smoothing will reduce the retrieval noise to a level close to the a 1 BT 11 noise.

Evidence for the Multi-Pixel Approach
To reduce the impact of radiometric noise on the retrieved SST, we explored a multipixel approach based on the assumption that the spatial variability of the atmosphere is lower than that of the surface [51,52].
In order to implement such an approach in SW retrieval, we introduce an additional assumption: the atmospheric effect is fully represented by the SW BT differences in the proposed SST retrieval functional form.As an example, Figure 2 shows for a SLSTR Sentinel 3 (S3)-A acquisition over the Kuroshio current area.

•
The S8 channel TOA BT field, as a proxy for the SST field; • The difference between the S8 and S9 channel-derived BTs corresponding to the SW difference (∆BT hereinafter); • The standard deviation over 3 × 3 pixels of the ∆BT;

•
The ratio between standard deviation of the ∆BT and the associated noise.The latter is estimated assuming uncorrelated random noise and using the TOA BT NEdT field for each pixel/band, included in the L2 product [53].
The S8 channel BT exhibits a signature of meso to sub-mesoscale activity, as expected for a western boundary current such as the Kuroshio current.The scene is indeed dominated by mesoscale to sub-mesoscale eddies and filaments, with clearly visible signatures in the SST and hence the TOA BT field (Figure 2a).The ∆BT term still shows large scale features but appears more noisy, given also the lower range of variability (Figure 2b).The 3 × 3 standard deviation of the ∆BT is on the order of tenths of mK (Figure 2c) and exhibits a local variability of the same order of magnitude of the expected variability due to radiometric noise (Figure 2d).This condition mainly fails in areas characterized by spatial gradients in the original BT signal.Over different water masses, it is likely that the air mass can also be different, particularly on the scale of the example, which partly justifies the occurrence of relatively high values of the standard deviation in the BT difference field.On the other hand, for strong SST gradients, such as the one in the example, spatial features can also be generated by the inter channel co-registration uncertainty that, in the system requirements for SLSTR, should not exceed 0.1 SSD (https://sentinel.esa.int/documents/247904/4755279/S3.PN-SLSTR-L1 .09+-+i1r1+-+SLSTR+L1+PB+SL__L1_.004.04.00.pdf).In the case shown in Figure 2, it is thus possible to assume atmospheric homogeneity over a 3 × 3 pixels box, corresponding roughly to an area of about 3 × 3 km 2 .The dimension (in pixels) and shape of the area used to estimate the atmospheric correction should, however, also be a function of the specific scene/sea surface conditions.Corresponding to intense seawater thermal gradients, modulations in the lower atmospheric responses are indeed expected [42,54,55].This additional investigation is however left for future studies.

Results
Based on the study cases detailed in Table 2, we quantified the validity of our assumptions via statistical and spectral analyses.
Figure 3a shows the histogram distribution of the ratio between the local variability of the ∆BT term (∆BT STD) and the associated noise (NEdT).As already mentioned in Section 2, the statistics uniquely rely on cloud-free sea pixels.However, to avoid artefacts due to cloud/land contamination, all sea pixels contiguous to a cloudy or land area have not been accounted for.The histograms, reported as normalized frequencies per AOI, peak around the unity value, supporting the hypothesis that the local variability of the atmospheric correction term (over the 3 × 3 km 2 area) is dominated by radiometric noise.
For some regions (e.g., EAC), the distribution still peaks around unity, although a slight shift toward larger values can occur.This is discussed on a test case in the Gulf Stream area, acquired by Sentinel 3B on 30 March 2021 (Figure 4c, file S3B_SL_2_WST____20210330T140133_20210330T154232_202104 01T003643_6059_050_352______MAR_O_NT_003), but it was found also for other AOIs such as the EAC.
For this specific scene, the ratio between the 3 × 3 variability of the ∆BT and the associated noise is around 1 north of the main SST front located between 39 • N and 40 • N, while it increases to 2 south of it (Figure 4a).A closer look at the ∆BT noise field (Figure 4b) [53], as well as the BT from the S8 channel (i.e., our best proxy for SST), evidenced a distinct behavior between the colder and warmer areas in the analyzed test case.In the presence of warmer SSTs, the overall ∆BT noise is reduced and, simultaneously, one can expect a higher atmospheric variability, i.e., an enhancement in the ∆BT standard deviation [42].These effects can thus contribute to an overall increase in the ∆BT standard deviation/∆BT noise ratio.This can also explain the peak displacement of the normalized histograms presented in Figure 3a with respect to unity.However, this is not expected, or at least, not to the same extent, in very high resolution missions such as, e.g., ECOSTRESS, where the NEdT at high temperatures is expected to be significantly larger than for SLSTR.The assumption of atmospheric homogeneity may fail if we account for areas larger than 3 × 3 km 2 .We thus analyzed a specific case study over the Kuroshio current on the 19 September 2021 (file S3A_SL_2_WST____20210919T232056_20210920T010155_20210921T1 02504_6059_076_272______MAR_O_NT_003).
By imposing homogeneity over a 5 × 5 km 2 area, the distribution of the ∆BT standard deviation divided by the ∆BT NEdT becomes indeed significantly broader and it is not centered around unity, as shown in Figure 3b.
In addition, to further support the assumption on the basis of the multi-pixel approach, we compared the power spectral density (PSD) of the S8 channel BT, the PSD of the ∆BT and the PSD of the ∆BT term averaged by means of a 3 × 3 moving window.The PSD was computed via fast Fourier transform, relying on a Blackman-Harris window to reduce spectral leakage, following [56].An example of S3A acquisition over the Kuroshio current area (already discussed in Section 3) is depicted in Figure 5.In particular, the PSD was computed over three cross-sections (white thick lines in Figure 5a).The cross-sections cover approximately 5 longitude degrees and 0.1 latitude degrees, and were chosen to perform the analyses over distinct dynamical regimes of the selected case study:

•
A region with predominance of large-scale patterns at the westernmost and easternmost zones of the scene (cross-section 1); • An area with tiny filaments characterized by sharp temperature differences, such as the one located at 149.7 • E−48.0 • N (cross-section 2); • An area with much smoother spatial variability (cross-section 3).
The results reported in Figure 5b-d were computed assuming a uniform pixel spacing of 0.015 degrees, the latter was determined as the mean longitudinal spacing in the areas interested by the spectral analysis.Apart from small differences, the three sets of transects yield the same results.The S8 BT (BT8) spectrum exhibits PSDs larger than 10 2 for wavenumbers smaller than 1 deg −1 (scales 100 km).The PSD progressively decreases for larger wavenumbers and eventually exhibits a flat spectrum, indicating noise, for scales 2.5 km.The ∆BT (or SW) term PSD is fully compliant with the BT8 one until spatial scales of 10 km; the noisy regime is observed immediately afterwards.Such a behavior suggests that the ∆BT term prevents a clear description of spatial features below 10 km, injecting noise even in the retrieval of SST based on SW-like algorithms.
However, considering the findings reported in Section 3 and summarized in Figure 3, we can perform the same spectral analysis after applying a moving 3 × 3 averaging window to the original SW field, with the aim of double checking whether this operation reduces the noise affecting the SW term.This operation makes the ∆BT field spectrally compliant with the S8 channel, reducing the amount of noise in the 10 km to 3 km range.For lower scales, the PSD of the 3 × 3 averaged ∆BT term (∆BT 3 × 3 hereinafter) exhibits an abrupt decrease and a subsequent increase, likely due to artefacts introduced by the 3 × 3 averaging.In future studies, this can be further investigated by choosing different averaging kernels.However, this is not expected to affect our analysis, since this occurs in a spectral region where noise dominates both the SW and the S8 field.

Discussion and Conclusions
This study approaches the issue of reducing the noise in the atmospheric correction term in multi-spectral regression IR-based SST retrieval algorithms.This approach has been presented and documented by several authors in the past [25][26][27][28][29][30] for state-of-the-art SST-dedicated satellite missions.Within the frame of present day or planned missions with a relatively high NEdT, this methodology can be useful to reduce the noise in the retrieved SST field with advantages for the effective resolution of the estimated SST. Characterization of the main SST frontal features may benefit from such an approach.
The study is based on the hypothesis that the spatial variability of relevant atmospheric variables (mostly temperature and water vapour distribution in the troposphere) is lower than that of the sea surface temperature.In this study, we documented the hypothesis behind the smoothing of the atmospheric correction term by using Sentinel 3 (A/B) SLSTR acquisitions over regions mainly dominated by mesoscale/sub-mesoscale turbulence (e.g., western boundary currents) or those corresponding to major upwelling systems.The aim of this choice was to test the methodology over SST scenes with complex spatial patterns and enhanced frontal activity.The dataset was chosen in order to sample the annual cycle and is characterized by a spatial resolution of 1 km.Relying on statistical and spectral analyses, we found the following main results:

•
The local variability of the atmospheric correction term generally employed in splitwindow algorithms (here referred to as ∆BT), if computed over 3 × 3 km 2 areas, is of the same order of magnitude of the estimated ∆BT radiometric noise; • Spectral analyses of the BT, ∆BT and 3 × 3 km 2 averaged ∆BT (performed over a specific study case) suggest that 3 × 3 smoothing on the SW term is helpful in reducing noise in the atmospheric correction procedure and obtaining a SW correction with similar spectral properties to the BT derived from channel S8.This is verified up to scales larger than 3-4 km.It is thus very likely that the SST obtained with the 3 × 3 averaged SW will benefit in terms of noise reduction/effective resolution; • The assumption of atmospheric homogeneity over the 3 × 3 km 2 area for the SLSTR case studies was thus supported by our findings.The horizontal extent of the atmospheric homogeneity assumption is, however, a function of the specific satellite mission and should be adapted as a function of the IR sensor SSD and radiometric characteristics.
As shown in some of the examples, there are cases where the local standard deviation of the SW exceeds the noise value expected for the difference.Besides the contamination by clouds and coastlines, which should be reduced by the applied masks, there are other two factors that could be causing this feature observed mostly corresponding to SST gradients:

•
A realistic difference due to the atmospheric response to SST gradients; • An artefact due to inter-channel co-registration issues; In our analyses, the adopted macro-pixel was always a 3 × 3 box independent of the effective spatial variability of the BT fields.The technique should be optimized, particularly for the HR satellite missions (e.g., ECOSTRESS, ASTER, SBG, THRISNA and LSTM).
ECOSTRESS is an example of a high spatial resolution instrument with a high noise level.We used ECOSTRESS L1B Geolocation [57] and L1B Radiance and L2 LST [58] products produced by the NASA Jet Propulsion Laboratory and distributed by the Land Processes Distributed Active Archive Center (LPDAAC, https://lpdaac.usgs.gov).Conversion of radiance to BTs and re-projection of the swath data onto a rectilinear longitude-latitude grid was carried out with the ECOSTRESS_swath2grid Python script provided by LP.DAAC.We analyzed a case study over the Alboran Sea (region ALB in Figure 1) (files listed below): We used a collocated VIIRS SST retrieval (file 20210830015000-OSPO-L2P_GHRSST-SSTsubskin-VIIRS_N20-ACSPO_V2.61-v02.0-fv01.0.nc) reprojected onto a rectilinear longitude-latitude grid with SeaDAS 8.2.0 software.The ECOSTRESS channel wavelengths were 10.49 µm and 12.09 µm, which we used for the BT 10 and BT 11 terms in the SST retrieval equation (3).To establish regression coefficients for the MC SST retrieval algorithm, the 60 pixel regional averages of ECOSTRESS BTs and view zenith angle fields were down-scaled to the VIIRS coordinates with a k-nearest neighbor algorithm.Smoothing was performed by convolving the SW (∆BT) fields with 3 × 3 ( 200 × 200 m), 11 × 11 ( 770 × 770 m) and 51 × 51 ( 3.57 × 3.57 km) averaging kernels.BT noise (robust standard deviation), retrieval equation terms, and the retrieved SST were measured with the von Neumann sequential difference method [59].The noise levels were calculated separately for all rows and columns of the data fields and the median of the row values and column values was used as per image noise estimates.The regression of VIIRS SST against ECOSTRESS predictors was highly significant (R 2 = 0.929, RMSE = 0.475, p < 0.001, df = 3.88327).The coefficients are listed in Table 3.The coefficients associated with the SW terms in the equation (a 2 and a 3 ) were >=1, indicating the potential for amplification of noise by the SW terms during retrieval of the SST [50].The noise level of the 12 µm BT was approximately 3 times greater than the noise level of the 11 µm BT (Table 1).Smoothing greatly decreased the noise level of the SW ∆BT (Table 4 and Figure 6), but had little effect on the noise level of the retrieved SST.This occurred because the noise level of BT11 µm was much higher than the noise level of the smoothed ∆BT.The high BT11 µm noise level in ECOSTRESS is largely systematic rather than Gaussian (see Figure 7), and derives from stripes in the images caused by calibration artefacts and non-uniformity of the focal plane array.3) and for the retrieved ECOSTRESS SST.The values refer to the case shown in Figure 6.

Equation Term
No However, depending upon the form of the SST retrieval algorithm and its coefficients, in sensors with high noise levels in the primary SST channel (10-11 µm), smoothing of the split-window terms may not always reduce the retrieval noise.In these cases, other noise mitigation approaches are necessary, including de-striping and direct smoothing of the retrieved SST.
For such missions, in order to obtain an effective noise reduction in the retrieved SST, smoothing should be applied not only to the SW term but also to the BT 11 term.However, based on the assumption of the different scales of variability for oceans and the atmosphere, BT 11 smoothing should be performed over smaller areas, with the objective of maintaining the highest effective resolution possible for the retrieved SST.
The scales of smoothing for each specific sensor and scene should be a trade-off between (i) the noise associated with the measurement, (ii) the sensor SSD and (iii) the target scene effective spatial variability.
Finally, in this work, we adopted a simple square smoothing box.The extent and shape of the smoothing area can be dynamically defined on the basis of the characteristics of the observed field (e.g., [60]).
Finally, it should be noted that for each pixel of the area used to average the SW difference, it is possible to retrieve an SST value.By moving the area for a given pixel, on can obtain different estimates of the SST.For example, in the 3 × 3 pixel cases, up to nine SST values can be obtained for a single pixel.The averages and standard deviations of the SST computed over the set of values relative to a single pixel can be further used as indicators of the robustness of the procedure.In principle, if the hypotheses are correct, the standard deviation of the estimated SST should be on the order of magnitude of the expected noise, as estimated by uncertainty propagation.Deviations from the expected behavior should be analyzed to optimize the procedure.

Figure 3 .
Figure 3. (a) Distribution of the ratio between the 3×3 pixel standard deviation (STD) of ∆BT and the associated noise.The analyses are reported per each AOI of Figure 1.(b) Distributions of the ratio between a 3 × 3 (green) and a 5 × 5 (blue) ∆BT STD and its associated noise.Example for the study case over the KUR AOI on 19 September 2021 reported in Figure 2.

Figure 4 .
Figure 4. (a) Ratio of the 3 × 3 pixel standard deviation (STD) of the ∆BT term and the associated noise.(b) ∆BT noise, expressed as NEdT [53].(c) BT acquired from the S8 Channel.The data refer to an acquisition from Sentinel 3B on the 30 March 2021 at 14:01 UTC.White areas stand for low quality pixels.

Figure 5 .
Figure 5. (a) TOA BT (S8 channel) field over the Kuroshio area on 19 September 2021 at 23:20 UTC, acquired by Sentinel 3A.(b-d) Power spectral densities of TOA BT S8 (blue), Split Window (SW) term (red) and SW averaged over a 3 × 3 pixel area for the cross-sections 1-3 reported in panel (a).The PSD is presented as a function of the wavenumber, the corresponding spatial scales are also reported.

Table 4 .
Noise levels (K) for the geometric and split-window (∆BT) terms of Equation (

Table 1 .
Examples of spectral and radiometric (NEdT) characteristics of IR channels commonly used for SW computations of interest for this study.

Table 2 .
Summary of good quality pixels used in the analyses for each AOI, including specific information of the data source: day (D), month (M), year (Y), Universal Time Clock hour (UTC-H), minute (m) and Sentinel 3 platform (A or B).