High-Resolution Sea Surface Temperature Retrieval from Landsat 8 OLI / TIRS Data at Coastal Regions

: High-resolution sea surface temperature (SST) images are essential to study the highly variable small-scale oceanic phenomena in a coastal region. Most previous SST algorithms are focused on the low or medium resolution SST from the near polar orbiting or geostationary satellites. The Landsat 8 Operational Land Imager and Thermal Infrared Sensor (OLI / TIRS) makes it possible to obtain high-resolution SST images of coastal regions. This study performed a matchup procedure between 276 Landsat 8 images and in-situ temperature measurements of buoys o ﬀ the coast of the Korean Peninsula from April 2013 to August 2017. Using the matchup database, we investigated SST errors for each formulation of the Multi-Channel SST (MCSST) and the Non-Linear SST (NLSST) by considering the satellite zenith angle (SZA) and the ﬁrst-guess SST. The retrieved SST equations showed a root-mean-square error (RMSE) from 0.59 to 0.72 ◦ C. The smallest errors were found for the NLSST equation that considers the SZA and uses the ﬁrst-guess SST, compared with the MCSST equations. The SST errors showed characteristic dependences on the atmospheric water vapor, the SZA, and the wind speed. In spite of the narrow swath width of the Landsat 8, the e ﬀ ect of the SZA on the errors was estimated to be signiﬁcant and considerable for all the formations. Although the coe ﬃ cients were calculated in the coastal regions around the Korean Peninsula, these coe ﬃ cients are expected to be feasible for SST retrieval applied to any other parts of the global ocean. This study also addressed the need for high-resolution coastal SST, by emphasizing the usefulness of the high-resolution Landsat 8 OLI / TIRS data for monitoring the small-scale oceanic phenomena in coastal regions.


Introduction
Sea surface temperature (SST) plays an important role in monitoring and understanding a variety of oceanic and atmospheric phenomena in space and time. In addition to the high-quality in-situ temperature measurements, including the data obtained by surface drifters, moored buoys, or ship of opportunity measurements, accurate satellite-observed SSTs are of importance for diverse purposes and have been operationally used with satisfaction [1,2]. The satellite-observed SST databases have been produced and extensively utilized for investigating SST fronts [3,4], mesoscale eddy dynamics [5], ocean surface current [6], air-sea interaction [7], and their effects on climate change [8]. They are also used as input data for atmospheric and oceanic circulation models as well as for data assimilation purposes [9,10].
Using brightness temperatures (BT) of infrared sensors on satellites, the SSTs have been retrieved based on multi-channel algorithms since the 1970s [11,12]. Multi-Channel SST (MCSST) [13][14][15][16] and Non-Linear SST (NLSST) [17][18][19][20] have been the most commonly used infrared-based SST retrieval algorithms as operational algorithms for the National Oceanic and Atmospheric Administration tidal currents in the Yellow Sea, however, the bathymetry of the coastal regions of the EJS is comparatively deep, amounting to a few thousand meters (Figure 1b). For the SST retrieval procedure, a wide range of SST variations is important for the stable use of the SST retrieval equations. The seasonal variations of SSTs are dominant with a wide range of SSTs between winter and summer as shown in the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) SST climatology for the period of 2002 to 2015 (Figure 1c,d). In addition to the range of SST variations, the accuracy of infrared channel-based SST is extensively affected by changes in the marine environment, such as bathymetry, land mass, tidal flat, tidal current, river discharge, and so on [56,57]. Therefore, this study region is suitable for the calculation of infrared SST coefficients and validation of the equations because of these diverse marine and atmospheric characteristics. Nevertheless, none of the previous studies have focused on deriving the SST coefficients or the Landsat 8 OLI/TIRS data for this region. Furthermore, SST coefficients for the Landsat 8 data have not been provided for coastal regions even in other parts of the global ocean. This might be due to insufficient matchups between the Landsat 8 OLI/TIRS data and in-situ temperature measurements in addition to cloud conditions and other local atmospheric conditions. (a) Location of the study area with contours of the water depth (m) in the seas around the Northeast Asia including China, Japan, Korea, and Russia, which red box indicates the study area, (b) a schematic current map with cold (blue) and warm (red) currents [58], and monthly mean of sea surface temperature (SST) climatology (°C) estimated from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) SST database in the study area from 2002 to 2015 in (c) winter (January) and (d) summer (July). (a) Location of the study area with contours of the water depth (m) in the seas around the Northeast Asia including China, Japan, Korea, and Russia, which red box indicates the study area, (b) a schematic current map with cold (blue) and warm (red) currents [58], and monthly mean of sea surface temperature (SST) climatology ( • C) estimated from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) SST database in the study area from 2002 to 2015 in (c) winter (January) and (d) summer (July).

Satellite Data
Landsat 8 has five channels in the visible spectrum (Band 1, Band 2, Band 3, Band 4, and Band 8) and 6 channels in the infrared spectrum (Band 5, Band 6, Band 7, Band 9, Band 10, and Band 11) [44]. The spatial resolution of Landsat 8 OLI/TIRS is 15 m × 15 m in Band 8 (Panchromatic band), 100 m × 100 m in Band 10 and Band 11, and 30 m × 30 m in other bands. The repeat cycle of the Landsat 8 is 16 days and the satellite observes the Korean Peninsula from 02 Universal Time Coordinated (UTC) to 03 UTC. In this study, we selected all of the Landsat 8 OLI/TIRS images over the coastal region of the Korean Peninsula as indicated by the rectangular edges of the images in Figure 2. In total, 276 images from the Landsat 8 OLI/TIRS were processed for the period from 1 April 2013 to 1 August 2017. All Landsat 8 OLI/TIRS data used in this study were the L1TP data obtained from NASA and the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/).

Satellite Data
Landsat 8 has five channels in the visible spectrum (Band 1, Band 2, Band 3, Band 4, and Band 8) and 6 channels in the infrared spectrum (Band 5, Band 6, Band 7, Band 9, Band 10, and Band 11) [44]. The spatial resolution of Landsat 8 OLI/TIRS is 15 m × 15 m in Band 8 (Panchromatic band), 100 m × 100 m in Band 10 and Band 11, and 30 m × 30 m in other bands. The repeat cycle of the Landsat 8 is 16 days and the satellite observes the Korean Peninsula from 02 Universal Time Coordinated (UTC) to 03 UTC. In this study, we selected all of the Landsat 8 OLI/TIRS images over the coastal region of the Korean Peninsula as indicated by the rectangular edges of the images in Figure 2. In total, 276 images from the Landsat 8 OLI/TIRS were processed for the period from 1 April 2013 to 1 August 2017. All Landsat 8 OLI/TIRS data used in this study were the L1TP data obtained from NASA and the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/). In order to address the advantage of high-resolution SST from Landsat 8 OLI/TIRS, we compared the high-resolution SST with the SST data produced from other satellites. We collected the Global Change Observation Mission 1st-Water (GCOM-W1) Advanced Microwave Scanning Radiometer 2 (AMSR-2), the Himawari-8 the Advanced Himawari Imager (AHI), and the NOAA/AVHRR collocated with the Landsat 8 OLI/TIRS data. GCOM-W1/AMSR-2 data used in this study were L2 data obtained from the Remote Sensing System (http://www.remss.com/missions/amsr/), Himawari-8 AHI data used in this study were L2 data obtained from the National Meteorological Satellite Center (NMSC), KMA, and NOAA/AVHRR data used in this study were L2 data obtained from the Research Institute of Oceanography (RIO), Seoul National University (SNU).

In-Situ Measurements
To derive the coefficients of SST algorithms for SST retrieval from the Landsat 8 OLI/TIRS data, we collected in-situ measurements. Since the Landsat observations have been concentrated at the coastal regions, the data from surface drifters were not sufficient to be used for the estimation and validation of SST coefficients because of the few observations. Korean Meteorological Administration (KMA) has been operating the marine meteorological buoy stations to monitor the coastal In order to address the advantage of high-resolution SST from Landsat 8 OLI/TIRS, we compared the high-resolution SST with the SST data produced from other satellites. We collected the Global Change Observation Mission 1st-Water (GCOM-W1) Advanced Microwave Scanning Radiometer 2 (AMSR-2), the Himawari-8 the Advanced Himawari Imager (AHI), and the NOAA/AVHRR collocated with the Landsat 8 OLI/TIRS data. GCOM-W1/AMSR-2 data used in this study were L2 data obtained from the Remote Sensing System (http://www.remss.com/missions/amsr/), Himawari-8 AHI data used in this study were L2 data obtained from the National Meteorological Satellite Center (NMSC), KMA, and NOAA/AVHRR data used in this study were L2 data obtained from the Research Institute of Oceanography (RIO), Seoul National University (SNU).

In-Situ Measurements
To derive the coefficients of SST algorithms for SST retrieval from the Landsat 8 OLI/TIRS data, we collected in-situ measurements. Since the Landsat observations have been concentrated at the coastal regions, the data from surface drifters were not sufficient to be used for the estimation and validation of SST coefficients because of the few observations. Korean Meteorological Administration (KMA) Remote Sens. 2019, 11, 2687 5 of 21 has been operating the marine meteorological buoy stations to monitor the coastal environments around the Korean peninsula as shown in Figure 2. In this study, we use data from 17 KMA marine meteorological buoys, covering the Landsat observational areas, around the coastal regions of the Yellow Sea, the EJS, and in the southern coastal region (https://kma.go.kr/). These KMA buoys observe SST at different depths from 0.1 m to 0.4 m and wind speeds at heights from 3.6 to 4.0 m every hour ( Table 1). The equipped thermistors of the buoys have a precision of 0.1 • C and an accuracy of ±0.1 • C for a temperature range from −2 to 40 • C (https://kma.go.kr/). To distinguish the locations of the buoys, we set the symbols of the buoys following the first alphabet of the names of the seas; Y1 to Y6 are buoys in the Yellow Sea, S1 to S6 are in the southern region, and E1 to E5 are in the EJS.

Daily SST Data
For estimating the coefficients using the NLSST formulation, the first-guess SST is needed and can be obtained from the estimated MCSST or climatological SST [18,59]. It is also possible to use other SST products from diverse SST databases. Considering the high-resolution Landsat data and complicated coastal environments for calculating high-resolution SST of less than 100 m, if possible, a high-resolution SST field is preferable as the first-guess SST. Therefore, this study used the OSTIA data provided by the Met Office with 6 km spatial resolution and the Multi-scale Ultra-High Resolution SST (MURSST) data provided by the Group for High-Resolution SST (GHRSST) with 1 km spatial resolution, as the first-guess SST data [60,61]. These SST products from OSTIA and MURSST provide optimally interpolated daily SST, using infrared sensor-based satellites, passive microwave-based satellite data, and in-situ measurements [61,62].

Removal of Cloud-Contaminated Pixels
Removing pixels contaminated by clouds is important for SST retrieval using infrared sensor-based satellite data. Cloud mask algorithms applicable to the Landsat 8 OLI/TIRS data have been developed by various studies, including the Automated Cloud Cover Assessment (ACCA), the Web Enabled Landsat Data (WELD) Cloud Cover Assessment (CCA) algorithm, the Function of Mask (Fmask) algorithm, and so on [63][64][65]. In this study, we applied the modified Moderate Resolution Imaging Spectroradiometer (MODIS) Cloud mask algorithm to the Landsat 8 OLI/TIRS data to remove the cloud pixels and the cloud contaminated pixels, prior to the SST retrieval [66,67]. Several tests of this Remote Sens. 2019, 11, 2687 6 of 21 algorithm have shown that it effectively detects snow and clouds pixels in the Landsat 8 OLI/TIRS images [68]. The Band 3, Band 1, Band 2, Band 6, Band 26, Band 31, and Band 32 of MODIS are substituted respectively with Band 2, Band 4, Band 5, Band 6, Band 9, Band 10, and Band 11 of the Landsat 8 OLI/TIRS. The cloud mask algorithm classifies the pixels of the Landsat 8 OLI/TIRS images into five types (non-vegetated land, vegetated land, snow/ice, water bodies, and clouds). However, in this study we masked the land pixels using the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data (https://earthexplorer.usgs.gov/) and applied the cloud mask algorithm to the pixels for classifying them into snow, ice, water bodies, or clouds ( Figure 3).
Removing pixels contaminated by clouds is important for SST retrieval using infrared sensorbased satellite data. Cloud mask algorithms applicable to the Landsat 8 OLI/TIRS data have been developed by various studies, including the Automated Cloud Cover Assessment (ACCA), the Web Enabled Landsat Data (WELD) Cloud Cover Assessment (CCA) algorithm, the Function of Mask (Fmask) algorithm, and so on [63][64][65]. In this study, we applied the modified Moderate Resolution Imaging Spectroradiometer (MODIS) Cloud mask algorithm to the Landsat 8 OLI/TIRS data to remove the cloud pixels and the cloud contaminated pixels, prior to the SST retrieval [66,67]. Several tests of this algorithm have shown that it effectively detects snow and clouds pixels in the Landsat 8 OLI/TIRS images [68]. The Band 3, Band 1, Band 2, Band 6, Band 26, Band 31, and Band 32 of MODIS are substituted respectively with Band 2, Band 4, Band 5, Band 6, Band 9, Band 10, and Band 11 of the Landsat 8 OLI/TIRS. The cloud mask algorithm classifies the pixels of the Landsat 8 OLI/TIRS images into five types (non-vegetated land, vegetated land, snow/ice, water bodies, and clouds). However, in this study we masked the land pixels using the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data (https://earthexplorer.usgs.gov/) and applied the cloud mask algorithm to the pixels for classifying them into snow, ice, water bodies, or clouds ( Figure 3).  Figure 3 shows a series of cloud removal procedures in which the algorithm was developed based on the spectral characteristics of water, snow/ice, and cloud in the visible and shortwave infrared spectrum [69]. The split window test and the cirrus test were added in the algorithm, and the threshold test used in the split window test was calculated by radiative transfer simulations [67]. The threshold of BT difference (BTD) at the 11 μm channel was used to remove large BTD values corresponding to the clouds. In addition, cirrus pixels were removed by applying a threshold of Band 9, which was similar to the cirrus band for the Landsat 8 OLI/TIRS [70,71].  Figure 3 shows a series of cloud removal procedures in which the algorithm was developed based on the spectral characteristics of water, snow/ice, and cloud in the visible and shortwave infrared spectrum [69]. The split window test and the cirrus test were added in the algorithm, and the threshold test used in the split window test was calculated by radiative transfer simulations [67]. The threshold of BT difference (BTD) at the 11 µm channel was used to remove large BTD values corresponding to the clouds. In addition, cirrus pixels were removed by applying a threshold of Band 9, which was similar to the cirrus band for the Landsat 8 OLI/TIRS [70,71].

SST Retrieval Algorithms
For SST retrieval from satellite data, MCSST and NLSST algorithms have been most commonly used, for the 11 and 12 µm channels. These MCSST and NLSST formulations have been used in various Remote Sens. 2019, 11, 2687 7 of 21 formula with respect to satellite data to various purposes for their uses [72][73][74][75]. In case of Landsat data, a simple formulation as the first equations of the MCSST (MCSST1) has been used. However, there is a high possibility that other formulations, containing a term of satellite zenith angle (SZA), would be much more appropriate than the simple MCSST formulation.
The swath of Landsat is relatively narrow with a SZA of up to 7.5 degrees on either side of the satellite orbit path as compared with that of the NOAA/AVHRR or other near-polar orbiting satellites [76][77][78][79]. Most of the previous studies have not considered the effect of SZAs in the estimation of SST errors for the Landsat data. Even studies on land surface temperature retrieval and SST retrieval have neglected SZAs under the assumption of nadir observation for all the pixels [50][51][52][53][54][55]80]. In this study, in order to analyze the effect of the SZA on the SST accuracy, we classified into the two cases, including considering and neglecting the SZA-related terms both in the MCSST and the NLSST formulations in Table 2. Therefore, we proposed eight SST equations, two cases for the MCSST and four cases for the NLSST with respect to the input data of first-guess SST and consideration of the SZA. For each formulation, we derived the SST coefficients and assessed the accuracy of each equation. Table 2. Formula based on MCSST and NLSST algorithms (T 11 and T 12 : Brightness temperature at 11 µm and 12 µm (in Celsius), T sfc1 , T sfc2 , and T sfc3 : MCSST, OSTIA SST, and MURSST as first-guess SST (in Celsius), θ: satellite zenith angle, a 1 , a 2 , a 3 , and a 4 : regression coefficients).

Algorithm
Symbol Equation

Calculation of Satellite Zenith Angle
The swath width of the Landsat 8 OLI/TIRS satellite is relatively small, hence, many studies usually assume that the viewing zenith angle α as the SZA, without considering the curvature of the Earth. In this case, the SZA is up to about 7.5 degrees, but it is essential to calculate the precise SZA for accurate SST retrieval. Therefore, in this study, in order to obtain the accurate SZA, we calculated it with by considering the curvature of the Earth [81]. Figure 4 shows the schematic diagram of the calculation of SZA of the Landsat 8 satellite. In the schematic diagram, the target pixel T for SZA retrieval is calculated using the distance d from the sub-satellite point S. Assuming the radius of the Earth R is 6371 km, the center angle of the Earth β can be calculated based on the law of the arc as follows: If the altitude h of the Landsat 8 OLI/TIRS is assumed to be 705 km, it is possible to calculate the distance b from the satellite to the target pixel according to the second law of cosine: Finally, we calculate γ based on the law of sine and the SZA, defined as 180 • -γ: Remote Sens. 2019, 11, 2687 8 of 21 In this way, the SZA value of each pixel was calculated for the ranges from satellite orbit path with zero degree to the end points on the left and the right side with the angles of about ±8.4 degrees.
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 21 In this way, the SZA value of each pixel was calculated for the ranges from satellite orbit path with zero degree to the end points on the left and the right side with the angles of about ±8.4 degrees.

Quality Control of In-Situ Measurements
The buoys, used in this study for calculation of the SST coefficients and their validation, have been providing marine environmental observations for a long period ( Table 1). The in-situ measurements from the KMA marine meteorological buoy are distributed primarily after the QC process, however, some anomalous observations are produced unexpectedly as well. Therefore, a QC process is necessary to accurately calculate and verify the SST data [82]. For this, a series of QC methods developed by the KMA were applied to the in-situ data. This method consists of sequential tests to remove low-quality data and outliers for each individual buoy data. The buoy data for QC should be composed of at least ten measurements within a day and the difference between the maximum and minimum SST within a day should be greater than 0 °C and less than 4 °C for the monitoring of daily SST fluctuations. For the time-continuity test, the difference of SST from the daily mean should be less than three standard deviations of the temperatures within 1 and 4 days. The standard deviation of SST within 4 days should be less than 2 °C. Figure 5 shows an example of SST QC results before and after applying the QC process to the buoy measurements. Prior to the QC process, the time-series of the buoy temperatures revealed abnormal lower peaks and higher peaks by deviating from general SST fluctuations (Figure 5a,b). After applying a QC process similar to a despiking procedure, the unusual peaks were completely eliminated as shown in Figure 5c-f.

Quality Control of In-Situ Measurements
The buoys, used in this study for calculation of the SST coefficients and their validation, have been providing marine environmental observations for a long period ( Table 1). The in-situ measurements from the KMA marine meteorological buoy are distributed primarily after the QC process, however, some anomalous observations are produced unexpectedly as well. Therefore, a QC process is necessary to accurately calculate and verify the SST data [82]. For this, a series of QC methods developed by the KMA were applied to the in-situ data. This method consists of sequential tests to remove low-quality data and outliers for each individual buoy data. The buoy data for QC should be composed of at least ten measurements within a day and the difference between the maximum and minimum SST within a day should be greater than 0 • C and less than 4 • C for the monitoring of daily SST fluctuations. For the time-continuity test, the difference of SST from the daily mean should be less than three standard deviations of the temperatures within 1 and 4 days. The standard deviation of SST within 4 days should be less than 2 • C. Figure 5 shows an example of SST QC results before and after applying the QC process to the buoy measurements. Prior to the QC process, the time-series of the buoy temperatures revealed abnormal lower peaks and higher peaks by deviating from general SST fluctuations (Figure 5a,b). After applying a QC process similar to a de-spiking procedure, the unusual peaks were completely eliminated as shown in Figure 5c Figure 6a shows the spatial distribution of the number of matchup data at the buoy locations off the coast of the Korean peninsula from April 1, 2013 to August 31, 2017. The total number of the matchup data between the Landsat 8 OLI/TIRS images and the KMA ocean buoy data amounted to 320 and the largest number of matchups was 43, for the buoy Y1. The smallest number of the matchups was 6, for the buoys Y2 and S3 (Figure 6a). During the study period over five years, there were no matchups between the Landsat 8 data and the buoy data in some months (June 2013, June 2014, August 2014, January 2015, and June 2016) due to the mask associated with clouds, as shown in Figure 6b. The month-year plot of the matchup numbers in Figure 6b exhibits the maximum number of the matchups of 14 on April 2017. The in-situ temperatures of the matchup database ranged from 3.2 to 30.2 °C and the maximum frequency of occurrence of in-situ temperatures was in the range of 14.5 to 15.5 °C (Figure 6c). The KMA ocean buoys measure not only the sea water temperature at a depth of less than 0.4 m but also the wind speed. The in-situ wind speed of the matchup database ranged from 0.1 to 12.2 m s -1 , and the maximum frequency of occurrence of the in-situ wind speed was located at the bin from 4.3 to 4.7 m s -1 (Figure 6d).  Figure 6a shows the spatial distribution of the number of matchup data at the buoy locations off the coast of the Korean peninsula from 1 April 2013 to 31 August 2017. The total number of the matchup data between the Landsat 8 OLI/TIRS images and the KMA ocean buoy data amounted to 320 and the largest number of matchups was 43, for the buoy Y1. The smallest number of the matchups was 6, for the buoys Y2 and S3 (Figure 6a). During the study period over five years, there were no matchups between the Landsat 8 data and the buoy data in some months (June 2013, June 2014, August 2014, January 2015, and June 2016) due to the mask associated with clouds, as shown in Figure 6b. The month-year plot of the matchup numbers in Figure 6b exhibits the maximum number of the matchups of 14 on April 2017. The in-situ temperatures of the matchup database ranged from 3.2 to 30.2 • C and the maximum frequency of occurrence of in-situ temperatures was in the range of 14.5 to 15.5 • C (Figure 6c). The KMA ocean buoys measure not only the sea water temperature at a depth of less than 0.4 m but also the wind speed. The in-situ wind speed of the matchup database ranged from 0.1 to 12.2 m s −1 , and the maximum frequency of occurrence of the in-situ wind speed was located at the bin from 4.3 to 4.7 m s −1 (Figure 6d).      Table 3 summarizes the accuracy and coefficients of the SST retrieval algorithm. Coefficients of the 11 µm BT term in the MCSST algorithms were about 0.97, closer to 1.0 than the NLSST algorithms. For the Landsat 8 OLI/TIRS, the coefficient of SZA term was about 32 to 35 for MCSST2, NLSST4, NLSST5, and NLSST6 algorithms, which was much higher than that of the MCSST and the NLSST algorithms using other algorithms for other satellites. Since the maximum value of SZA is up to about 8.4, the SZA term is close to 0. The high coefficients of these algorithms are believed to reflect the role of the SZA term in the SST retrieval equations. If the SZA term was neglected from the formulations, the accuracy of the NLSST with an RMSE of 0.61-0.66 • C was better than that of the MCSST with an RMSE of 0.72 • C. The cases where the SZA terms were considered, such as MCSST2, NLSST4, NLSST5, and NLSST6, were evaluated to be more accurate than those without the SZA effect considered in the formulation. For the first-guess SST, the input of OSTIA SST data was more accurate than the use of the estimated SST from the MCSST equation and MURSST. The accuracy of the SST derived from the Landsat 8 OLI/TIRS varied depending on the algorithms as well as the first-guess SST used as input data for the SST formulations. The RMSE values of the Landsat 8 OLI/TIRS given in Table 3 were quite similar to those from the numerous previous studies using other satellites [25,26,83].   (Figure 7). All the algorithms showed that the biases are close to zero without any noticeable trend (Figure 7, Table 3). As can be seen from the distribution of the number frequencies in the comparison plots, the scattering index (SI) of each algorithm (MCSST1, MCSST2, NLSST1, NLSST2, NLSST3, NLSST4, NLSST5, and NLSST6) was 0.0450, 0.0444, 0.0416, 0.0382, 0.0397, 0.0407, 0.0374, and 0.0387, respectively. Their correlation coefficients (R) were relatively high amounting to 0. 9933, 0.9935, 0.9943, 0.9952, 0.9948, 0.9945, 0.9954, and 0

High-Resolution SST
The Landsat 8 OLI/TIRS has an exceptionally high spatial resolution compared with other satellites used for SST retrieval. Using the coefficients for Landsat 8 derived in this study, SSTs were estimated for the southeastern coastal region of the Korean peninsula. These SSTs were compared with SSTs from the OSTIA, the MURSST, the GCOM-W1/AMSR-2, the Himawari-8 AHI, and the NOAA/AVHRR (Figure 8). The measured SSTs in S05 and S06 buoys were 13.7 • C and 14.8 • C at 02 UTC on 19 April 2019, respectively. Figure 8a shows the SST distribution from the OSTIA on 19 April 2016 and the collocated SSTs with each buoy were 14.1 • C and 14.6 • C, respectively. Figure 8b shows the SST distribution from the MURSST on 19 April 2016 and the collocated SSTs with each buoy were 14.0 • C and 14.2 • C, respectively. Figure 8c shows the SST distribution from the GCOM-W1/AMSR-2 at 05UTC on 19 April 2016. Due to its low spatial resolution of about 25 km and the effect of land on microwave measurements, there were only a few pixels with SSTs without any observations in the coastal area from 50 km to 100 km distance from the shoreline [84,85]. This implies that the microwave measurements are inappropriate to observe the SSTs of the coastal region. On the contrary, the Himawari-8 AHI produced SSTs with relatively high spatial resolution of about 2 km as shown in Figure 8d and the collocated SSTs with each buoy were 14.4 • C and 14.9 • C at 02 UTC on 19 April 2016, respectively. The NOAA/AVHRR data with a spatial resolution of about 1.1 km revealed more detailed structure of SST distribution including the contrast differences of temperatures at the frontal region between the onshore and the offshore regions and the collocated SSTs with each buoy were 14.3 • C and 14.8 • C at 05 UTC on 19 April 2016, respectively (Figure 8e). retrieval and the collocated SSTs with each buoy were 14.3 °C and 14.8 °C, respectively. Coastal SST patterns of the Landsat 8 (Figure 8f) illustrated more detailed spatial structures compared to other SST products (Figure 8a-e). Since the SSTs from the databases (AMSR-2, Himawari-8, NOAA-19, and Landsat 8) were observed at different times from 02 UTC to 05 UTC, they might be different due to diurnal variations of surface temperatures depending on the observation times during a day. Moreover, the composite products of OSTIA and MURSST represent average values of SSTs at that day. Nevertheless, most of the SSTs presented similar SSTs with relatively small differences of less than 1.6 °C. It is notable that the SSTs of Landsat 8 were well produced in the offshore region as well as along the coastline and around many complicated islands as shown in Figure 8f. This demonstrated the capability and usefulness of Landsat 8 data at the coastal regions including the estuarine regions. Thus, the comparison addresses the importance of using high-resolution SST data for coastal regions with very complicated temperature variations.

Effect of Atmospheric Moisture
The 11 μm channel is regarded to be a clean infrared atmospheric-window channel because it has limited absorption by atmospheric water vapor. On the contrary, the radiance at the 12 μm channel is relatively higher absorbed by atmospheric water vapor. Therefore, the BTD values between these two channels can reflect the absorption effect of the water vapor in the air [86]. The SST retrieval algorithms used these BTD terms to consider the effects of atmospheric water vapor on the estimation of SST. Figure 9 shows the SST errors; satellite SST minus in-situ temperature, as a The land pixels in Landsat 8 OLI/TIRS data were masked using SRTM DEM data prior to SST retrieval and the collocated SSTs with each buoy were 14.3 • C and 14.8 • C, respectively. Coastal SST patterns of the Landsat 8 ( Figure 8f) illustrated more detailed spatial structures compared to other SST products (Figure 8a-e). Since the SSTs from the databases (AMSR-2, Himawari-8, NOAA-19, and Landsat 8) were observed at different times from 02 UTC to 05 UTC, they might be different due to diurnal variations of surface temperatures depending on the observation times during a day. Moreover, the composite products of OSTIA and MURSST represent average values of SSTs at that day. Nevertheless, most of the SSTs presented similar SSTs with relatively small differences of less than 1.6 • C. It is notable that the SSTs of Landsat 8 were well produced in the offshore region as well as along the coastline and around many complicated islands as shown in Figure 8f. This demonstrated the capability and usefulness of Landsat 8 data at the coastal regions including the estuarine regions. Thus, the comparison addresses the importance of using high-resolution SST data for coastal regions with very complicated temperature variations.

Effect of Atmospheric Moisture
The 11 µm channel is regarded to be a clean infrared atmospheric-window channel because it has limited absorption by atmospheric water vapor. On the contrary, the radiance at the 12 µm channel is relatively higher absorbed by atmospheric water vapor. Therefore, the BTD values between these two channels can reflect the absorption effect of the water vapor in the air [86]. The SST retrieval algorithms used these BTD terms to consider the effects of atmospheric water vapor on the estimation of SST. Figure 9 shows the SST errors; satellite SST minus in-situ temperature, as a function of the BTD. One of conspicuous features is that both the MCSST1 and the MCSST2 algorithms resulted in obvious biases with significant underestimation of lower than −2.0 • C in the range of BTD greater than 2.0 K (Figure 9a,b). By contrast, these errors were considerably reduced in cases of the NLSST algorithms as shown in the algorithms (NLSST1, NLSST2, NLSST3, and NLSST4) from Figure 9c to Figure 9f. It is expected that the NLSST formulations would reduce the SST errors arising from the amount of water vapor in the air. Thus, the non-linear approaches based on the NLSST formulations are concluded to be better than the MCSST algorithms.
Remote Sens. 2019, 11, x FOR PEER REVIEW 13 of 21 function of the BTD. One of conspicuous features is that both the MCSST1 and the MCSST2 algorithms resulted in obvious biases with significant underestimation of lower than -2.0 °C in the range of BTD greater than 2.0 K (Figure 9a,b). By contrast, these errors were considerably reduced in cases of the NLSST algorithms as shown in the algorithms (NLSST1, NLSST2, NLSST3, and NLSST4) from Figure 9c to Figure 9f. It is expected that the NLSST formulations would reduce the SST errors arising from the amount of water vapor in the air. Thus, the non-linear approaches based on the NLSST formulations are concluded to be better than the MCSST algorithms.

Effect of SZA
As SZA increases, the atmospheric layer thickness through which the signal emitted from sea surface pass increases, affecting the intensity of the signal reaching the satellite. Accordingly, the SZA is one of the most important elements to accurately derive SST, as the atmospheric absorption effect is proportional to the SZA [87]. Most of the wide-scanning sensors such as the NOAA/AVHRR with a swath angle of 55 degrees in both sides have considered the effect of the SZA [25]. Since the Landsat sensors have very narrow viewing angles from -8.4 degrees to the left and 8.4 degrees to the right based on the satellite orbit path, the effect of viewing angle has been neglected in the SST retrieval procedures for a long time [46,[51][52][53][54]. However, the accuracy assessment of this study showed that the algorithms including the SZA term estimated SST with better performance. Regardless of the algorithms or the different use of the first-guess SST, the inclusion of the SZA reduced the SST errors.
To evaluate the performance rate of the improvement in SST when the SZA term is included, we estimated the percentage of the SST differences with respect to the SZA ( Figure 10). In case of the MCSST errors, the percentage was estimated with the errors, MCSST2 with the SZA term -MCSST1 without the SZA term divided by the MCSST2. Positive (Negative) angles of the SZA values correspond to the SZA on the right (left) side from the central ground track of the satellites. The enhanced percentages of the MCSST case showed a symmetric structure as indicated in Figure 10a. The maximum decrease was 0.75% in the range from -1 to 1 degree, and the maximum increase was 1.49% in the range from 7 to 9 degree. By adding the SZA term, the estimated SST was improved in the opposite direction based on ±4 degrees, resulting in an increase of about 0.1 °C in RMSE. The accuracies of the SSTs on the edges with high SZA of about -8 or 8 degrees improved over 1% by reaching 3%. Other cases of the NLSST formulations also revealed improved capability of the algorithms in reducing the errors when the SST formulations included the SZA term. This implies

Effect of SZA
As SZA increases, the atmospheric layer thickness through which the signal emitted from sea surface pass increases, affecting the intensity of the signal reaching the satellite. Accordingly, the SZA is one of the most important elements to accurately derive SST, as the atmospheric absorption effect is proportional to the SZA [87]. Most of the wide-scanning sensors such as the NOAA/AVHRR with a swath angle of 55 degrees in both sides have considered the effect of the SZA [25]. Since the Landsat sensors have very narrow viewing angles from −8.4 degrees to the left and 8.4 degrees to the right based on the satellite orbit path, the effect of viewing angle has been neglected in the SST retrieval procedures for a long time [46,[51][52][53][54]. However, the accuracy assessment of this study showed that the algorithms including the SZA term estimated SST with better performance. Regardless of the algorithms or the different use of the first-guess SST, the inclusion of the SZA reduced the SST errors.
To evaluate the performance rate of the improvement in SST when the SZA term is included, we estimated the percentage of the SST differences with respect to the SZA ( Figure 10). In case of the MCSST errors, the percentage was estimated with the errors, MCSST2 with the SZA term -MCSST1 without the SZA term divided by the MCSST2. Positive (Negative) angles of the SZA values correspond to the SZA on the right (left) side from the central ground track of the satellites. The enhanced percentages of the MCSST case showed a symmetric structure as indicated in Figure 10a. The maximum decrease was 0.75% in the range from −1 to 1 degree, and the maximum increase was 1.49% in the range from 7 to 9 degree. By adding the SZA term, the estimated SST was improved in the opposite direction based on ±4 degrees, resulting in an increase of about 0.1 • C in RMSE. The accuracies of the SSTs on the edges with high SZA of about −8 or 8 degrees improved over 1% by reaching 3%. Other cases of the NLSST formulations also revealed improved capability of the algorithms in reducing the errors when the SST formulations included the SZA term. This implies that the SST formulation should include the SZA to estimate more accurate SST from the Landsat 8 data in spite of its narrow swath width.
Remote Sens. 2019, 11, x FOR PEER REVIEW 14 of 21 that the SST formulation should include the SZA to estimate more accurate SST from the Landsat 8 data in spite of its narrow swath width.

Effect of Wind Speed
Although we derived the SST coefficients by regressing the satellite BT values to the buoy temperatures, the SST is different from the skin SST measured at the thin layer that is affected by the heat transfer and diffusion from the atmosphere ocean boundary layer [88]. Since the temperature sensors of the buoys are located at different depths from 0.2 to 0.4 m, the retrieved SST might have a possibility to be affected by the winds, as pointed out by the previous studies [89][90][91][92]. Therefore, we tested if the derived SSTs from the Landsat 8 have any relation to the wind speed.

Effect of Wind Speed
Although we derived the SST coefficients by regressing the satellite BT values to the buoy temperatures, the SST is different from the skin SST measured at the thin layer that is affected by the heat transfer and diffusion from the atmosphere ocean boundary layer [88]. Since the temperature sensors of the buoys are located at different depths from 0.2 to 0.4 m, the retrieved SST might have a possibility to be affected by the winds, as pointed out by the previous studies [89][90][91][92]. Therefore, we tested if the derived SSTs from the Landsat 8 have any relation to the wind speed. Figure 11 shows the difference between buoy temperature and SST derived from the Landsat 8 OLI/TIRS according to the wind speed, for the eight cases from MCSST1 to NLSST6. Regardless of the algorithms or the different first-guess SSTs, no distinct bias errors were observed in the range of wind speeds higher than 2 m s −1 . However, the Landsat 8 OLI/TIRS SSTs tended to be overestimated in the range of wind speeds below 2 m s −1 as indicated in all plots in Figure 11 [92,93]. Overall, an average positive bias of 0.63 • C was found in the range of wind speeds below 2 m s −1 .
Since the Landsat satellite observes the sea surface at 11 to 12 AM on local time, the difference between the bulk SST and the skin SST could be directly affected by the solar thermal energy, which causes the increase in the skin SST at this time [94]. At low wind speeds, the difference between the bulk SST and the skin SST is large because of the weak mixing in the upper marine layer, however at high wind speeds, the mixing of the upper layer is active and the difference between the skin SST and the bulk SST decreases [95]. This trend is evident in the previous validation studies of other satellite-derived SSTs [26,59,96]. In this study, it was confirmed that the error characteristics of SST derived from the Landsat 8 OLI/TIRS were clearly shown at daytime with low wind speed. This addresses the significant imprint of low wind condition on satellite SSTs, which produces the overestimation, especially in a range of less than 2 m s −1 .

Effect of Wind Speed
Although we derived the SST coefficients by regressing the satellite BT values to the buoy temperatures, the SST is different from the skin SST measured at the thin layer that is affected by the heat transfer and diffusion from the atmosphere ocean boundary layer [88]. Since the temperature sensors of the buoys are located at different depths from 0.2 to 0.4 m, the retrieved SST might have a possibility to be affected by the winds, as pointed out by the previous studies [89][90][91][92]. Therefore, we tested if the derived SSTs from the Landsat 8 have any relation to the wind speed.

Summary and Conclusions
In order to calculate the SST coefficients for the Landsat 8 OLI/TIRS, this study produced the matchup database between the Landsat 8 data and the KMA marine meteorological buoy measurements for the coastal region around the Korean Peninsula. A total number of 276 Landsat 8 OLI/TIRS images were used to produce 320 matchup data for the period from April 2013 to August 2017. The curvature of the Earth was taken into account for estimating more accurate SZA with a range of about ±8.4 degrees, which was a little wider than the angles calculated by the previous studies under the assumption of a flat sea surface, within the relatively narrow swath of the Landsat 8. The SST algorithms were classified into MCSST1, MCSST2, NLSST1, NLSST2, NLSST3, NLSST4, NLSST5, and NLSST6 depending on the uses of SZA terms or the first-guess SST. These algorithms produced RMSEs of 0.72 • C, 0.71 • C, 0.66 • C, 0.61 • C, 0.63 • C, 0.65 • C, 0.59 • C, and 0.62 • C, respectively. Among the algorithms, the smallest error in terms of RMSE and bias was found for the NLSST5 using SZA, with OSTIA SST as the first-guess SST rather than MCSST and MURSST.
Characteristics of the retrieved SSTs were analyzed by investigating the effects of the SZA, the absorption of atmospheric water vapor, and the wind speeds on the SST. The satellite-derived MCSSTs had a tendency to be underestimated at relatively high moist condition of greater than 2 K as inferred from the differences between the 11 µm and 12 µm channel BTs. Such errors of MCSSTs were reduced in the case of the NLSST formulations by applying the SZA and the use of the proper first-guess SST. The SST errors revealed the quadratic dependence on the SZA values from the nadir to the edge of the swath. After the inclusion of SZA term in the SST estimation, the SST results were moderately improved by less than 3% as compared with prior equations without the angle term. The error analysis due to the wind speeds showed that, the SST errors were amplified as the wind speeds became weak at a range of less than 2 m s −1 . This has been reported to be associated with the skin-bulk temperature differences and solar insolation under weak wind condition [86,89]. When the SST coefficients calculated from the matchup data from April 2013 to August 2016 were applied to the data from September 2016 to August 2017, the RMSE was less than about 0.7 • C. Therefore, it is expected that similar accuracy would be maintained even if the SST coefficients calculated in this study were applied to the Landsat 8 OLI/TIRS data; subsequent to other periods or other coastal regions.
The previous Landsat satellites observed using one thermal infrared band, but the Landsat 8 OLI/TIRS was classified into 11 and 12 µm bands in the thermal infrared band. Therefore, the most commonly used MCSST and NLSST for SST retrieval are available, but SST coefficients optimized for the Landsat 8 OLI/TIRS data have not been developed. In addition, this study considered the effect of the SZA in the SST retrieval formulations. The repeat cycle of the Landsat with a small swath width is much longer than other polar-orbiting satellites (e.g., MODIS and VIIRS) with complete global coverage every day. Moreover, it could be difficult for SST retrieval using Landsat to be applied to the areas with regional characteristics of intense cloud cover such as the coastal regions off California, Peru, or Chile. Nevertheless, by this study, when the SST algorithm suitable for Landsat 8 OLI/TIRS is established, we expect to be able to study sub-kilometer small scale ocean phenomena in the coastal region.
Since the Korean coastal region has a variety of oceanic phenomena in terms of time and space scales, with fronts, eddies, tidal currents, river discharge, and estuary, the SST coefficients derived in this study might be applicable to SST derivation for other coastal regions of the global ocean as well, especially for the regions without any local SST coefficients for the Landsat 8 data. For more extensive and operational use of the Landsat 8 OLI/TIRS derived SSTs, it is important to continuously monitor and understand the error characteristics of the SST using in-situ measurements in diverse ocean regions. If the validation studies in various sea areas are conducted with satisfactory results, it is expected that the present SST algorithms and coefficients could be further improved in the future and be extensively used for studying small-scale ocean phenomena in the coastal regions including estuarine regions. Based on the local error characteristics in these studies, it is possible to reduce the regional SST errors through regional optimization for accurate SST retrieval procedures.