Improving the combined use of reflectance and thermal channels for ocean and coastal cloud detection for the Sea and Land Surface Temperature Radiometer (SLSTR)

Reflectance imagery is used to aid daytime cloud detection in thermal remote sensing. This paper presents new approaches to utilization of reflectance for sea surface temperature (SST) remote sensing for the Sea and Land Surface Temperature Radiometer (SLSTR). SLSTR is an along-track scanning sensor with a complex instrumental viewing geometry, and no co-registration of the fields-of-view of reflectance and thermal pixels. Reflectance channels have twice the spatial resolution of thermal channels, and observations are placed on compatible “ image grids ” for the convenience of users of Level-1 data. We highlight limitations of simple methods, based on these image grids, of using reflectance imagery to inform cloud detection at the thermal resolution. We present improvements from averaging the N-nearest reflectance observations directly to the infrared instrument geometry, where N = 10 is chosen in this study when using the A and B stripes together. We show that the standard deviation of the N-nearest reflectance observations is another calculable quantity of use to improve the discrimination of clouds in the infrared image, beneficially reducing the weight placed on coarser-scale thermal spatial variability. The developments are illustrated by case studies in coastal zones over optically bright waters and around strong ocean fronts, and the benefit for SST products is quantified by the impacts on coverage and validation statistics. In a case study over optically bright waters, the clear-sky fraction increases from 46.2% to 93.1%. Coastal zone validation shows a 23.8-31.6% reduction in the false alarm rate and a corresponding 27.1-33.3% increase in the cloud detection true skill score. Globally the robust standard deviation for clear-sky matches between SLSTR-A and drifting buoys reduces from 0.3 to 0.29 with a 6% reduction in data due to improved screening of scattered cloud.


Introduction
The Sea and Land Surface Temperature Radiometers (SLSTRs) provide high-resolution thermal-infrared observations of the Earth's oceans and are essential to the continuity of high-quality sea-surface temperature (SST) climate data records (CDRs) (Merchant et al., 2019).SLSTR is a dual-view instrument flown aboard the Sentinel 3 satellites, making observations in a nadir scan and a rear-facing scan (ESA, 2022).The SLSTR instruments extend the data record of the Along-Track Scanning Radiometer (ATSR) sensor series, from which sea-surface temperatures can be retrieved with biases <0.1 K (Embury and Merchant, 2012).Two Sentinel 3 satellites are currently operational, flown in tandem with a revisit time of 0.9 days at the equator (0.8 days for latitudes >30 • ) (Donlon et al., 2012).Global ocean coverage at daily resolution is important for operational oceanographic monitoring and meteorological applications (Bonekamp et al., 2016;Coppo et al., 2010).Near real time (NRT) data products with a timeliness of three hours can be used operationally for oceanographic forecasts and numerical weather prediction (Bonekamp et al., 2016;Good et al., 2020).
Retrieval of SST from satellite data requires a pre-processing step to identify and remove cloud-contaminated observations that otherwise cause erroneous retrieved-surface-temperature (Merchant et al., 2005(Merchant et al., , 2019;;Závody et al., 2000;Simpson et al., 2001).Cloud detection algorithms work on the premise that clouds are typically 'brighter' and 'colder' than the underlying ocean surface (Merchant et al., 2019;Závody et al., 2000;Frey et al., 2008;Rossow and Garder, 1993).Clouds have different spectral characteristics in both reflectance and emissivity to the ocean and different spatial and temporal variability.For daytime SST retrievals, the best cloud detection results are often achieved when using a combination of observations at reflectance and infrared wavelengths (Merchant et al., 2019;Frey et al., 2008;Rossow and Garder, 1993;Bulgin et al., 2018).In coastal regions, cloud detection needs to account for conditions not seen over the open ocean, such as water that is optically brighter from suspended sediment or biological matter, as well as more turbid (Wang et al., 2020;Lu et al., 2021).
The SLSTR instrument reflectance channels have a higher spatial resolution (500 m at nadir) than the channels at infrared wavelengths (1 km).Measurements are mapped from the measurement grid defined by the conical scan of the SLSTR instrument to an image grid which is regularly spaced with respect to the satellite track (ESA, 2022).While the image grid definition is in common, the transformation of the infrared and reflectance channels is undertaken independently.Some applications such as cloud detection require the use of the reflectance data on the same image grid as the infrared data.It is commonly assumed that the higher-resolution reflectance data, maps directly to the infrared image grid, using a simple averaging of every 2 × 2 pixels (Coppo et al., 2010).However, due to the instrument viewing geometry and prior translation of the data from the measurement to image grids, this approach often results in sub-optimal spatial reconciliation of the reflectance and infrared observations, as shown below.
Any inconsistency in the mapping to the image grid of reflectance and infrared observations will have a direct impact on cloud detection skill.This is particularly important in regions with strong gradients e.g.ocean fronts, land-sea boundaries or cloud edges, where inconsistently mapped pixels may view different features.In this paper we present a new SLSTR preprocessor to map reflectance channel observations more closely to the location of the infrared data, and assess the impact of this improved co-location on cloud detection skill.In the process of mapping the reflectance data to the infrared observations, we also calculate the sub-pixel reflectance variability, which can be used as an additional metric in the cloud detection scheme.
The remainder of the paper is structured as follows: in section 2 we describe the characteristics of the SLSTR instrument and in section 3 we present the new SLSTR preprocessor for co-locating infrared and reflectance channel information.In section 4 we describe the Bayesian cloud detection and introduce the new sub-pixel reflectance variability metric.The benefit of the new preprocessor and sub-pixel reflectance variability to cloud detection is assessed in section 5 along with the presentation of global validation statistics.We discuss the outcomes and conclude the paper in section 6.

Imagery from SLSTR
SLSTR is a dual-view instrument flown aboard the Sentinel 3a and Sentinel 3b platforms.Observations are made at nine different wavelengths covering the visible (VIS, S1-3), shortwave infrared (SWIR, S4-6) and infrared (IR, S7-9) parts of the electromagnetic spectrum (Table 1).Two additional channels are available at infrared wavelengths, with an extended brightness temperature range for fire monitoring (F1-2) (ESA, 2022).VIS and SWIR channels have a higher spatial resolution of 500 m compared with 1 km for IR observations.SLSTR has a wider nadir swath than the earlier ATSR instruments (ATSR-1, ATSR-2 and AATSR), at 1470 km compared with 512 km in order to increase data coverage, and the dual-view swath width is 740 km (ESA, 2022).The IR channels are sensitive to surface temperatures in the absence of clouds and the VIS and SWIR channels are present to aid cloud detection during the day.
SLSTRs have synchronised conical scanners making observations in both a nadir and rear-facing (oblique) view (with satellite viewing angles of no greater than 55 • ) (Coppo et al., 2010).Every point within the oblique view swath is therefore observed twice with two different atmospheric paths.This "dual-view" capability is useful in reducing uncertainty when retrieving SST (Coppo et al., 2010).This dual-view concept is specified to enable SST retrievals with zero bias and an accuracy of <0.3 K (Donlon et al., 2012).
Fig. 1a shows the across-track viewing geometry for the nadir and oblique views.The nadir swath is asymmetric around the sub-satellite point extending 472 km on one side and 935 km on the other (Coppo et al., 2010).This offset swath provides full overlap with observations from the Ocean and Land Colour Instrument (OLCI), which is also carried on Sentinel 3 platforms.The oblique view swath is narrower with a constant satellite zenith angle.
SLSTR observations made in each of the two views follow the conical scan of the instrument and form the measurement grid (where observations are irregularly spaced).These are then mapped onto an image grid, regularly spaced with respect to the satellite track, and form the basis of the Level-1 data product.Where there are no corresponding observations on the measurement grid, or the data are invalid, the image grid pixel is cosmetically filled with duplicate data from the nearest of the eight surrounding pixels (ESA, 2022).There are also some data points on the measurement grid that are not used in the mapping to the image grid (due to the grid alignment they may not be the closest point to any image pixel); these are provided in the SLSTR level-1 product as an additional "orphan" pixel array.
In the IR (S7-S9), data are provided on a regular 1 km grid known as the 'I stripe' (Table 1).The higher-resolution VIS and SWIR channels have two detectors providing stripe A and stripe B data.Channels S1-S3 have stripe A data only, with nominally two observations corresponding to each I stripe pixel.Channels S4-6 have both A and B stripe data with four observations corresponding to each I stripe pixel.A, B and I stripe data are all mapped independently from the measurement to image grids.
Fig. 1 shows the differing alignment of the of the I, A, and B stripe observations on the measurement grid as a function of instrument view and swath location.The footprints of the SLSTR pixels on the ground are defined by the instantaneous field of view (IFOV) of the detectors, propagated using a model along the instrument line of sight (Cox et al., 2021).The IFOV map for each detector, with the scanning mirror stationary, was measured during the ground testing of SLSTR before launch with a resolution of 10 ′ ′ (corresponds to approximately 40 m).The projected footprint depends on both the position and motion of the scanning mirror during the detector integration time for each pixel.Due to the conical nature of the SLSTR scan, the footprint shape is stretched and increases in area towards the edge of the swath, and so must be calculated separately for every pixel.The footprint map was calculated from the 2D autocorrelation of the projected IFOV during the detector integration at each pixel position.
In the SLSTR Level-1 (L1) product, each pixel has an associated pixel number determined by its position in the swath.To convert between this pixel number and a footprint shape, the shape at four reference pixels across the swath is approximated using 6 vertices on the 10% contour of the footprint map.The positions for these vertices were chosen to give a good visual fit for all the detectors in the channel.The shape of the remaining pixels is found by interpolating between the vertex positions for the four reference pixels.
Footprint vertices were generated for channels S6 and S7 and these were assumed to be representative of all other VIS, SWIR and IR channels respectively.Channel S6 was chosen to represent the VIS and SWIR channels as it has a slightly broader footprint than the VIS channels, and thus the 10% contour for channels S1-S3 is enclosed within the set of vertices picked for S6.To determine the position of the vertices on the ground in latitude and longitude coordinates, it is necessary to find the angle between the along-track direction and north as this sets the orientation of the footprint.This was done by calculating the bearing between points thirty pixels above and thirty pixels below the pixel in the instrument grid.The latitude and longitude of the centre of each pixel are recorded in the L1 product.Combining this information with the orientation of the pixel and the shape of the footprint (as determined by the pixel number) enables calculation of the latitude and longitude of the vertices on the ground.
Subplot 1b shows the nadir view, close to the sub-satellite point.The large blue dot denotes the centre of the infrared pixel, and the blue line indicates the extent of the I stripe pixel IFOV.Pixel centres for observations from the A and B detectors are shown in orange and green respectively.Under this viewing geometry, around six pixel centres from both the A and B detectors fall within the IR IFOV, and the A and B observations are well aligned.Subplot 1c shows the oblique view, for the same part of the satellite swath.The increased satellite zenith angle enlarges the IR pixel footprint, with 12-15 observations from the A and B detectors falling within the IR IFOV.At the edge of swath for the nadir view (subplot 1d) there is a similar enlargement of the IR IFOV and more irregularity in the spacing of the A and B detector pixel centres.
Given the nature of the measurement grid as shown in Fig. 1, independent nearest-neighbour mapping of the A, B and I observations to two compatible image grids at different resolutions can lead to observations appearing in the same location that only partially overlap in reality.Using a simple averaging to downscale the VIS channel observations to the IR resolution can give rise to inconsistencies in the geolocation between channels, particularly along feature edges e.g.coastlines, cloud edges.This is demonstrated in Fig. 2 in imagery along the Norwegian coastline.Subplot 2a shows the S7 infrared channel, where the leftmost coastline appears to have a saw-tooth edge which is caused by cosmetic fill to locations over the ocean from nearest neighbours that are really over the land.The higher resolution S6 channel imagery (subplot 2b) verifies that the coastline is reasonably straight in the direction of the satellite track.Subplot 2c shows the outcome of coarsening the S6 data in the obvious manner using the image grids: 2 × 2 averaging of the reflectance measurements that in the image grid correspond to the infrared location.Along the coastline we see disparity in the signal and apparent shape of the coastline between the S7 and coarsened S6 data due to the independent mapping of each channel from the measurement to image grids, prior to coarsening the S6 data.While 2 × 2 averaging is the obvious approach given the image grids provided, the results are not optimal for joint use of infrared and reflectance imagery.No knowledge of cosmetically filled pixels in the IR is taken into account, and orphan pixels are excluded.

Reflectance imagery optimised for joint use with thermal channels
The spatial mismatch between the VIS and IR footprints that arises when coarsening the VIS data on the image grid can be demonstrated using the pixel IFOVs.Fig. 3 shows the two nadir view IR IFOVs from Fig. 1, one at the sub-satellite point (3a) and one at the edge of swath (3d).The corresponding VIS pixel IFOVs obtained when selecting 2 × 2 pixels from the VIS image grid to be averaged, are shown in orange.For the nadir pixel at the sub-satellite point, we see that the VIS channel data is offset from the centre of the IR IFOV, observing parts of neighbouring IR pixels and not fully observing the spatial extent of the matching IR IFOV.For the pixel close to the edge of swath, only two corresponding VIS pixels are available.Both fall within the IR IFOV but don't cover the full spatial extent.
To address these limitations, we have developed a new SLSTR preprocessor (McCarroll and Embury, 2022), which remaps the VIS and SWIR data provided on the image grids to the IR IFOV, taking into account the individual pixel locations and orphan pixels and avoiding the need to coarsen the VIS data with 2 × 2 averaging.The preprocessor first identifies a local neighbourhood including the N-nearest VIS or SWIR pixels to the IR pixel centre, where N is configurable.For all the results presented in this manuscript a value of N = 5 has been used for channels with A-stripe data only, and N = 10 for channels with A and B stripe data (where the ten nearest pixels are chosen independently of the detector).Searches for the N-nearest neighbours are made within a defined region of the image, bounded by a distance threshold, D. The location of each pixel centre is defined by the along-track and across-track coordinates.For each VIS or SWIR pixel in the defined region, the Cartesian distance to the centre of the IR pixel is calculated.Orphan pixels are included in the search and any duplicates are used only once.Where stripe A and B data are available (S4-S6), separate neighbourhoods can be calculated for each detector and then merged, so that the N-nearest neighbours can be identified from all available observations.The VIS or SWIR data are then mapped to the IR pixel by taking the arithmetic mean of the Nnearest neighbours.The preprocessor also supports calculation of the maximum, maximum-minimum and standard deviation of the N-nearest neighbours.
Subplots 3b and 3e show the improvement in co-registration of the VIS and IR data when using the new preprocessor.In this example we use only the A stripe and 5 nearest neighbours.This gives a much closer match compared to Fig. 2c between the spatial extent of the IR and VIS pixel IFOVs.At the swath edge, the VIS pixels cover most of the IR IFOV, whilst at the sub-satellite point, the VIS domain slightly exceeds the IR IFOV.Subplots 3c and 3f show the co-registration when using 10 nearest neighbours with both stripes A and B. The A stripe detector centres are shown in indigo and the B stripe detector centres in yellow.Using the B stripe where available doubles the number of observations within the IR pixel IFOV, where the A and B detector centres are closely aligned.The new preprocessor effectively calculates 'C stripe' data at 1 km resolution for channels S4-S6, combining the information from the A and B stripes and is better optimised to benefit any application using VIS data co- registered with IR channels.
The improved co-registration of the VIS and IR data can also be seen when applying the new preprocessor to the SLSTR image of the Norwegian coastline as shown in Fig. 2. When comparing the S7 brightness temperatures (Fig. 4a) to the S6 reflectance (4b) using the new preprocessor, we see a much closer match between the shape of the coastline in the data from the two channels.Fig. 4c shows the standard deviation of the S6 observations within the IR IFOV.This additional metric is also found to be useful in the cloud detection process, as described in section 4. The reflectance variability is elevated over land (as shown here) but also along reflectance gradients such as those seen at cloud edges over a dark ocean.

Bayesian cloud detection
Operational SLSTR SST products are generated in near real time by EUMETSAT.A Bayesian cloud detection algorithm is used to screen the data prior to retrieving the surface temperature.This algorithm has been widely described elsewhere (Merchant et al., 2005(Merchant et al., , 2019;;Bulgin et al., 2018) so we provide only a brief overview here, focusing on details specific to SLSTR processing.Bayes theorem as applied to cloud detection can be used to calculate the probability of clear-sky conditions as shown in Eq. (1).
The clear-sky probability P ( c|y o , x b ) is calculated given knowledge of the satellite observations (y o ) and the prior background state (x b ).We use the term 'cloud' here to describe all non-clear conditions where SST cannot be retrieved.P(c) and P(c) are the prior probabilities of clear-sky and cloud respectively.The prior probability of cloud is determined by the numerical weather prediction (NWP) clear-sky fraction, constrained to within the range 0.5-0.95 in order to avoid strongly preconditioning the cloud detection by the prior.The prior clear-sky probability is 1 − P(c).
The probability of the observations given the background state P ( y o |x b , c ) can be split into two terms, with a spectral component denoted by the subscript 's', and a textural component denoted by the subscript 't' (Merchant et al., 2005).
For clear-sky observations, P ( ) is simulated using the fastforward model RTTOV v11.3 (Hocking et al., 2015) constrained by NWP data from the European Centre for Medium Wave Forecasting (ECMWF) (Cox et al., 2021), provided with the SLSTR data.For cloudy observations, an empirical probability density function (PDF) is used to describe the prior background state.This is a practical choice given the timeliness of the data and the computational expense of simulating all possible cloud configurations.
The empirical PDFs are constructed using ten years of data from the Metop-A Advanced Very High Resolution Radiometer (AVHRR) spanning 2007-2016.PDFs were required prior to the launch of SLSTR to commence operational data production.AVHRR data were used in preference to ATSR due to the limited view angle ranges of ATSR.Cloudy-sky data for inclusion in the PDFs were bootstrapped using the EUMETSAT AVHRR cloud mask (EUMETSAT, 2011) and then iterated through one cycle of the Bayesian cloud detection (Merchant et al., 2019).
The textural component P ( ) is represented by an empirical PDF for both clear-sky and cloudy conditions.The texture is calculated as the standard deviation of the observation in a given channel over the pixel of interest and the eight surrounding pixels.These PDFs are derived using data from the Advanced Along-Track Scanning Radiometer (AATSR).
To use empirical PDFs from another sensor requires application of a shift to the SLSTR data prior to look-up to account for the differences in the spectral response functions of the two sensors (Bulgin et al., 2018).Essentially we make the SLSTR data 'look-like' the data used in the PDF, just prior to indexing the PDF.The spectral shifts are cubic functions, applied as a function of total column water vapour (TCWV) to the data at infrared wavelengths (Bulgin et al., 2018).Coefficients are calculated for atmospheric path lengths of 1 and 1.8 (Table 2) corresponding to viewing zenith angles between 0 and 56.25 • with linear interpolation applied to the coefficients between the two atmospheric path lengths.

Operational configuration
The current operational SLSTR daytime cloud detection uses three reflectance channels (S2, S3 and S5) and two infrared channels (S8 and S9) in the calculation of the spectral probability of the observations given the prior background state vector (P(y s o | x b , c) (Pearson et al., 2017).For cloudy observations, two empirical PDFs are used as described in Table 3.Over the open ocean and at sea-ice boundaries, it is typically found that adding channels in the visible part of the spectrum aids Bayesian cloud detection (Bulgin et al., 2015(Bulgin et al., , 2018)).In these regions, the reflectance wavelengths help to distinguish between brighter ice or cloud surfaces and the darker ocean.They are particularly useful where the temperature differences between surfaces are small such as at sea-ice  edges or in regions of low-lying fog.The trade-off for using these wavelengths occurs in coastal areas, where optically bright waters can be erroneously flagged as cloud.An example of this is presented in section 5.This arises when regions of turbid water or suspended matter in coastal regions cause enhanced reflectance in the S2 and S3 channels, which is not captured by the radiative transfer model.The textural probability P ( ) is defined using the standard deviation of the S8 channel over 3 × 3 pixels, centred on the pixel of interest (Fig. 5 (left)).The clear-sky PDF peaks at low local standard deviations (LSD), <0.2 K.For LSD values >0.2 K, cloudy conditions are more probable.This probability is particularly useful for identifying cloud edges and scattered cloud, where there is greater heterogeneity in the brightness temperature than over the open ocean (or uniform cloud decks).The trade-off when using this metric, is the sensitivity to thermal gradients in the ocean caused by SST fronts, which are typically masked by the Bayesian cloud detection (Bulgin et al., 2018).This metric may also over-flag scattered cloud and waters along the coast where strong thermal contrasts are found.An example of this is shown in section 5.

Use of optimised imagery and standard deviation
We propose here a new configuration for the cloud detection which is of benefit to cloud detection in coastal regions.Three primary changes are made with respect to the operational configuration: 1) We use the new preprocessor to match the reflectance channel observations to the infrared image grid.2) We use only the S5, S8 and S9 channels in the spectral component of the Bayesian calculation during the day and 3) we use the sub-pixel variability in the S3 channel as the texture metric during the daytime.
Changes 2) and 3) are designed to address the shortcomings in the current cloud detection configuration as identified in section 4.2.The first of these reduces the dependence of the cloud detection on the VIS channels where optically bright or turbid waters coastal waters are not well modelled.The second is expected to be beneficial to the cloud detection as clear-sky oceans are generally invariant with respect to reflectance on the kilometres scale of infrared pixels.The exceptions to this are sharp fronts separating water masses with contrasting turbidity, and near-specular conditions where there may be a sun-glint-related reflectance gradient modulated by short-scale changes in wind-driven surface roughness.The reduced spatial footprint of the textural probability should also improve edge detection of cloudy features.
We use the sub-pixel variability in S3 (s) in the cloud detection by generating a look-up table of the ratio, ϕ, between the probability of s in the cloudy and clear-sky cases.ϕ = P(s|c)P(s|c) (3) P(s|c) reflects the distribution of the sensor noise, which given a sample of N reflectances can be assumed to be drawn from a normal distribution over a uniform scene, where the standard deviation of that population (σ) is channel dependent.In principle, the reflectance uncertainty in the product should provide σ, but in the present operational data this is not available.Therefore we use SLSTR-A clear-sky extracts to calculate a typical fixed value of σ, which is 0.00065 for S3.Using the standard result that the probability density function of x = Ns 2 /σ 2 is described by the function χ 2 we derive P(s|c).
Using the parameter ρ = s/σ collapses the independent variations of s and σ onto a single parameter, suitable for use in a look-up table as shown in Fig. 5 (right).
The distribution of s for non-clear samples will include significant geophysical variability in addition to sensor noise, the distribution of which cannot be determined from first principles.Using non-clear samples of SLSTR-A data where the clear-sky probability is <0.1, for values of s much greater than the noise level, the distribution is well represented by an exponential distribution, where ν = 0.02273.
This fit was derived over a range beyond the influence of reflectance noise, so this exponential distribution is then convolved with the sensor noise to give the full distribution of P(s|c).This is stored as a numerical look up table, since the convolved distribution has no convenient analytical form.We introduce the parameter ϖ = ν/σ to express this distribution also in terms of ρ.
The result is shown in Fig. 5 (right).All values of ρ > 2 represent cloudy conditions.

Case study examples
Cloud detection performance in coastal zones is assessed qualitatively using 17 SLSTR 3-min granules.These are globally distributed, including the Yellow Sea, Red Sea, Gulf of Finland, Cape of Good Hope, Chile, Gulf of Venice, Portugal, the Mediterranean, Namibia and the Amazon delta.Two case studies are presented here in detail that illustrate the main differences between the operational and new cloud detection configurations (Fig. 6).The top panel shows the Bohai and Yellow Seas as observed on the 31st Jan 2020.The S2 channel reflectance (6a) shows optically bright waters over much of the Bohai Sea, around the entrance to the Yellow Sea and along the coast of Lianyungang and Yancheng (towards the bottom of the images).These areas are indicated by the red boxes in Fig. 6.Most of these areas are C.E. Bulgin et al. erroneously flagged as cloud when using the operational configuration of the cloud detection.Fig. 6b shows the clear-sky probability from the operational Bayesian calculation, to which a threshold of 0.9 would typically be applied to generate a binary cloud mask (Merchant et al., 2019).Therefore only regions that are dark blue over the ocean in these plots would be classified as clear-sky.Results for the new configuration are shown in 6c.Much of the optically bright waters previously flagged as cloud are correctly classified as clear-sky in this configuration (an C.E. Bulgin et al. increase from 46.2% to 93.1% in the clear-sky fraction inside the red boxes), primarily resulting from replacing the spectral S2 and S3 information (where optically bright waters are poorly modelled) with the S3 sub-pixel texture.The one area still falsely flagged as cloud is Liadong Bay (top, centre of the image).The waters in this region are cold, and the sub-pixel S3 variability is elevated, suggesting high turbidity.
The second example (Fig. 6, second row) is around the South Korean coast on the 21st January 2020 where there are a number of thermal SST fronts to the south and west of the mainland as shown in the S8 infrared image (6d, magenta).In the operational configuration (6e) the cloud detection is sensitive to many of these fronts along both the south-west and south coasts as they raise the local variability in the S8 brightness temperature used as the textural metric.In the new configuration (6f), replacing the S8 textural metric with the S3 sub-pixel variability decreases the sensitivity to thermal fronts.To the north-west of the image there is an area of cloud, in which a number of contrails are visible in the S8 thermal imagery (Fig. 6, cyan box in the second row, and an expanded view in the bottom row).These features are better defined in the new cloud detection configuration due to the reduced spatial footprint of the S3 textural metric.This metric also increases the sensitivity of the cloud detection to scattered cloud fields often located at cloud edges (orange boxes).

Coastal zone performance metrics
Quantitative assessment of cloud detection performance requires a 'truth' against which the mask can be compared.This 'truth' is often constructed by expert inspection of satellite imagery or semi-automated labelling of pixel clusters (Skakun et al., 2022;Bulgin et al., 2018), but these datasets are typically limited due to the time-consuming nature of manual classification (Bulgin et al., 2014).In this study we make use of two methods to validate coastal cloud detection performance.The first uses the S7 (3.7 μm) brightness temperature to identify clouds in cases where they are markedly warmer than the underlying ocean surface.This arises when the cloud top properties and sun-satellite geometry are such that significant solar irradiance is scattered towards the sensor by the cloud tops (Hunt, 1973).The S7 channel is particularly sensitive to fields of broken cloud, where cloud edges enhance reflection (Coakley and Davies, 1986) and to ice particle size in cirrus cloud (smaller particles absorb less solar radiation increasing reflectance) (Baran et al., 2003).In these cases, the scattering increases the observed brightness temperature in the S7 channel.When evaluating both the S7 and S8 brightness temperatures a scene-dependent threshold can be readily found by inspection that enables cloud discrimination at the full spatial resolution of the thermal image.Warm S7 clouds occur in 47% of the 17 worldwide granules, as identified by visual inspection (including the Yellow Sea, Chile, Gulf of Finland, Portugal, Namibia, Red Sea and the Amazon Delta).From these granules, 12 areas are extracted by inspection, in which clouds are warm in the S7 channel and masked in order to evaluate the Bayesian cloud detection (Fig. 7).
An example of this scattering effect and the threshold-based cloud 'truth' is shown in Fig. 8.We extract a halo of 20 pixels (corresponding to a distance of 20 km at nadir) around the coastline for each of these segments (Fig. 8a).From these pixels we generate a 2D density plot of the observations in the S7 and S8 channels (Fig. 8c) with a bin size of 1 K 2 .Clear-sky pixels lie close to the 1:1 line, with a scene-dependent offset in the S8 channel due to the greater sensitivity of the brightness temperature to atmospheric water vapour at this wavelength (Sobrino et al., 2003).To set the threshold for classification of cloudy pixels we identify the bin with the highest count and use the distance of this bin from the 1:1 line to determine the offset (c) in the threshold along the yaxis (the S7 channel).
To this offset we add an additional 2 K representing the uncertainty in the observations from sensor noise and water vapour variability and assume a slope of one for the threshold.This is illustrated in Fig. 8c, where all pixels in the top-left of the plot above the dotted line would be classified as cloud and those in the bottom right as clear-sky.The resultant cloud mask is shown by the blue shading in Fig. 8b.Using this form of validation provides a high-level of objectivity and benefits from the inclusion of pixels in the clear-sky to cloud transition zones, which are commonly the hardest to classify.In Table 4 we provide a full list of the SLSTR-A granules used in this validation along with the x and y limits of the extracts and the scene-dependent threshold offset, c.
The second validation method is used in areas where clouds are not uniformly warm in the S7 channel.In these regions, clusters of clear-sky pixels and clusters of cloudy-pixels are identified by expert inspection,

Table 4
SLSTR-A segments, extraction limits and offset used for coastal cloud detection performance analysis.The unique beginning of the L1 filename can be constructed using the prefix 'S3A_SL_1_RBT____' followed by the date and time in the following format: 20200121 T012840.where the classification is unambiguous (as shown in Fig. 7).This validation is included to provide confidence that the algorithm works across different atmospheric regimes but doesn't include the regions of transition between clear-sky and cloud (where classification is arguably more difficult), by its methodology.We calculate four performance metrics for the operational and new configurations of the cloud mask (Bengtsson and Hodges, 2005;Mackie et al., 2010): the percentage of perfect classification (PP), hit rate (HR), false alarm rate (FAR) and true skill score (TSS).The percentage of perfect classification is defined as the percentage of pixels correctly classified as 'cloud' or 'clearsky over water' with reference to the 'truth'.The hit rate is the percentage of cloudy pixels correctly identified as cloud and the false alarm rate is the number of clear-sky pixels falsely flagged as cloud.The true skill score is defined as the hit rate minus the false alarm rate.
The performance metrics presented in Table 5 are calculated for the two validation methods.For the regions where S7 clouds are warm (referred to as scene extracts in Table 5) there are a total of 567,439 coastal pixels, of which 161,625 are cloudy and 405,814 are clear-sky.The new configuration increases the hit rate by 1.7% compared to the operational configuration, whilst simultaneously reducing the false alarm rate by 31.6%.This leads to a net gain in the true skill score of 33.3% for the new configuration (74.2% compared to 40.9% for the operational configuration).The improvement in the percentage of perfect classification is considerable (84.8% for the new configuration compared to 61.7% for the operational configuration).The new configuration significantly improves identification of clear-sky pixels in coastal regions by reducing the false alarm rate, whilst also improving the hit rate, both of which are advantageous in the retrieval of surface properties such as sea surface temperature.
For the pixel cluster analysis, a total of 39,518 pixels were evaluated, of which 24,070 are cloud and 15,448 clear-sky.For all metrics, both algorithms perform better for the cluster analysis than their corresponding performance in the scene extracts, likely because the clusters only include pixels that are unambiguous in their classification.A similar increase in performance is seen for the new configuration in comparison to the operational configuration in the pixel cluster analysis.The hit rate increases by 3.3% to 94.3%, whilst the false alarm rate falls by 23.8%.The net gain in the true skill score between the new and operational configuration is 27.1% and the corresponding increase in the percentage of perfect classification is 11.3%.

Impact on global validation statistics
Both case study analysis and calculation of performance metrics have demonstrated the positive impact of the new cloud detection configuration in coastal zones in the detection of scattered cloud, correct classification of thermal fronts and reduced false-flagging of optically bright waters.The primary application of the Bayesian cloud detection algorithm operationally is in the identification of clear-sky over ocean pixels for SST retrieval.Having demonstrated a clear improvement in cloud detection in coastal regions, global, open-ocean validation is required to ensure that the proposed modifications are also beneficial (or at least not detrimental) over open water and across the globe.Validating SST retrievals from satellite data against drifting buoys is a good test of cloud screening performance.Any reduction in cloud detection skill will introduce cold biases and worsen the validation statistics, whilst an increase in false-flagging of cloud would significantly reduce the number of clear-sky matches.Table 6 shows global validation statistics for a nadir-only, two-channel (S8 and S9, 11 and 12 μm), coefficient-based retrieval (N2) (Embury et al., 2012a) against global drifting buoy data, packaged and distributed by the Copernicus Marine Environment Monitoring Service (CMEMS) for all daytime SLSTR-A observations in 2021.Matches are constrained by the spatial and temporal differences between the satellite and buoy measurements with a maximum separation 1 km and two hours.
Comparisons show the satellite minus buoy SST's with skin-to-depth comparisons of the raw data and skin-to-skin comparisons where the Fairall and Kantha-Clayson (FKC) models (Kantha and Clayson, 1994;Fairall et al., 1996), driven by ERA5 fluxes (C3S, 2021) are used to adjust buoy temperatures (taken at depth) to skin temperatures.This process occurs in two stages: the buoy temperature is first converted to a skin temperature and then this skin temperature is adjusted for the time difference between the satellite and buoy observations.The Kantha-Clayson model accounts for diurnal stratification and time-related changes to the water temperature at the depth of the buoy, whilst the Fairall model accounts for the physics of the cool skin layer (Embury et al., 2012b).These models are used together for each step.
The non-robust mean and standard deviations are identical for the operational and new configurations with a reduction of 6% in the number of observations classified as clear-sky when using the new configuration.There is a small reduction in the total number of observations in the skin-to-skin case compared to skin-to-depth (for both nonrobust and robust statistics), that result from instances where the FKC model predictions are not available.The robust statistics show an improvement in the robust standard deviation of the SST differences (0.30 for the operational configuration compared to 0.29 for the new configuration), where the new configuration removes more cloudy observations.This indicates that the new configuration does not degrade cloud detection performance over the open ocean.Fig. 9 compares the distribution of the SST differences for the operational and new configurations.There is a cold tail to the distribution (represented by the median difference of − 0.04 K in the skin-to-skin comparisons) that represents cloudy observations mis-classified as clear-sky and we see a small reduction in this when using the new configuration.Although this improvement is present also in the non-robust statistics, it is less obvious in the comparison statistics as outliers are typically driven by other aspects of the data (poor quality in-situ, the FKC model and aerosol).

Discussion and conclusion
The three main improvements expected from using the new cloud detection configuration during the day are 1) a reduction in the falseflagging of thermal fronts, 2) a reduction in the false-flagging of optically bright waters and 3) an improvement in the screening of scattered cloud.The first two improvements will not significantly impact the global statistics when comparing against drifting buoys as the main gains in clear-sky information will be close to the coast, where buoy

Table 5
Coastal zone performance metrics for the operational and new cloud detection configurations.Metrics are the percentage of perfect classification (PP), hit rate (HR), false alarm rate (FAR) and the true skill score (TSS) (Bengtsson and Hodges, 2005;Mackie et al., 2010).coverage is limited.The third improvement on the cloud detection will have the biggest impact on the global statistics by improving the identification of small scattered clouds, which are harder to differentiate using the infrared texture metric with the larger spatial footprint.It is likely that these features contribute to the majority of the reduction in clear-sky matches between the operational and new configurations.
In this paper we have demonstrated the benefits of the daytime cloud screening developments to data coverage in coastal zones, which can be achieved without detrimental impact to global validation statistics.This increased coverage in coastal regions is beneficial to SST operations and climate data records, the importance of which is increasingly recognised because a large fraction of the world's population live close to the coast.The higher resolution reflectance channel data from SLSTR is a powerful tool for detecting scattered cloud but is also sensitive to highly turbid waters as seen in the Bohai Sea case study.Further work is required to successfully distinguish clear-sky observations under these conditions.
More broadly, the SLSTR preprocessor facilitates re-mapping the VIS and SWIR data to the IR IFOV and would benefit any application in which the joint use of reflectance and infrared data is required.It also provides additional metrics including the maximum, the range and the sub-pixel variability of the VIS or SWIR data corresponding to any given IR pixel, which may have benefits for other applications, in addition to cloud detection.The code has been made freely and publicly available for download (McCarroll and Embury, 2022).

Fig. 1 .
Fig. 1.SLSTR viewing geometry showing the satellite zenith angle across-swath for the nadir and oblique views (a), and the IR pixel field of view relative to the A and B detector pixel sensors for the nadir view at the sub-satellite point (b), the oblique view centre-swath (c) and the nadir view edge of swath (d).

Fig. 2 .
Fig. 2. SLSTR S7 brightness temperatures (a) and S6 radiance (b) along part of the Norwegian coastline at 1 km and 500 m resolution respectively.(c) shows S6 data coarsened to the resolution of the S7 data using a simple 2 × 2 averaging of pixels.Data are extracted from the following L1 file: S3B_SL_1_RBT____20200422T110405_20200422T110705_20200423T161953_0179_038_094_1980_LN2_O_NT_004.SEN3.

Fig. 4 .
Fig. 4. (a) SLSTR S7 brightness temperatures along part of the Norwegian coastline at 1 km resolution.(b) shows S6 radiance remapped to the resolution of the infrared observations using the new preprocessor, with stripe A data only and five nearest neighbours.(c) shows the standard deviation of the S6 radiance over the footprint of each corresponding infrared observation.Data are extracted from the following L1 file: S3B_SL_1_RBT____20200422T110405_20200422T110705_ 20200423T161953_0179_038_094_1980_LN2_O_NT_004.SEN3.

Fig. 6 .
Fig. 6.Cloud detection for operational (centre) and new (right) configurations for optically bright waters (top) and thermal SST fronts (second row).a) S2 (0.6 μm) reflectance over the Bohai Sea (31/01/2020) and corresponding operational (b) and new (c) clear-sky probability.d) S7 (3.7 μm) brightness temperature around South Korea (21/01/2020) and corresponding operational (e) and new (f) clear-sky probability.Highlighted regions show optically bright waters (red), thermal SST fronts (magenta), contrails (cyan) and cloud edges (orange).The bottom row (panels g, h, i) show the region of contrails marked by the cyan box in the second row in more detail: (g) S7 brightness temperature, (h) operational clear-sky probability and (i) new clear-sky probability.The first example scene is processed from the following L1 file: S3A_SL_1_RBT____20200131T0209_20200131T02124_20200131T04080_0179_ 054_217_2340_LN2_O_NR_004.SEN3 and the second from: S3A_SL_1_RBT____20200121T012840_ 20200121T013140_20200121T032605_0179_054_074_2340_LN2_O_NR_004.SEN3.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. SLSTR scenes used for coastal validation of the new cloud detection algorithm.Figures show the full SLSTR segment: coastal zones display S7 brightness temperatures, open ocean is blue and land is green.Cyan boxes show the regions selected for validation using warm S7 temperatures to identify clouds.Pixel clusters used as an alternative validation metric under regimes where the S7 temperature of clouds is not uniformly warmer than the underlying ocean are shown for clear pixels (yellow) and cloudy pixels (red).The details of the SLSTR filenames for these scenes and extract locations are provided in Table 4.The rightmost cyan box in subplot 'b)' is the extract detailed in Fig. 8. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .
Fig. 8. Threshold-based cloud masking of test-scenes with warm S7 clouds.a) S7 brightness temperatures for a halo of 20 pixels around the coastline.Land is shown in green, open ocean excluded from the analysis in white.b) Resultant threshold-based cloud mask overlaid in blue on the S7 brightness temperature field.c) 2D density plot using the S8 and S7 brightness temperatures used to determine the threshold for identifying cloudy pixels.The example shown here is extract 6 from Table 4. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Probability distribution function (PDF) for N2 satellite minus drifter SST for clear-sky scenes, defined using a threshold of 90% on the clear-sky probability.The PDF for the operational configuration is shown in orange and the new configuration in blue.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
SLSTR channel characteristics including number, wavelength, spatial resolution and detectors.

Table 2
Spectral shift coefficients applied as a function of TCWV for SLSTR infrared channels.
C.E.Bulgin et al.

Table 3
Empirical PDF structure used to calculate the spectral probability of cloud.
Fig. 5. Daytime S8 local standard deviation (LSD) PDF (left) and S3 Rho PDF (right).Clear-sky conditions are shown in blue, cloud in orange.(Forinterpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)C.E.Bulgin et al.

Table 6
Validation statistics for N2 two-channel SST retrievals against drifting buoys for all SLSTR-A data in 2021.Clear-sky data are identified using the operational and new configurations of the Bayesian cloud detection.Non-robust statistics show number of matches, mean (standard deviation) and robust statistics show number of matches, median (robust standard deviation).