*2.2. Methodology*

The GK2A LST retrieval algorithm is mainly composed of three parts, as shown in Figure 1. The first part is the development of the LST retrieval algorithm using a simulated pseudo match-up database with the RTM (shown in Figure 1 on the right side). First, we constructed the pseudo match-up database through the RTM simulation called MODTRAN4 under various atmospheric and land surface

conditions. In this simulation, to improve the accuracy of the LST retrieval algorithm, we considered various impacting factors, such as the diurnal variations of air temperature and LST, spatial-temporal variations of spectral emissivity, and the non-linear effect of water vapor. To develop LST retrieval algorithms, we used the reference LST and calculated the brightness temperatures of the two SW channels with the inverse Planck function from the simulated radiances. The LST retrieval algorithms were separately developed using a statistical regression method according to the lapse rate (two types: day/night) and atmospheric conditions (three types: wet/normal/dry).

**Figure 1.** Flow chart of GK2A's land surface temperature retrieval process.

The split-window LST retrieval algorithm assumes that the LSE is prior known. In general, LSE is affected not only by the surface type and status, but also depends on the wavelength. Because the spatial resolution of the GK2A/AMI infrared channel is 2 km, most pixels are composed of a mixture of various vegetation and soil. Moreover, most vegetation not only undergoes seasonal variations, but also, soil conditions are affected by the presence or absence of precipitation (soil moisture and snow cover). The process of retrieving the GK2A LSE is briefly shown on the left side of Figure 1. The GK2A LSE algorithm was developed based on the VCM method and calculates the emissivity of a given pixel as a weighted average according to the fractional coverage of vegetation (FVC) and soil [46,68]. The FVC of a given pixel was derived by the method of Carlson and Ripley (1997) [69] using the 8-days cycle vegetation index data derived from GK2A/AMI. The snow cover fraction was calculated for the pixels covered with snow [70], and this was considered in retrieving the daily GK2A LSE.

The central column of Figure 1 shows the process of calculating LST every 10 min using GK2A/AMI data and pre-calculated emissivity. The GK2A LST was retrieved for clear and land pixels according to day/night and atmospheric conditions. In the process of calculating the GK2A LST, the SZA of each pixel was used for the division of daytime, nighttime, and dawn/twilight, as shown in the center of Figure 1. At this time, the threshold SZA values from the previous study were used [60]. Similarly, the atmospheric conditions were divided into dry/normal/wet according to the brightness temperature difference (BTD) threshold to calculate the GK2A LST.

Unlike the sea surface temperature, the available match-up data (ground-level observational data) for LST over the GK2A/AMI observation region are severely limited [71–73]. Therefore, to develop the GK2A LST algorithm from satellite data, we needed to prepare a similar on-site match-up database. For this purpose, the output from the RTM simulation could be used. As in many previous studies, the pseudo match-up database was generated through simulations using the RTM under various atmospheric and surface conditions. We developed an LST retrieval algorithm using the pseudo match-up database. SeeBor version 5 data have been used as an input profile for RTM simulations in many studies for the generation of a pseudo match-up database. This dataset contains 101 pressure levels of temperature, mixing ratio, and trace gases. For each profile in the dataset, a physically-based characterization of the surface skin temperature and surface emissivity is assigned [74]. The diurnal variation of LST is not provided in atmospheric profiles, so many studies have considered the diverse diurnal variation of LST in the RTM simulations according to the land cover, season, and weather conditions. In the research of [26,34], a large range of temperatures between the near-surface temperature and ground surface, which consists of five surface temperatures, Ta − 5 K, Ta, Ta + 5 K, Ta + 10 K, and Ta + 20 K, was considered for each atmospheric profile. In the research of [27,75], daily variations of LST from Ta − 16 K to Ta + 16 K were considered for each profile. In other studies [29,76], LST was prescribed for each profile in a range as Ta − 15 K < LST < Ta + 15 K with different increments of 1.5 K and 1 K, respectively. In [35], LST was considered asymmetrically from Ta − 6 K to Ta + 14 K in increments of 2 K. In [55], they used six values from Ta − 5 K to Ta + 20 K at 5 K intervals. In this study, the diurnal variation of LST was asymmetric according to the day (LSTday = Ta − 2 K ∼ Ta + 18 K) and night (LSTngt = Ta − 6 K ∼ Ta + 2 K) at 2 K intervals. In addition, we included the diurnal variations of air temperature profiles under the planetary boundary layer using SeeBor data and reference LST, as shown in Figure 2.

**Figure 2.** Conceptual atmospheric vertical profiles according to diurnal variation of land surface temperature (LST) and air temperature (bold red line: SeeBor profile; other lines: diurnal variation of air temperature profiles).

To generate the pseudo match-up database of the LST and measured spectral radiance, radiative transfer model simulations under various atmospheric and surface conditions were performed (Table 2). A total of 2694 sets of SeeBor data used were located at a satellite with a VZA of less than 50◦ . To consider the value of LSE in the RTM simulation, the range of LSE of channels 13 and 15 in the GK2A region was used, as shown in Table 2. Therefore, the total number of simulations for daytime and nighttime were 3,585,714 (2694 (atm) × 13 (LST) × 11 (ε*ch*13) × 11 (∆ε)) and 1,629,870 (2694 (atm) × 5 (LST) × 11 (ε*ch*13) × 11 (∆ε)), respectively.


**Table 2.** The input conditions of the radiative transfer model simulation according to the various atmospheric and surface conditions.

Several LST retrieval algorithms from satellites have been proposed to suit the characteristics of various sensors mounted on each satellite and retrieve LSTs through different approximations and assumptions. Among the various LST retrieval methods, the SW method was used to retrieve LST using different atmospheric absorptions of two adjacent infrared window channels in this study. The SW method assumes that the spectral emissivity of the land surface is a priori known, and the atmospheric effects can be corrected using thermal infrared channels. The LST was estimated through a regression analysis on the relationship between the brightness temperature of two infrared channels and the factors affecting the LST calculation [33,34,36,60,77]. The equation for retrieving LST from GK2A data is given in Equation (1) as follows:

$$LST = c\_0 + c\_1 BT\_{\text{cl13}} + c\_2 (BT\_{\text{cl15}} - BT\_{\text{cl15}}) + c\_3 (BT\_{\text{cl13}} - BT\_{\text{cl15}})^2 + c\_4 (\sec \theta - 1) + c\_5 (1 - \overline{\varepsilon}) + c\_6 \Delta \varepsilon \tag{1}$$

where *BTch*<sup>13</sup> and *BTch*<sup>15</sup> are the brightness temperatures of GK2A channels 13 and 15, respectively, θ is the satellite VZA, ε is the average LSE of GK2A channels 13 and 15, and ∆ε is the LSE difference of GK2A channels 13 and 15. The coefficients of the regression equations (*c*0~*c*6) were derived by simple regression.

As described in Section 2.1.2, MODIS Collection 6 LST data and BSRN Tateno station data were used for the validation of GK2A/AMI LST. The spatial resolution of the MODIS LST was 1 km. The MODIS LST data were provided at 5 min intervals as Terra/Aqua satellites moved along their orbits. Due to the characteristics of LST with large spatial-temporal variability, strict spatial-temporal collocations are needed. For temporal collocation, the MODIS LST data within ±5 min were used based on the observation time of GK2A, as in previous studies [60]. For the spatial collocation, a simple mean of the nearest 9 MODIS pixels (an average of 3 × 3) surrounding the closest MODIS pixel to the GK2A pixel was used. Validation with MODIS LST data was performed when the pixels had more than five pixels of clear sky and land.

In addition, the in situ measurement data, which were observed every minute from the upward longwave radiation at the Tateno station of the BSRN, were used as the validation data. In the GK2A observation area, there were few upward longwave radiation data or field-observed LST data, so Tateno station data were used as validation data. Validation was performed only when the observation times at GK2A and Tateno station were the same. Spatial collocation was performed by averaging the four GK2A pixels closest to the Tateno station. The measured upward longwave radiation was converted to LST using the Stefan–Boltzmann law. The surface type of Tateno station is grass, so we assumed the average value of the emissivity (ε = 0.986) value corresponding to the grassland from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) spectral library. The process of converting the observed upward longwave radiation to LST using the Stefan–Boltzmann law is shown in Equation (2):

$$\begin{aligned} L\_{\text{sensor}} &= \overline{\varepsilon} \sigma T^4\\ T &= \left( \frac{L\_{\text{escor}}}{\overline{\varepsilon} \sigma} \right)^{\frac{1}{4}} \end{aligned} \tag{2}$$

where *Lsensor* is the value of the measured upward longwave radiation by the CGR4 pyrgeometer from the Tateno station of the BSRN; ε is the average value of the emissivity corresponding to the grassland; σ is the Stefan-Boltzmann constant (σ = 5.670374 × 10−<sup>8</sup> Wm−2K −4 ); *T* is the value of ground temperature converted from the upward longwave radiation at Tateno station.

#### **3. Results**

#### *3.1. Results of the Radiative Transfer Model Simulation*

To set the threshold value of the LST algorithm, the distribution of brightness temperature differences (BTD) between channel 13 and channel 15 of GK2A was analyzed. This time, only the land pixels without clouds were analyzed. The analysis dates were selected as the same day each month from July 2019 to October 2019, and the count of BTDs accumulated was observed every 10 min. Figure 3 shows the frequency of BTDs between channel 13 and channel 15 of GK2A. BTD values were distributed from −3 K to 15 K. Most of the BTDs were counted at over 10 million from 1 K to 7 K. The results showed a similar distribution to that of the RTM simulated results, which are shown in Figure 4. Figure 4 shows the BTD distribution between channel 13 and channel 15 of GK2A from the simulated pseudo match-up database with the RTM. The coefficients of the LST equation (Equation (1)) were obtained from the simulated pseudo match-up database by regression analysis. Using these coefficients, the root mean square error (RMSE) values between the reference LST and the estimated LST were calculated. Overall, the distribution of simulated BTDs was similar to that of the actual GK2A BTDs, except for the −3 K to 0 K range. The RMSE value was relatively large when the simulated BTD value was less than −1 K or greater than 6 K. The difference between the actual GK2A BTD and the RTM simulation BTD was less than 0 K. This difference could have been caused by the small amount of aerosols given the atmospheric profile (SeeBor\_v5), which was used for the RTM simulation [74]. Using these results, the threshold values of BTD for dry/normal/wet conditions were determined to be 0 K and 6 K, respectively.

**Figure 3.** Distribution of brightness temperature difference between channel 13 and channel 15 of GK2A for selected days (25 July 2019; 25 August 2019; 22 September 2019; and 4 October 2019).

**Figure 4.** Distribution of radiative transfer model simulated brightness temperature differences between channel 13 and channel 15 of GK2A (green histogram) and their root mean square error (RMSE) values (blue line).

The coefficients of the GK2A/AMI LST equation (Equation (1)) for day/night and dry/normal/wet conditions were obtained from the RTM simulation results through regression analysis, as shown in Table 3.


**Table 3.** The coefficients of the GK2A/AMI equation (Equation (1)) according to the day/night and dry/normal/wet conditions.

To evaluate the RTM simulation results and the GK2A LST retrieval algorithm, we analyzed the relationship between the reference LST value and estimated LST values, as shown in Figure 5.

**Figure 5.** (**a**) Scatter plot (left) and (**b**) bias distribution (right) between the reference LST and estimated LST values using the GK2A/AMI LST retrieval algorithm (gray dotted line in (**a**) represents the 1:1 line).

The estimated LST generally matched well with the reference LST over a wide range from 245 K to 330 K. The comparison of the retrieved LST with the reference LST showed that the LST retrieval algorithm had a high retrieval accuracy in terms of correlation (0.998), bias (0.010 K), and RMSE (0.767 K). The LST retrieval algorithm, most of the reference LSTs, and the estimated LSTs were distributed based on a 1:1 line. In addition, most of the bias was less than ±3 K, with a peak value at 0 K. However, when the LST was greater than 300 K, the LST retrieval algorithm slightly over/underestimated the LSTs compared to the reference LST.

Figure 6 shows the distribution of RMSEs between the estimated LSTs and reference LSTs based on the impact factors (lapse rate, BTD, surface emissivity differences, and satellite VZAs that affect the retrieval accuracy of the LST). Most of the RMSEs are less than 2.1 K, excluding some ranges where the VZA is larger than 40◦ and the BTD value is larger than 7 K. The reason that the RMSE was large when the VZA was large appears to be due to the fact that the effect of the emissivity becomes significant when the VZA is large, as suggested in many studies [78,79]. One of the most distinctive features of all the cases is that the RMSE was slightly larger when the BTD values were larger than 6 K. The magnitude of BTDs is mainly caused by the different sensitivities of the two window channels to

aerosol (channel 13) and water vapor (channel 15). In addition, the GK2A/AMI LST retrieval algorithm is significantly influenced by the VZA and the surface lapse rate among several impacting factors. Compared to previous studies that retrieved LST from Himawari-8, RMSE decreased in all ranges [60]. In particular, when BTD was a positive large value in the previous study, most of the RMSEs were significantly reduced. These results seem to be related to the consideration of the diurnal variation of the boundary layer temperature in the RTM simulations and the consideration of the effect of nonlinearity of water vapor in the LST calculation formula.

**Figure 6.** Distribution of RMSEs for the LST retrieval algorithms based on the different impacting factors: (**a**) brightness temperature difference (BTD) and surface lapse rate; (**b**) emissivity difference (∆ε) and BTD; (**c**) ∆ε and surface lapse rate; (**d**) satellite viewing zenith angle (VZA) and BTD; (**e**) VZA and surface lapse rate.

#### *3.2. Comparison of Retrieved GK2A LST and MODIS LST Products*

To evaluate the GK2A/AMI LST retrieval algorithm, GK2A LSTs were retrieved for 4 months, from July 2019 to October 2019, when GK2A started operational observation. Unlike the RTM simulations, when retrieving LSTs from GK2A/AMI data, the criterion for dividing day and night was the SZA of each pixel, which is the same as that in previous studies [36,60]. Only clear sky and land pixels that satisfied the strict spatial-temporal matching conditions were compared and validated, as mentioned in Section 2.1.2. Figures 7–10 show the spatial distribution of the LSTs retrieved from GK2A/AMI data and the MODIS LSTs, along with the spatial distribution of their differences between two datasets for the selected days.

**Figure 7.** Spatial distribution of (**a**) GK2A LST, (**b**) MODIS (MYD11\_L2) LST, (**c**) differences between GK2A and MODIS LSTs, and (**d**) scatter plot between GK2A and MODIS (MYD11\_L2) LSTs for 30 August 2019, 1700 and 1705 UTC (gray dotted line in (**d**) represents the 1:1 line).

**Figure 8.** Spatial distribution of (**a**) GK2A LST, (**b**) MODIS (MOD11\_L2) LST, (**c**) differences between the GK2A and MODIS LSTs, and (**d**) scatter plot between the GK2A and MODIS (MOD11\_L2) LSTs for 30 August 2019, 0030 UTC (gray dotted line in (**d**) represents the 1:1 line).

**Figure 9.** Spatial distribution of (**a**) GK2A LST, (**b**) MODIS (MOD11\_L2) LST, (**c**) differences between GK2A and MODIS LSTs, and (**d**) scatter plot between GK2A and MODIS (MOD11\_L2) LSTs for 21 September 2019, 1330 and 1335 UTC (gray dotted line in (**d**) represents the 1:1 line).

**Figure 10.** Spatial distribution of (**a**) GK2A LST, (**b**) MODIS (MYD11\_L2) LST, (**c**) differences between the GK2A and MODIS LSTs, and (**d**) scatter plot between the GK2A and MODIS (MYD11\_L2) LSTs for 20 October 2019, 0540 UTC (gray dotted line in (**d**) represents the 1:1 line).

Figure 7 shows the spatial distribution of the GK2A LST and MODIS LST and their differences at 1700 and 1705 UTC on 30 August 2019. The two LSTs show similar spatial distributions in large areas of the Korean Peninsula and northeast China. The difference between the two LSTs was within ±2 K. In the scatter plot of the two LSTs, the distribution spread slightly to both sides around the 1:1 line from 270 K to 305 K. The bias and RMSE of this scene are −0.113 K and 1.05 K, respectively.

Figure 8 shows the spatial distribution of the two LSTs at 0030 UTC on 30 August 2019. The area in which the LST was retrieved is a region between central and southern Australia consisting of desert, bare soil, and grassland. The difference between the two LSTs in the spatial distribution showed, where the GK2A LST was greater than the MODIS LST. In particular, the warm bias of the GK2A LST was more significant in the desert regions of south Australia than in other regions. For this region, even in the scatter plot, the GK2A LST was higher than the MODIS LST over a wide range (275–310 K). In this case, the bias and RMSE are 0.81 K and 1.264 K, respectively.

The spatial distribution of the two LSTs and their differences, along with the scatter plot on 21 September 2019, are shown in Figure 9. The observation time of the two LSTs differs by 5 min. The spatial distribution of the two LSTs is generally similar, but the spatial distribution of the difference between the two LSTs showed a warm bias in the desert region of Australia's inland but a cold bias in south Australia. As shown in the scatter plot, the overall GK2A LST was higher than the MODIS LST in all ranges (270–305 K). As a result, the GK2A LST shows a warm bias of 0.924 K and an RMSE of 1.554 K.

The spatial distribution of the GK2A LST, the MODIS LST, and their differences at 0540 UTC on 20 October 2019, are shown in Figure 10. The two LSTs have similar spatial distributions in large areas of Mongolia and China. The differences between the two LSTs were within ±3 K. In the scatter plot of the two LSTs, the distribution showed a 1:1 line from 270 K to 295 K, but some pixels were spread around the 1:1 line. The bias and RMSE of this scene are 0.137 K and 1.411 K, respectively.

The indirect validation results of the GK2A LST with the MODIS LST for 4 months (July, August, September, and October 2019) are shown in Table 4. The correlation coefficient between GK2A LST and MODIS LST was greater than 0.96 regardless of the satellite (Terra/Aqua) and months. In addition, the bias showed a warm bias of less than 1.6 K every month. The combined results of the two validation satellites show a high correlation coefficient (0.969), slightly warm bias (1.227 K), and RMSE (2.281 K).

Because the diurnal variation pattern of LST and the retrieval algorithm are very different according to day and night, the performance of the GK2A LST algorithm was evaluated accordingly. Table 5 shows the comparison results between the GK2A and MODIS LSTs according to the daytime and nighttime. The performance of the GK2A LST algorithm is dependent on the validation time, day, and night. In general, the skill of the GK2A LST retrieval algorithm for daytime was about two times worse than that of nighttime, in terms of bias and RMSE, irrespective of the months. The relatively strong warm biases of the GK2A LST during daytime can be attributed to the current status of the MODIS LST [80,81]. The biggest feature is that the GK2A LST overestimated more than the MODIS LST during the daytime, showing a large warm bias (above 1.8 K) and RMSE (above 2.8 K). Moreover, for nighttime, the RMSE was small due to a relatively high correlation and a small bias compared to the daytime.


**Table 4.** Comparison results of the GK2A LST data and the MODIS LST product (Collection 6) for the selected four months (July, August, September, and October 2019).
