Assessment and Improvement of Global Gridded Sea Surface Temperature Datasets in the Yellow Sea Using In Situ Ocean Buoy and Research Vessel Observations

: The sea surface temperature (SST) is essential data for the ocean and atmospheric prediction systems and climate change studies. Five global gridded sea surface temperature products were evaluated with independent in situ SST data of the Yellow Sea (YS) from 2010 to 2013 and the sources of SST error were identiﬁed. On average, SST from the gridded optimally interpolated level 4 (L4) datasets had a root mean square di ﬀ erence (RMSD) of less than 1 ◦ C compared to the in situ observation data of the YS. However, the RMSD was relatively high (2.3 ◦ C) in the shallow coastal region in June and July and this RMSD was mostly attributed to the large warm bias ( > 2 ◦ C). The level 3 (L3) SST data were frequently missing in early summer because of frequent sea fog formation and a strong ( > 1.2 ◦ C / 12 km) spatial temperature gradient across the tidal mixing front in the eastern YS. The missing data were optimally interpolated from the SST observation in o ﬀ shore warm water and warm biased SST climatology in the region. To fundamentally improve the accuracy of the L4 gridded SST data, it is necessary to increase the number of SST observation data in the tidally well mixed region. As an interim solution to the warm bias in the gridded SST datasets in the eastern YS, the SST climatology for the optimal interpolation can be improved based on long-term in situ observation data. To reduce the warm bias in the gridded SST products, two bias correction methods were suggested and compared. Bias correction methods using a simple analytical function and using climatological observation data reduced the RMSD by 19–29% and 37–49%, respectively, in June.


Introduction
Sea surface temperature (SST) is one of the main physical variables that provides information regarding the current state of the ocean [1]. SST data have been used for data assimilation in ocean circulation models and as a bottom boundary condition for atmospheric prediction models [2][3][4]. They are also essential for climate modeling and ocean-atmosphere interaction studies [5]. Satellite-based SST datasets have been produced by various organizations and are available at large spatial scales and at small temporal intervals [6][7][8][9]. Several previous studies have compared these global SST products with in situ observation data. Evaluations of the global SST datasets have mostly been performed at

Satellite SST Products
The main characteristics of the five L4 gridded SST datasets are summarized in Tables 1 and 2. The datasets have different horizontal grid spacing and represent different SST types. The OISST data were generated by the National Climatic Data Center (NCDC) of National Oceanic and Atmospheric Administration (NOAA) with a spatial grid spacing of 0.25 • and a daily temporal resolution [27]. The OISST product uses operational Advanced Very High Resolution Radiometer (AVHRR) satellite data. Both Infrared (IR)-Only and Infrared-Microwave (IR-MW) products are available for the OISST datasets. The IR-Only OISST dataset was used in this study. The MGDSST data were provided by the Japan Meteorological Agency (JMA) using an optimal interpolation (OI) technique. The MGDSST dataset also has a spatial grid resolution of 0.25 • and a daily temporal resolution. The product uses operational AVHRR and AMSR-E satellite data [28]. The OSTIA dataset was produced by the UK Met Office and uses satellite observation data together with in situ observations to determine the SST [29]. The SST analysis was performed using a variant of the OI technique described by Reynolds et al. [27]. The OSTIA SST has a daily temporal scale and a resolution of 0.05 • . The MWIR SST dataset uses the OI technique on SST data from the operational Tropical Rainfall Measuring Mission Microwave Imager (TMI), AMSR-E and WindSAT microwave sensors as well as the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra and MODIS Aqua IR sensors [30,31]. The MWIR SST dataset has a spatial grid resolution of 9 km and a daily temporal resolution. The GMPE system was developed and is operated by the UK Met Office. The dataset consists of a median and standard deviation of a daily ensemble at a resolution of 0.25 • . Twelve SST datasets contributed to the ensemble median [6]. Both daytime and nighttime SST measurements were used for the L4 gridded SST datasets, and the diurnal warming was removed from the daytime SST using different methods [27][28][29][30]. The AVHRR Pathfinder Version 5.3 (PFV53) Level 3 (L3C) SST dataset is a collection of global, twice-daily (Day and Night) 4 km SST data produced by the NOAA National Centers for Environmental Information (NCEI). The PFV53 dataset was computed with data from the AVHRR instruments on board NOAA's polar orbiting satellite series using an entirely modernized system based on SeaDAS. The PFV53 data were collected through the operational periods of the NOAA-7 through NOAA-19 Polar Operational Environmental Satellites, and are available from 1981 to the present. The L3C SST dataset for the period from 2010 to 2013 was examined in this study.
The Multifunctional Transport Satellites (MTSAT) is a Japanese geostationary dual-function satellite program, procured by the Japan Civil Aviation Bureau and the JMA. The satellite data from the MTSAT has short-wave infrared (SWIR, 3.8 µm) and long-wave infrared (LWIR, 10.8 µm) channels, which were used to monitor the area of fog. The MTSAT has five channels that consist of four IR window channels and one visible window channel, and makes observations using a 4 km IR resolution, which is sufficient for monitoring sea fog. The MTSAT data consist of snapshots taken at 30 min intervals. Images from the loaded 3.8 µm sensor and the IR differential images (3.8-10.8 µm) from the MTSAT can detect nighttime fog and show its distribution over a wide area. The sea fog images from the MTSAT for the years of 2010-2013 were used to detect sea fog in this study.

In Situ SST Data in the Yellow Sea
The in situ temperature observation data were obtained from ship CTD profiles and ocean buoys in the YS to evaluate the five L4 gridded SST datasets. The CTD observation data were provided by the National Institute of Fisheries Science (NIFS) of Korea; the data were obtained along 25 lines and 207 stations in the Eastern YS ( Figure 1). The surface temperature was measured six times per year. Another set of CTD observation data was provided by the Korea Marine Environment Management Corporation (KOEM). These CTD data were obtained from 417 stations off the Korean coastlines, measured four times per year. The surface observation depth of the CTD profiles was 1.5-2.0 m and these measurements were used as the CTD SSTs in this study.  Surface buoy D is located near Deokjeok Island and the surface buoy O is located near Oeyeon Island. The SST data measured by buoys D and O were provided by the Korea Meteorological Administration (KMA). The buoys measured the SST every hour and the daily mean SSTs were calculated using the hourly data. Buoy KB was located in the middle of the YS. The SST data at KB was measured every 10 min and daily mean SSTs were calculated from the 10 min data. Surface buoy SB was located near the southern boundary of the YS and hourly mean data from this buoy were used to produce the daily mean data (Figure 1). The buoys measured the temperatures at a depth of 0.5-1.0 m. A total of 1249 and 2031 SST observations were collected from the CTD profiles and ocean buoys, respectively. The SST data from the ocean buoys (KB and SB) and the bi-monthly routine hydrographic observation stations were independent data and they were not used to produce the global gridded SST datasets, whereas the OSTIA and MGDSST datasets were produced using SST data from buoys D and O. The validations were performed by in situ datasets, and confidence intervals for statistical results were evaluated by the bootstrapping procedure [32].

Bias Correction Method of the L4 SST Products
To reduce the high RMSD in the coastal region of the Eastern YS, two spatial-temporal bias correction methods were designed in this study. In the first method, the bias correction function (∆T) consisted of an exponential function depending on distance from the land and a cosine function depending on time as follows: where t represents the number of days from January 1 of each year and r is the distance from the land in km. The SST bias correction function was designed to fit the five L4 SST datasets with the in situ surface temperature data. A 0 is the maximum bias near the coast in June and was calculated to be 3.64, 3.21, 2.65, 5.33, and 4.10 • C for the OISST, MGDSST, OSTIA, MWIR and GMPE SST, respectively. In the second bias correction method, the bi-monthly horizontal distributions of the in situ CTD observation data from 2010 to 2013 were used to construct daily bias correction functions from June to July. First, the monthly climatologies were constructed using the CTD observation and the original L4 gridded SST data. Daily climatologies were then calculated from the monthly climatologies by linear interpolation in time. The daily bias correction function (∆T) was defined as the difference between the two daily climatology datasets on day t: ∆T(lon, lat, t) = SST Level4 (lon, lat, t) − SST CTD (lon, lat, t), 120 < t < 240 (2) where t represents the number of days from January 1 of each year. The overbar represents the average of the four years at (lon, lat) on day t and (lon, lat) represents the horizontal location. The daily bias correction function at time t was calculated based on the quantitative comparison of the daily OISST, MGDSST, OSTIA, MWIR and GMPE SST climatologies with the daily in situ CTD SST climatology.

OSTIA SST Interpolation Scheme
The following is a short description of the OSTIA SST analysis procedure [29]. The background field x b i, k at grid point i and time k is defined as: where λ i, k is a scalar less than 1, x a i,k−1 is the previous analysis, and x c i,k is a reference climatology valid for the same time of year at time k. Relaxation time scale at grid point i and time k is 30 days for the determination of λ i, k . The background field was calculated using Equation (3) and the bias corrected measurements were then used to produce an analysis using a multi-scale OI type scheme. The OI equation is given by: where x k is a vector containing all values x i, k , i = 1, . . . , N on the analysis grids at time k, y k is a vector containing the observations, B is the background error covariance matrix, R k is the observation error covariance matrix, and H k is the observation operator that interpolates from the model grid to the observation locations. This equation was solved using the analysis correction method [33,34]. The background error covariance matrix B was calculated using model outputs from a three-year numerical integration using the Forecasting Ocean Assimilation Model (FOAM) of the UK Met Office [35,36]. Covariances were split into two components, and two error correlation length scales of 10 km and 100 km were used for mesoscale and large scale errors, respectively. The data in the observation error covariance matrix R were assumed to be uncorrelated with each other and were obtained from the information supplied by the GHRSST single sensor error statistic values. A simple quadratic bottom friction (C d ) value of 0.00125 was used to crudely parametrize tidal mixing in the FOAM [37].

Comparison with In Situ SST Observations
To evaluate the accuracy of L4 gridded SST datasets (OISST, MGDSST, OSTIA, MWIR, and GMPE) from five operational systems, the RMSD and bias of the gridded SST datasets were calculated relative to the in situ SST obtained from the ocean buoys and routine hydrographic observation stations in the Eastern YS. The correlation coefficients between the buoy and L4 gridded SST data were 0.97-0.99 ( Figure 2). The RMSD of the gridded SST data was the lowest (0.74 • C) for the GMPE and the highest (1.04 • C) for the OISST. The blue dots lying above the center line of the OISST and GMPE datasets in August imply that the L4 gridded SSTs were higher on some days than the in situ SSTs. The L4 gridded SSTs in June (red dots) were aligned with the in situ ocean buoy SSTs from the deep stratified region.
Compared with the CTD SST data, some of the red and blue dots were above the center line ( Figure 3), which indicated that the SST from the satellite-based observations were warmer than that from the in situ CTD data in June and August. The correlation coefficients were 0.93-0.98. The RMSD of the satellite observations was the lowest (0.96 • C) for the OSTIA SST and the highest (1.39 • C) for the MWIR SST. These comparisons suggest that the SST from the OSTIA, GMPE, and MGDSST datasets had a RMSD of less than 1.00 • C compared with the in situ CTD observations in the Eastern YS (Table 3). Compared with the CTD SST data, some of the red and blue dots were above the center line ( Figure 3), which indicated that the SST from the satellite-based observations were warmer than that from the in situ CTD data in June and August. The correlation coefficients were 0.93-0.98. The RMSD of the satellite observations was the lowest (0.96 °C) for the OSTIA SST and the highest (1.39 °C) for the MWIR SST. These comparisons suggest that the SST from the OSTIA, GMPE, and MGDSST datasets had a RMSD of less than 1.00 °C compared with the in situ CTD observations in the Eastern YS (Table 3). from the in situ CTD data in June and August. The correlation coefficients were 0.93-0.98. The RMSD of the satellite observations was the lowest (0.96 °C) for the OSTIA SST and the highest (1.39 °C) for the MWIR SST. These comparisons suggest that the SST from the OSTIA, GMPE, and MGDSST datasets had a RMSD of less than 1.00 °C compared with the in situ CTD observations in the Eastern YS (Table 3).   Table 3. Root mean square difference (RMSD) and bias of SST data obtained from in situ (the ocean buoys and conductivity-temperature-depth (CTD) measurements and remote sensing (satellite observation). The smallest RMSD and bias relative to KB, SB, and NIFS observations are given in bold face. The number of matchup data is given in parentheses.

Spatial Distribution of the SST RMSD
The spatial distribution of the RMSD between the SST from the global gridded products and those from the in situ research vessel observations was calculated at each observation station over the four years ( Figure 4). About 24 values of SST data were compared at each station (black dot in Figure 4). The RMSD between the satellite-based SSTs and the CTD SSTs was found to be comparatively higher in the shallow and vertically well-mixed region or in the shallow coastal region within 50 km from the land and lower in the offshore region. The RMSD of the MWIR SST (Figure 4d) was the highest and that of the OSTIA SST was the lowest (Figure 4c). The annual mean RMSDs for the shallow coastal region within 50 km from the land were 1.50, 1.17, 1.14, 1.72, and 1.27 • C for the OISST, MGDSST, OSTIA, MWR, and GMPE datasets, respectively ( Figure 4 and Table 4). The RMSDs were larger in June-August than in other months (Table 4). Spatial distribution of the four-years mean CTD SST in June showed relatively cooler coastal water (T < 18 • C) near the coast (Figure 4f). The cool water bounded by approximately 18 • C isotherm is located in the shallow coastal region within 50 km from the land due to strong tidal stirring, which forms a tidal mixing front along the 30-50 m isobaths. This region coincided with the high RMSD area of the L4 gridded SST datasets. The annual mean RMSDs for the offshore region were 0.88, 0.73, 0.64, 0.74, and 0.59 • C for the OISST, MGDSST, OSTIA, MWIR, and GMPE, respectively ( Figure 4 and Table 4). The RMSD of the GMPE SST was the lowest and that of the OISST was the highest. Table 4. Seasonal variations in the RMSD between the L4 gridded SST products (mostly based on satellite observation) and the in situ measurements (ocean buoys and research vessels) in the shallow and vertically well-mixed region. The lowest RMSD in each column is given in bold face. The annual mean RMSDs for the offshore region were 0.88, 0.73, 0.64, 0.74, and 0.59 °C for the OISST, MGDSST, OSTIA, MWIR, and GMPE, respectively ( Figure 4 and Table 4). The RMSD of the GMPE SST was the lowest and that of the OISST was the highest.

Seasonal Variation of the RMSD and Bias
There were relatively large differences between the L4 gridded SSTs and the in situ SSTs in June and August ( Figure 3). To examine seasonal the variations in the RMSD, monthly mean RMSD of the L4 gridded SST data were calculated at the ocean buoy stations (KB, SB, D, and O), the NIFS observation stations in the deep stratified region (NIFS-S) where the bottom depth was deeper than 50 m, and the NIFS observation stations in the shallow vertically mixed region (NIFS-M) where the bottom depth was shallower than 50 m, for four years ( Figure 5). The ocean buoys at KB and SB were located in the relatively deep stratified region while the ocean buoy at D was located in the relatively shallow and vertically well-mixed region. The ocean buoy at O was located near the boundary between the stratified and the mixed regions ( Figure 1). In general, the RMSDs of the OISST, MGDSST, OSTIA, and GMPE SST datasets were higher from June to August in the Eastern YS. However, the RMSD of the MWIR dataset did not have a clear seasonal variation in the deep stratified region. The RMSD was higher in the shallow and vertically well-mixed region (D and NIFS-M) from June to July (Figure 5e,f). In June and July, a cold bias in the global gridded SST data appeared in the deep offshore region such as at the KB, SB, and NIFS-S stations, while a warm bias (1-4 • C) appeared in the shallow coastal region, such as at the D and NIFS-M stations ( Figure 6). The relatively large RMSD in the shallow and vertically well-mixed region was related to the warm bias in June-July. In the shallow mixed region, the annual mean biases of the OISST, MGDSST, OSTIA, MWIR, and GMPE SST data compared with the CTD SST data in June were 1.96, 1.73, 1.17, 2.28, and 2.27 • C, respectively (Figure 6f).

Seasonal Variation of the RMSD and Bias
There were relatively large differences between the L4 gridded SSTs and the in situ SSTs in June and August (Figure 3). To examine seasonal the variations in the RMSD, monthly mean RMSD of the L4 gridded SST data were calculated at the ocean buoy stations (KB, SB, D, and O), the NIFS observation stations in the deep stratified region (NIFS-S) where the bottom depth was deeper than 50 m, and the NIFS observation stations in the shallow vertically mixed region (NIFS-M) where the bottom depth was shallower than 50 m, for four years ( Figure 5). The ocean buoys at KB and SB were located in the relatively deep stratified region while the ocean buoy at D was located in the relatively shallow and vertically well-mixed region. The ocean buoy at O was located near the boundary between the stratified and the mixed regions ( Figure 1). In general, the RMSDs of the OISST, MGDSST, OSTIA, and GMPE SST datasets were higher from June to August in the Eastern YS. However, the RMSD of the MWIR dataset did not have a clear seasonal variation in the deep stratified region. The RMSD was higher in the shallow and vertically well-mixed region (D and NIFS-M) from June to July (Figure 5e,f). In June and July, a cold bias in the global gridded SST data appeared in the deep offshore region such as at the KB, SB, and NIFS-S stations, while a warm bias (1-4 °C) appeared in the shallow coastal region, such as at the D and NIFS-M stations ( Figure 6). The relatively large RMSD in the shallow and vertically well-mixed region was related to the warm bias in June-July. In the shallow mixed region, the annual mean biases of the OISST, MGDSST, OSTIA, MWIR, and GMPE SST data compared with the CTD SST data in June were 1.   To examine the seasonal characteristics of each SST dataset, the RMSDs between the L4 gridded SST datasets and the in situ data were divided for MM (March-May), JA (June-August), SN (September-November), and DF (December-February), and were examined in the vertically wellmixed shallow region and the stratified region (Table 4 and Table 5). In the shallow and vertically well-mixed region, the RMSD between the MGDSST and the in situ CTD SST was the lowest in MM and DF while the RMSD between the OSTIA SST and the in situ CTD SST was the lowest in JA and SN (Table 4). This indicates that the MGDSST and OSTIA SST datasets are the best for the shallow mixed region of the Eastern YS throughout the year.   To examine the seasonal characteristics of each SST dataset, the RMSDs between the L4 gridded SST datasets and the in situ data were divided for MM (March-May), JA (June-August), SN (September-November), and DF (December-February), and were examined in the vertically well-mixed shallow region and the stratified region (Tables 4 and 5). In the shallow and vertically well-mixed region, the RMSD between the MGDSST and the in situ CTD SST was the lowest in MM and DF while the RMSD between the OSTIA SST and the in situ CTD SST was the lowest in JA and SN (Table 4). This indicates that the MGDSST and OSTIA SST datasets are the best for the shallow mixed region of the Eastern YS throughout the year. In the stratified offshore region, the RMSD between the GMPE SST and the in situ CTD SST was the lowest in MM, JA, and DF, and that between the OSTIA SST and the in situ CTD SST was the lowest in SN (Table 5). When the OISST (IR band product) and the MWIR (microwave band product) were compared with the in situ data in June-August, the RMSDs of the OISST and MWIR were 1.14-1.45 • C and 0.84-0.85 • C, respectively, in the deep stratified region (Table 5). This implies that the microwave band observations are more effective in the stratified region during the summer than the IR band observation (Figure 5a-d).

Bias Correction of the L4 SST Products in the YS
The RMSDs of the OSTIA SST data were 1.14 and 0.64 • C in the shallow mixed and deep stratified regions, respectively. The relatively high RMSD in the shallow and vertically well-mixed region in June and July was largely contributed by the warm bias (Figures 5e,f and 6e,f). To improve the accuracy of the L4 gridded SST data in June and July, two bias correction methods were developed as described in Section 2.3. After applying the first bias correction method (Equation (1)) based on the analytical functions of the OISST dataset, the bias corrected OISST (OISSTabc) data were compared with the SST from ocean buoy D and NIFS-M (Figure 7). The OSTIA dataset was also corrected using the bias correction function and the bias corrected OSTIA (OSTIAabc) data were compared with NIFS-M SST data. The OSTIA SST dataset was produced using SST data from ocean buoy D, therefore, it was not compared with the SST data from this buoy. The RMSDs of the OISST and OISSTabc datasets were compared with the in situ SST data from buoy D and showed that the RMSD was reduced by 54% in May-August. The RMSD of the OISST (OSTIA) dataset was decreased by 28% (19%) with the bias correction at NIFS-M stations in June (Figure 7). Although not shown in the figure, the RMSD of the MGDSST, MWIR and GMPE datasets decreased by 23%, 26% and 29%, respectively, at NIFS-M stations in June. After the second bias correction method (Equation (2)) was applied to the L4 gridded SST datasets, the bias corrected SST (OISSTrbc and OSTIArbc) data were compared with the in situ SST from the NIFS-M and ocean buoy D (Figure 8). Comparing the RMSDs of the original OISST and bias corrected OISSTrbc data with the in situ SST data from buoy D, the RMSD was reduced by 58% in May-August. The RMSD of the OISST (OSTIA) dataset decreased by 41% (37%) with the bias correction at NIFS-M stations in June. Although not shown in the figure, the RMSD of the MGDSST, MWIR and GMPE datasets decreased by 43%, 47% and 49%, respectively, at NIFS-M stations in June.
in situ SST from the NIFS-M and ocean buoy D (Figure 8). Comparing the RMSDs of the original OISST and bias corrected OISSTrbc data with the in situ SST data from buoy D, the RMSD was reduced by 58% in May-August. The RMSD of the OISST (OSTIA) dataset decreased by 41% (37%) with the bias correction at NIFS-M stations in June. Although not shown in the figure, the RMSD of the MGDSST, MWIR and GMPE datasets decreased by 43, 47 and 49%, respectively, at NIFS-M stations in June.  When the second bias correction method was applied to the L4 gridded OSTIA SST dataset in June, the RMSD reduced in both the shallow coastal region and the deep stratified region (Figure 9). This suggests that there were large warm biases in the background field of the L4 gridded SST datasets in June (Figure 4 and Figure 9) and these could be corrected using the long-term climatology  When the second bias correction method was applied to the L4 gridded OSTIA SST dataset in June, the RMSD reduced in both the shallow coastal region and the deep stratified region (Figure 9). This suggests that there were large warm biases in the background field of the L4 gridded SST datasets in June (Figure 4 and Figure 9) and these could be corrected using the long-term climatology When the second bias correction method was applied to the L4 gridded OSTIA SST dataset in June, the RMSD reduced in both the shallow coastal region and the deep stratified region (Figure 9). This suggests that there were large warm biases in the background field of the L4 gridded SST datasets in June (Figures 4 and 9) and these could be corrected using the long-term climatology of in situ CTD SST observation data for the Eastern YS. The bias correction reduced not only bias but also the RMSD for the Eastern YS (Figures 7-9).
Remote Sens. 2020, 12, 759 15 of 25 of in situ CTD SST observation data for the Eastern YS. The bias correction reduced not only bias but also the RMSD for the Eastern YS (Figure 7, Figure 8 and Figure 9).

Discussion
The RMSD between the L4 gridded SST and in situ SST had a seasonal variation and it was relatively large in June and July. The large RMSD in June and July ( Figure 5) was contributed by the warm bias in the shallow and vertically well-mixed region ( Figure 6). Since the bias contributes to a large portion of the RMSD in the YS during summer, it is important to find the causes that induce the SST bias in the L4 gridded SST datasets. There were several error (the RMSD and bias) sources of the L4 gridded SST products in the coastal mixed regions.
To identify oceanographic conditions that affect the relatively large RMSD and bias in the Eastern YS, (i) spatial distribution of the warm bias in the L4 gridded SST datasets, (ii) frequent sea fog formation, and (iii) the number of available L3 SST data were examined in this section. (iv) How

Discussion
The RMSD between the L4 gridded SST and in situ SST had a seasonal variation and it was relatively large in June and July. The large RMSD in June and July ( Figure 5) was contributed by the warm bias in the shallow and vertically well-mixed region ( Figure 6). Since the bias contributes to a large portion of the RMSD in the YS during summer, it is important to find the causes that induce the SST bias in the L4 gridded SST datasets. There were several error (the RMSD and bias) sources of the L4 gridded SST products in the coastal mixed regions.
To identify oceanographic conditions that affect the relatively large RMSD and bias in the Eastern YS, (i) spatial distribution of the warm bias in the L4 gridded SST datasets, (ii) frequent sea fog formation, and (iii) the number of available L3 SST data were examined in this section. (iv) How to reduce the bias during the optimal interpolation procedure of the L4 gridded SST is sought for the Eastern YS.

Spatial Pattern of Warm Bias in the Shallow Region
To find the characteristics of the warm bias in the L4 gridded SST data, spatial distributions of the SST from the OSTIA, OISST, and MWIR dataset were compared with that of the SST from CTD profiles in the Eastern YS for June from 2010 to 2013 ( Figure 10). The SST distributions from the OISST and MWIR dataset were distinctly different and had a warm bias in the tidally well-mixed and cool water region. However, the OSTIA SST had a smaller bias than the other SSTs because it uses the in situ data from the ocean buoys during the OI process (Figure 10b). The original SST rather than the bias was plotted in Figure 10 to show both spatial pattern and horizontal gradient of the SST from the L4 gridded datasets. The horizontal gradient was higher in the in situ SST than the L4 gridded SSTs. The spatial pattern of bias were similar to those of RMSD in Figure 4.
Remote Sens. 2020, 12, 759 16 of 25 to reduce the bias during the optimal interpolation procedure of the L4 gridded SST is sought for the Eastern YS.

Spatial Pattern of Warm Bias in the Shallow Region
To find the characteristics of the warm bias in the L4 gridded SST data, spatial distributions of the SST from the OSTIA, OISST, and MWIR dataset were compared with that of the SST from CTD profiles in the Eastern YS for June from 2010 to 2013 ( Figure 10). The SST distributions from the OISST and MWIR dataset were distinctly different and had a warm bias in the tidally well-mixed and cool water region. However, the OSTIA SST had a smaller bias than the other SSTs because it uses the in situ data from the ocean buoys during the OI process (Figure 10b). The original SST rather than the bias was plotted in Figure 10 to show both spatial pattern and horizontal gradient of the SST from the L4 gridded datasets. The horizontal gradient was higher in the in situ SST than the L4 gridded SSTs. The spatial pattern of bias were similar to those of RMSD in Figure 4.  The combined effect of tidal mixing and stratification produces the tidal mixing front along the 30-50 m depth isobaths and relatively cool water occupies the shallow coastal region in summer. The OISST data are based on satellite IR band observations, which are prone to cloud and fog shielding, and underestimate tidal cooling over the shallow mixed region [38]. Using microwave sensors, the AMSR-E can measure SST through clouds and capture the tidal cooling effect well. However, the microwave sensors have a large observation error near land and islands. The MWIR SST also had a warm bias in the shallow mixed region (Figure 10d). This was because of the presence of land inference errors due to the characteristic natures of the microwave sensors [31]. Xie et al. [8] accessed satellite-based SST data in the shelf and coastal seas around China and suggested a depth-dependent bias correction method. Using the method, the bias in the SST data in the region shallower than 40 m depth was significantly corrected but the reduction in the RMSD was relatively small (6-14%). In their bias correction method the seasonal variation of the SST bias was not considered.

Sea Fog Formation in the Eastern YS
Sea fog was frequently observed as the cold patches along the Kuril Islands from July to August due to a tidal cooling effect, and the AVHRR SST was lower than AMSR-E SST [38]. The SST based on the satellite IR sensors in the Kashevarov Bank has cloud and fog shielding, which results in a warm bias larger than 5 • C and an underestimation of tidal cooling. Satellites with a microwave sensor can measure the SST through the clouds and also in the regions under the influence of strong tidal mixing. The satellites with IR sensors could not observe SST in the tidally well-mixed and cool water region because of lower clouds and sea fog in the Okhotsk Sea [38].
Inter-comparison of the in situ (ship and buoy) SST data with satellite-based SST products for the shelf and coastal seas around China and in the Northwest Pacific over a one-year period indicate that the RMSD increases sharply in the coastal region [3]. A comparison of the OISST data with the CTD SST observation data over 30 years for the Eastern YS revealed that the RMSD was 2-3 • C in the coastal area of the Eastern YS [26]. However, a reason for the large RMSD near the coastal shallow region has not yet been explained. The sea fog frequently forms in June and July and hinders measurement of SST by satellites IR-based sensors in the Eastern YS [14,39]. Sea fog forms frequently in June and July ( Figure 11). To investigate why fog forms frequently in June and July in the Southeastern YS, the monthly mean wind speed, relative humidity, and wind vector measured by the ocean buoys D and O were compared for the four years ( Figure 11). The monthly mean relative humidity in June and July were 88% and 91%, respectively. The wind speeds in June and July were 3.16 and 4.96 m/s, respectively. The sea fog forms frequently because of the temperature difference between the warm and humid air, and the relatively cold sea surface with some aid from the southerly wind in June and July. The ocean buoy D was located where the temperature difference between the SAT and SST was 2.1 • C. Since ocean buoy O was located at the boundary between the stratified and mixed regions, the temperature difference was 0.8 • C. In contrast, the ocean buoys KB and SB were located in the relatively deep stratified region and the temperature differences at KB and SB were −0.4 and −0.3 • C, respectively ( Figure 11).
Although the relative humidity, wind speed, and wind vector in August were similar to those in July, the number of foggy days (triangles and circles in Figure 11) drastically decreased in August. Zhang et al. [40] studied the reasons for the occurrence and disappearance of sea fog in the YS and East China Sea. The YS and East China Sea are situated between the subtropical low and mid-latitude high and therefore experience an easterly wind shift in the seasonal march from July to August. As the intensification of atmospheric convection over the subtropical Northwestern Pacific region excites a barotropic wave train from July to August, the meridional dipole pattern of a geopotential decrease over the subtropics and an increase in the mid-latitudes causes the prevailing winds to shift from southerly to easterly over the YS. This wind shift terminates the southerly wind and the northward advection of warm and moist air that sustains the YS fog from April to July. The atmospheric environment is therefore foggy in June and July, which obstructs the satellite-based observation of SST for the Eastern YS.
Remote Sens. 2020, 12, 759 18 of 25 Figure 11. Monthly evolution of differences (°C) between the surface air temperature and the SST measured by the ocean buoys D (cross), O (diamond), KB (square), and SB (closed dot). Temperature differences were calculated from the hourly data. The monthly variation of sea foggy days from March to October from Cho et al. [14] and Heo et al. [39]. The monthly variation of relative humidity (green star) and wind speed (black cross) are in the upper panel. Variations of monthly mean wind vector are in the lower panel.
Although the relative humidity, wind speed, and wind vector in August were similar to those in July, the number of foggy days (triangles and circles in Figure 11) drastically decreased in August. Zhang et al. [40] studied the reasons for the occurrence and disappearance of sea fog in the YS and East China Sea. The YS and East China Sea are situated between the subtropical low and mid-latitude high and therefore experience an easterly wind shift in the seasonal march from July to August. As the intensification of atmospheric convection over the subtropical Northwestern Pacific region excites a barotropic wave train from July to August, the meridional dipole pattern of a geopotential decrease over the subtropics and an increase in the mid-latitudes causes the prevailing winds to shift from southerly to easterly over the YS. This wind shift terminates the southerly wind and the northward advection of warm and moist air that sustains the YS fog from April to July. The atmospheric environment is therefore foggy in June and July, which obstructs the satellite-based observation of SST for the Eastern YS.

Availabily of a L3 SST in the YS
The AVHRR PFV53 L3 SST data are one of the most widely used data and they have five internal quality flags, such as best, acceptable, low, bad, and worst. The L4 gridded SST data were generated by the optimal interpolation of available L3 SST data with the best and acceptable quality [29]. For the Eastern YS, the number of AVHRR PFV53 L3 SST data with each quality flag was counted for June and July from 2010 to 2013 (Figure 12). In general, in June and July, the number of SST observations with the best and acceptable flags was relatively smaller in the shallow and vertically well-mixed regions. The SST data with the best and acceptable quality flags represented 8.2% of the Figure 11. Monthly evolution of differences ( • C) between the surface air temperature and the SST measured by the ocean buoys D (cross), O (diamond), KB (square), and SB (closed dot). Temperature differences were calculated from the hourly data. The monthly variation of sea foggy days from March to October from Cho et al. [14] and Heo et al. [39]. The monthly variation of relative humidity (green star) and wind speed (black cross) are in the upper panel. Variations of monthly mean wind vector are in the lower panel.

Availabily of a L3 SST in the YS
The AVHRR PFV53 L3 SST data are one of the most widely used data and they have five internal quality flags, such as best, acceptable, low, bad, and worst. The L4 gridded SST data were generated by the optimal interpolation of available L3 SST data with the best and acceptable quality [29]. For the Eastern YS, the number of AVHRR PFV53 L3 SST data with each quality flag was counted for June and July from 2010 to 2013 ( Figure 12). In general, in June and July, the number of SST observations with the best and acceptable flags was relatively smaller in the shallow and vertically well-mixed regions. The SST data with the best and acceptable quality flags represented 8.2% of the total observations in the tidally well-mixed region (Figure 12a,b) and 70.9% of the data had the worst quality flags.
Remote Sens. 2020, 12, 759 19 of 25 total observations in the tidally well-mixed region (Figure 12a, b) and 70.9% of the data had the worst quality flags. To evaluate the satellite-based L3 SST product in the shallow and vertically well-mixed regions, the RMSD and bias of the L3 SST data were calculated at ocean buoys D and O in June and July from 2010 to 2013 ( Figure 13). There were 17, 12, 48, 22, and 346 observational data available with the best, acceptable, low, bad, and worst quality flags, respectively. None of the satellite-based L3 SST data at buoy D were assigned as best quality data in June and July for the four years (Figure 13a). The RMSDs between the L3 SST data with the best, acceptable, low, bad, and worst quality flags and in situ SST data were found to be 0.60, 1.06, 2.37, 4.05, and 15.26 °C, respectively. The biases were found to be -0.38, 0.07, -0.22, -0.09, and -11.90 °C, respectively. Only the SST data classified as the best and acceptable quality were found to constitute relatively reliable data with an RMSD lower than 0.83 °C and the corresponding bias was smaller than -0.16 °C. To evaluate the satellite-based L3 SST product in the shallow and vertically well-mixed regions, the RMSD and bias of the L3 SST data were calculated at ocean buoys D and O in June and July from 2010 to 2013 ( Figure 13). There were 17, 12, 48, 22, and 346 observational data available with the best, acceptable, low, bad, and worst quality flags, respectively. None of the satellite-based L3 SST data at buoy D were assigned as best quality data in June and July for the four years (Figure 13a). The RMSDs between the L3 SST data with the best, acceptable, low, bad, and worst quality flags and in situ SST data were found to be 0.60, 1.06, 2.37, 4.05, and 15.26 • C, respectively. The biases were found to be −0.38, 0.07, −0.22, −0.09, and −11.90 • C, respectively. Only the SST data classified as the best and acceptable quality were found to constitute relatively reliable data with an RMSD lower than 0.83 • C and the corresponding bias was smaller than −0.16 • C. Some of the satellite-based L3 SST data plotted near the ideal regression line were classified as the worst quality data (Figure 13e). It is necessary to know how the quality flags of the AVHRR pathfinder L3 SST data are assigned to recover the data with incorrectly assigned flags. The AVHRR pathfinder L3 SST data processing performs a series of tests to determine whether a pixel contains SST data of suspicious quality including cloud contamination and uniformity tests [41]. The maximum and minimum brightness temperature values for channels 4 and 5 were calculated for 3 × 3 boxes centered around the pixel being classified. The difference between the maximum and minimum brightness temperatures for both channels must be less than 1.2 °C [31]. This uniformity test or SST gradient test seeks to identify contamination by small clouds and is based on the assumption that SSTs are relatively uniform at small scales (e.g., 3 × 3 pixels). However, pixels in the areas of sharp frontal features on scales smaller than 12 km might be erroneously identified as suspicious data by this test. The temperature gradient in the Eastern YS during the summer was calculated using the in situ SST data in June from 2010 to 2013 ( Figure 14). The horizontal temperature gradient was higher than 1.2 °C across 12 km along the tidal mixing front region in the Eastern YS Some of the satellite-based L3 SST data plotted near the ideal regression line were classified as the worst quality data (Figure 13e). It is necessary to know how the quality flags of the AVHRR pathfinder L3 SST data are assigned to recover the data with incorrectly assigned flags. The AVHRR pathfinder L3 SST data processing performs a series of tests to determine whether a pixel contains SST data of suspicious quality including cloud contamination and uniformity tests [41]. The maximum and minimum brightness temperature values for channels 4 and 5 were calculated for 3 × 3 boxes centered around the pixel being classified. The difference between the maximum and minimum brightness temperatures for both channels must be less than 1.2 • C [31]. This uniformity test or SST gradient test seeks to identify contamination by small clouds and is based on the assumption that SSTs are relatively uniform at small scales (e.g., 3 × 3 pixels). However, pixels in the areas of sharp frontal features on scales smaller than 12 km might be erroneously identified as suspicious data by this test. The temperature gradient in the Eastern YS during the summer was calculated using the in situ SST data in June from 2010 to 2013 ( Figure 14). The horizontal temperature gradient was higher than 1.2 • C across 12 km along the tidal mixing front region in the Eastern YS (Figures 10a and 14). This strong SST gradient could lead the uniformity test in the Pathfinder algorithm to classify good quality SST data in the frontal region as worst quality data.
Remote Sens. 2020, 12, 759 21 of 25 ( Figure 10a and Figure 14). This strong SST gradient could lead the uniformity test in the Pathfinder algorithm to classify good quality SST data in the frontal region as worst quality data. To investigate the relationship between sea fog formation and the L3 SST data quality flags, the L3 SST data flags were compared with the detected sea fog events from the MTSAT images at the locations of buoys D and O in June and July from 2010 to 2013. There were 445 in situ SST observations from the buoys and 488 values from the matching satellite-based L3 SST data. The number of buoy SST data was smaller than that of the L3 SST data because there were some days when the buoys could not measure SST. When the L3 SST data had the best and acceptable quality flags, no sea fog was found in the MTSAT images (Table 6). When sea fog formed in the lower column of the atmosphere, 37% of the L3 SST data had low or bad quality flags and 63% had the worst quality flags. This implied that sea fog formation degrades the quality of the satellite-based L3 SST data. Most (64%) of the low and bad quality data were related to sea fog and 36% of the low and bad quality data were not related to sea fog. However, most (78%) of the worst quality L3 SST data were not related to sea fog formation and only 22% of the L3 SST data with the worst quality flag were related to sea fog formation. This indicated that there were other causes, such as the strong horizontal SST gradient, which fails the uniformity test in the Pathfinder algorithm, which incorrectly assigned the worst quality flag to the good or acceptable quality L3 SST data (Figure 13e). To investigate the relationship between sea fog formation and the L3 SST data quality flags, the L3 SST data flags were compared with the detected sea fog events from the MTSAT images at the locations of buoys D and O in June and July from 2010 to 2013. There were 445 in situ SST observations from the buoys and 488 values from the matching satellite-based L3 SST data. The number of buoy SST data was smaller than that of the L3 SST data because there were some days when the buoys could not measure SST. When the L3 SST data had the best and acceptable quality flags, no sea fog was found in the MTSAT images (Table 6). When sea fog formed in the lower column of the atmosphere, 37% of the L3 SST data had low or bad quality flags and 63% had the worst quality flags. This implied that sea fog formation degrades the quality of the satellite-based L3 SST data. Most (64%) of the low and bad quality data were related to sea fog and 36% of the low and bad quality data were not related to sea fog. However, most (78%) of the worst quality L3 SST data were not related to sea fog formation and only 22% of the L3 SST data with the worst quality flag were related to sea fog formation. This indicated that there were other causes, such as the strong horizontal SST gradient, which fails the uniformity test in the Pathfinder algorithm, which incorrectly assigned the worst quality flag to the good or acceptable quality L3 SST data (Figure 13e). In the coastal region of the Southeastern YS, it was difficult to observe SSTs owing to poor atmospheric conditions in June and July ( Figure 11 and Table 6). Even if the SSTs are observed by satellite sensors in this coastal region, good or acceptable SST data may be classified as the worst quality data because of the strong horizontal temperature gradient higher than 1.2 • C/12 km across the tidal mixing front (Figure 14). While the L4 gridded SST data were produced using the OI method, there were insufficient good and acceptable L3 SST data in this coastal region. Then, the missing data points in the L4 gridded SST dataset were filled with a combination of SST climatology at the point and offshore SST data in warm water using the OI algorithm [29].

Background and Covariances Used for the Optimal Interpolation
Since the L3 SST data were sparse in the Eastern YS during summer as shown in the previous section, it was necessary to understand the process of optimal interpolation for the missing data for the L4 gridded SST dataset. As the OSTIA SST dataset was found to be the best among the five gridded SST datasets when compared to the in situ data, in this section, we used this dataset as an example of the L4 gridded SST data and examined its production procedure.
The L4 gridded SST is less reliable where the observation data are sparse because the missing SST at a grid point is interpolated from the data at the surrounding grid points. Especially in the coastal region of the YS, the tidal front develops along the coast and induces a strong horizontal SST gradient, causing more observational data errors (Figures 10 and 14). In June and July, most (93%) of the L3 SST data in the coastal region of the Southeastern YS were assigned as low, bad or worst quality data (Table 6), which resulted in the interpolation of SST data in the coastal region during the production of the L4 gridded SST datasets. The missing SST data in the coastal region were not interpolated from the nearby grid points but they were interpolated using the offshore observation data far away (10-70 km) in warm water and the climatology based on Equationa (3) and (4).
The L4 gridded SST datasets had the warm bias in the Eastern YS and the suggested two bias correction methods might be useful. However, it is not a fundamental solution. The background (x b i, k ) or climatology (x c i,k ) in Equation (3) needs to be improved using the long-term in situ observation data in this region. The background error covariance matrix B in Equation (4) has limited skill in the coastal regions, such as the Eastern YS, because the FOAM model does not include the tides, which are important forcing to produce the well-mixed cold water region in the Eastern YS. Near the tidal front in the Eastern YS, covariance B can be another error source during the interpolation process for the L4 gridded SST datasets.
Currently, the regularly obtained in situ SST data near the coast of YS is sparse. In the Southeastern YS where a strong SST gradient develops and satellite observations are technically difficult, the errors in the L4 gridded SST observation can be fundamentally reduced by installing more ocean buoys to measure SST.

Conclusions
The SST is an essential data for data assimilation of the ocean and atmospheric prediction systems. To evaluate SST datasets that cover the YS and to identify the error sources of the SST data, five L4 gridded SST datasets were selected and the spatio-temporal characteristics of the L4 gridded SST datasets were analyzed from January 2010 to December 2013. The SST from the OSTIA, GMPE and MGDSST datasets had a RMSD lower than 1 • C. The spatial distribution of RMSDs between the CTD SST and the five L4 gridded SST datasets revealed that the RMSD was higher in the vertically mixed coastal region within 50 km of the coastline or islands and was relatively lower in the offshore stratified region. In the mixed region, the RMSD was lowest for the OSTIA SST dataset and highest for the MWIR SST dataset. In the stratified region, the RMSDs of all L4 gridded SST datasets was less than 1 • C.
The RMSD of the L4 gridded SST data compared to in situ SST data was higher in June and July because of the warm bias in the shallow mixed region. To identify the cause of warm bias in the Southeastern YS, the satellite-based L3 SST data were examined. The L3 SST data were found to be rare, i.e., the number of high quality data were less than 30 in June and July for the four years at a grid point in the shallow mixed region. The L3 SST data were frequently missing because of sea fog formation in the coastal region and the large SST gradient in the region. The missing SST data were interpolated from the SST observation data in offshore warm water and warm biased SST climatology in the shallow coastal region. The optimal interpolation procedure for the L4 gridded SST was the main cause of the warm bias. To fundamentally improve the accuracy of the L4 gridded SST data, it is necessary to increase the number of SST observation data in the tidally well-mixed region. One way to improve the L4 gridded SST data in the Eastern YS is to loosen the horizontal gradient test threshold (1.2 • C/12 km) during the L3 data processing for June and July. Another way to achieve this would be to add more in situ observation data from ocean surface observation buoys in the tidal front zone and the shallow mixed region within 50 km from the land and islands. As a temporary solution to the warm bias in the gridded SST datasets in the Eastern YS, it is possible to improve the SST climatology for the optimal interpolation based on long-term field observation data.
Although the OSTIA SST dataset was the best in the YS, the RMSDs in the shallow mixed and stratified regions were 1.24 and 0.68 • C, respectively. To cut back the RMSD in the shallow mixed region, two bias correction methods were suggested. The bias correction reduced the RMSDs of the OISST, MGDSST, OSTIA, MWIR, and GMPE datasets by 41%, 43%, 37%, 47%, and 49%, respectively, in June.