Advanced Dual-Satellite Method for Detection of Low Stratus and Fog near Japan at Dawn from FY-4A and Himawari-8

: The detection of low stratus and fog (LSF) at dawn remains limited because of their optical features and weak solar radiation. LSF could be better identiﬁed by simultaneous observations of two geostationary satellites from different viewing angles. The present study developed an advanced dual-satellite method (DSM) using FY-4A and Himawari-8 for LSF detection at dawn in terms of probability indices. Optimal thresholds for identifying the LSF from the spectral tests in DSM were determined by the comparison with ground observations of fog and clear sky in/around Japan between April to November of 2018. Then the validation of these thresholds was carried out for the same months of 2019. The DSM essentially used two traditional single-satellite tests for daytime such as the 0.65- µ m reﬂectance (R 0.65 ), and the brightness temperature difference between 3.7 µ m and 11 µ m (BTD 3.7-11 ); in addition to four more tests such as Himawari-8 R 0.65 and BTD 13.5-8.5 , the dual-satellite stereoscopic difference in BTD 3.7-11 ( ∆ BTD 3.7-11 ), and that in the Normalized Difference Snow Index ( ∆ NDSI). The four were found to show very high skill scores (POD: 0.82 ± 0.04; FAR, 0.10 ± 0.04). The radiative transfer simulation supported optical characteristics of LSF in observations. The LSF probability indices (average POD: 0.83, FAR: 0.10) were constructed by a statistical combination of the four to derive the ﬁve-class probability values of LSF occurrence in a grid. The indices provided more details and useful results in LSF spatial distribution, compared to the single satellite observations (i.e., R 0.65 and/or BTD 3.7-11 ) of either LSF or no LSF. The present DSM could apply for remote sensing of environmental phenomena if the stereoscopic viewing angle between two satellites is appropriate.


Introduction
Reliable monitoring of low stratus and fog (LSF) is of vital importance because restricted visibility can cause hazards in transportation and navigation including LSF-related accidents (e.g., [1,2]). Thick LSF at dawn is a major hazard during the morning rush hours. The persistence and frequency of LSF on land can be associated with various factors: air pollution and quality [3][4][5] and climate variabilities such as Arctic Oscillation [6] and El Niño/Southern Oscillation [7,8].
According to Yoo et al. [9], 71-76% of fog phenomena in surface synoptic observations (SYNOP) tend to accompany the low stratus over the Korean Peninsula at summertime dawn. This indicates that the better the LSF detection, the higher the possibility for fog detection. Low stratus and fog are often decoupled from the upper part of fog at night. Thus, low stratus and fog are examined in terms of LSF without individual distinctions in this study. LSF is also an appropriate term in satellite remote sensing since the distinction between low stratus and fog is fundamentally difficult without additional ground information such as surface elevation and cloud base height, etc. [2,3,[10][11][12]. The only difference between the two weather phenomena is whether their bases are in contact with the ground [13].  Table (LUT) in advance, were comprehensively used for an approximate diagnostic check of the LSF optical properties. In Section 2, the data and method were described. Here spectral properties of newly developed tests were explained. The results were shown in terms of the LSF and clear-sky classification using optimal thresholds with radiative transfer simulations in Section 3. The case study for validation was shown with the PI formulation. Discussion and conclusions were presented in Sections 4 and 5, respectively.

Materials and Methods
To derive and validate the threshold values for LSF detection at dawn near Japan, the following three datasets were used: geostationary satellite observations from Himawari-8 and FY-4A, ground observations from Meteorological Aerodrome Report (METAR) for the weather phenomena of fog and clear sky, and the LUT from RTM. For this study, "dawn" was defined as within 2 h after sunrise (67° < SZA < 86°). The detailed information for the datasets is described below. For SZA, the average of LSF cases in 2018 was ~79° and that of clear-sky conditions was ~78°. In 2019, both SZA values were ~80°. Also, the meanings of acronyms used in this study are noted in Table A1.

Satellite Observations
The geostationary satellites of Himawari-8 and FY-4A are located at different longitudes over the equator, and they show a VZA difference of 40.4˚ near the nadir of Himawari-8 ( Figure 1). The satellite information and application in this study with respect to the channel are described in Table 1. The primary payloads of FY-4A and Himawari-8 are the Advanced Geostationary Radiation Imager (AGRI; [38,39]) and AHI; [40], respectively. The instruments provided significant advancements in terms of the number of bands, spatial resolution, and temporal frequency. Their solar and thermal channels are separately presented in the table. There were two visible (VIS) channels of solar radiation  Table (LUT) in advance, were comprehensively used for an approximate diagnostic check of the LSF optical properties. In Section 2, the data and method were described. Here spectral properties of newly developed tests were explained. The results were shown in terms of the LSF and clear-sky classification using optimal thresholds with radiative transfer simulations in Section 3. The case study for validation was shown with the PI formulation. Discussion and conclusions were presented in Sections 4 and 5, respectively.

Materials and Methods
To derive and validate the threshold values for LSF detection at dawn near Japan, the following three datasets were used: geostationary satellite observations from Himawari-8 and FY-4A, ground observations from Meteorological Aerodrome Report (METAR) for the weather phenomena of fog and clear sky, and the LUT from RTM. For this study, "dawn" was defined as within 2 h after sunrise (67 • < SZA < 86 • ). The detailed information for the datasets is described below. For SZA, the average of LSF cases in 2018 was~79 • and that of clear-sky conditions was~78 • . In 2019, both SZA values were~80 • . Also, the meanings of acronyms used in this study are noted in Table A1.

Satellite Observations
The geostationary satellites of Himawari-8 and FY-4A are located at different longitudes over the equator, and they show a VZA difference of 40.4 • near the nadir of Himawari-8 ( Figure 1). The satellite information and application in this study with respect to the channel are described in Table 1. The primary payloads of FY-4A and Himawari-8 are the Advanced Geostationary Radiation Imager (AGRI; [38,39]) and AHI; [40], respectively. The instruments provided significant advancements in terms of the number of bands, spatial resolution, and temporal frequency. Their solar and thermal channels are separately presented in the table. There were two visible (VIS) channels of solar radiation centered at 0.65 µm and~1.6 µm. In the 0.65 µm channel, clouds appeared bright due to their high reflectance, while the signals from the land/sea surface were generally weak. The 1.6 µm channel was in the infrared (IR) region and allowed distinction between snow/ice-covered land and signals from other sources. This channel was sensitive to the cloud phase, and water clouds reflected more energy in this spectral range than the ice clouds. Of the six channels used, the four non-solar channels were classified as thermal parts. The 3.7 µm lies in the spectrum of two outgoing energy sources (e.g., solar reflectance and thermally emitted radiation). This channel is not only sensitive to the cloud phase but also useful for the detection of water clouds at night. In the region around 8.5 µm, there is a little water vapor absorption, but it is still in the window region. The IR window channel around 11 µm is affected primarily by radiation from the Earth's surface or cloud tops. The 13.3-13.5 µm band is used as cloud masks and cloud-top height.
Level 1B data as digital count values were used to convert the VIS and near-infrared (NIR) data to visible reflectance and the others (mid-IR, longwave-IR1, longwave-IR2, longwave-IR3) to brightness temperature. In order to observe the two satellite observations on the same-sized grid, collocation was performed by interpolating the Himawari-8 observations to FY-4A along with the ground observations shown in Table 2 and Figure 2. The analysis of this study is not so sensitive to the grid-size of either 4 km by 4 km or the original pixel (2 km by 2 km), but the data smoothing effect tends to increase in the former grid. To compare the dual-satellite observations by pixel, the spatial resolution at the channels noted above was set to 4 km by 4 km. Thus, the Himawari-8 data in a two-dimensional array of longitude and latitude near Japan (122.5-144E, 23.5-44N) were interpolated into the FY-4A array. This method keeps the FY-4A accuracy without degradation. The observations for seven variables (BTD 13.5-8.5 , R 0.65 , BTD 3.7-11 , NDSI, ∆R 0.65 , ∆NDSI, and ∆BTD 3.7-11 ) were analyzed for LSF detection. Additionally, the difference value between Land Surface Temperature (LST) and 11 µm brightness temperature (BT 11 ) was used to analyze cloud-top height as an acronym of LST-BT 11 . In this study, BTD 13.5-8.5 was the brightness temperature difference between 13.5 µm and 8.5 µm. Based on the main weighting function peak in the Satellite Infrared Spectrometer (SIRS-B) channel at 750 cm −1 [42], the peak of AHI at 13.28 µm (753 cm −1 ) was located near the 850-900 hPa level. Meanwhile, the 8.5 µm AHI surface channel was used for the analysis of the cloud phase. Reflectance values (R 0.6 and R 1.6 ) were normalized by dividing by cos(SZA) and set to unity when the values were greater than 1 due to normalization (e.g., large SZA ≥~80 • ). The difference in values between the two collocated satellites is defined with Equations (1)-(3) as follows: Four variables (BTD 13.5-8.5 , R 0.65 , BTD 3.7-11 , and NDSI) were obtained from a single satellite's (i.e., Himawari-8 only) observations, while three variables (∆R 0.65 , ∆NDSI, and ∆BTD 3.7-11 ) were derived from both satellites. Also, LST data were given from ground observations.
The NDSI index has been widely used to distinguish between snow and clouds since snow has a lower reflectance at 1.6 µm [43]. The studies [14,44] presented the NDSI values for various surface conditions (e.g., water clouds, ice clouds, and vegetation,) to eliminate snow pixels, which have relatively large R 0.65 values. Wu and Li [45] determined both the NDSI and R 0.67 thresholds for daytime sea fog detection using MODIS data. In their study, the NDSI was used for LSF detection over inland China to divide vegetation and clouds during the warm season, because vegetation has larger R 1.6 and smaller R 0.65 than clouds. Ryu and Hong [46] proposed a regression-based method for sea fog detection from a combination of NDSI and reflectance in the AHI green band (0.51 µm).
The BTD13.5-8.5 value which is the difference in brightness temperature between 13.5 and 8.5 μm was used to detect higher ice clouds [47,48], because the 13.5 μm channel belongs to a CO2 absorption band that affects atmospheric transmission for the calculation of cloud altitudes. Thus, this value is sensitive to atmospheric altitude. This result indicates that the BTD13.5-8.5 value can be utilized as an indicator to classify LSF among the weather phenomena of clear sky, LSF, and high clouds.

Ground Observations
In this study, ground observation data from METAR were regarded as the ground truth for LSF and clear-sky cases from March to November in 2018 and 2019. The METAR is a weather status report designed for aviation based on information from an Automatic Weather Station (AWS), such as visibility, cloud amount, cloud base height present weather, ambient temperature, and relative humidity at hourly intervals. Figure 2 shows the locations of METAR stations in Japan including the "Ulleung" island station in South Korea. The information for a total of 31 stations over the study region is also described in Table 2. The LSF occurrences were selected only in cases where the current weather conditions were reported as "FG" (i.e., fog). In METAR data, fog appears in various forms depending on accompanying weather phenomena and full or partial spatial covering. Only cases of "FG" that covered the entire observation area without precipitation were selected in this study. The clear-sky cases as a contrast group of LSF cases were adopted under conditions with the abbreviation "CAVOK" or cloud amount with "FEW" (1 to 2 oktas). CAVOK stands for Ceiling (or clouds) And Visibility OK, indicating no clouds below 1500 m or the highest minimum sector altitude and no cumulonimbus or towering cumulus at any level, with a visibility of 10 km or more and no significant weather change [49]. The total observation numbers for LSF and clear-sky occurrences used in this study were 347 and 229 during 2018 and 2019, respectively.
Weather maps (i.e., SYNOP chart) provided by the Korea Meteorological Administration (KMA) were used to verify the presence of LSF and upper cloud information detected by satellites and to obtain extensive information on air pressure The BTD 13.5-8.5 value which is the difference in brightness temperature between 13.5 and 8.5 µm was used to detect higher ice clouds [47,48], because the 13.5 µm channel belongs to a CO 2 absorption band that affects atmospheric transmission for the calculation of cloud altitudes. Thus, this value is sensitive to atmospheric altitude. This result indicates that the BTD 13.5-8.5 value can be utilized as an indicator to classify LSF among the weather phenomena of clear sky, LSF, and high clouds.

Ground Observations
In this study, ground observation data from METAR were regarded as the ground truth for LSF and clear-sky cases from March to November in 2018 and 2019. The METAR is a weather status report designed for aviation based on information from an Automatic Weather Station (AWS), such as visibility, cloud amount, cloud base height present weather, ambient temperature, and relative humidity at hourly intervals. Figure 2 shows the locations of METAR stations in Japan including the "Ulleung" island station in South Korea. The information for a total of 31 stations over the study region is also described in Table 2. The LSF occurrences were selected only in cases where the current weather conditions were reported as "FG" (i.e., fog). In METAR data, fog appears in various forms depending on accompanying weather phenomena and full or partial spatial covering. Only cases of "FG" that covered the entire observation area without precipitation were selected in this study. The clear-sky cases as a contrast group of LSF cases were adopted under conditions with the abbreviation "CAVOK" or cloud amount with "FEW" (1 to 2 oktas). CAVOK stands for Ceiling (or clouds) And Visibility OK, indicating no clouds below 1500 m or the highest minimum sector altitude and no cumulonimbus or towering cumulus at any level, with a visibility of 10 km or more and no significant weather change [49]. The total observation numbers for LSF and clear-sky occurrences used in this study were 347 and 229 during 2018 and 2019, respectively.
Weather maps (i.e., SYNOP chart) provided by the Korea Meteorological Administration (KMA) were used to verify the presence of LSF and upper cloud information detected by satellites and to obtain extensive information on air pressure patterns and wet areas. However, since the weather maps were produced at three-hour intervals, a direct comparative analysis was not possible due to the ordinary rapid generation and dissipation of fog. Therefore, weather map analysis was performed only when the map time coincided with the ground or satellite observation time. If there was a time difference between fog occurrence time and weather map, the BT 11 distribution, which is approximately related to cloud top temperature, was used instead.

Radiative Transfer Model
Optical characteristics of LSF can be estimated using the RTM simulation with satellite observations (e.g., [50]). In this study, the Santo Barbara DISORT (SBDART) of a planeparallel radiative transfer model was used [51] and the LUT was calculated from RTM in advance for estimation. The cloud model in the SBDART yielded a radiance value relating to cloud optical properties (e.g., Fog Optical Thickness, FOT; Cloud Optical Thickness, COT; Effective Radius, ER) based on the Mie theory. This allowed us to understand satelliteobserved spectral values in the LSF environment and their optical difference between the two satellites. In addition, the simulation provided insight into the effects of higher clouds over the LSF layer.
The Spectral Response Function (SRF) difference between the two satellites was less significant than that of the VZA [9]. Therefore, the values at the central wavelength of each channel (https://www.data.jma.go.jp accessed on 3 June 2019 for Himawari-8; http://satellite.nsmc.org.cn accessed on 5 June 2019 for FY-4A) were used instead of the SRF integrated over its entire wavelength range, of which the calculation process was time-consuming. Table A2 shows the input variables used for the RTM run in this study with their acronyms. The LUT was set up for five conditions similar to those in the study by [36] (their Figure A1). LSF1 at 0-1 km and 0-2 km without higher clouds above the LSF layer shown in Cases A and B. Cases C to E were with LSF2 accompanying the clouds. Case C and D were the clouds at 4-6 km above the 2 km fog layer, composed of either water particles or ice, respectively. Also, in Case E, clouds of ice particles exist 8-10 km above LSF1. We obtained simulated values from the LUT based on the input data of three sun or satellite orbit angles (SZA, VZA, and RAA) for each LSF case. Figure 3 shows the flow diagram for LSF detection near Japan at dawn during the period from March to November of 2018-2019, based on the near-simultaneous satellite observations of Himawari-8 and FY-4A. After collocating 31 ground station sites onto satellite observation grids, satellite data were collected when fog occurred at the ground stations. The optimal thresholds of eight satellite-derived variables were derived from iteration processes to discriminate between LSF and clear-sky cases based on a skill-score test. The variables were BTD 13.5-8.5 , ∆NDSI, R 0.65 , ∆BTD 3.7-11 , LST-BT 11 , ∆R 0.65 , BTD 3.7-11 , and NDSI. Using the top four skill-score variables, a total of 16 cases could be extracted from their statistical combinations (Tables 3 and 4). In addition, the weights were given and merged into a total of five classes for the probability of LSF occurrence for the 16 cases. To calculate PI in real time, calculation of five PI classes and the AP (i.e., assigned probability) value in each class were conducted using the long-term data of the other training period (i.e., the 2018 control period in this current study). The introduction of the AP values is necessary to quantify the LSF occurrence by its probability. The magnitude order of AP was set based on the top four "POD minus FAR" values in the skill score test of the pre-processing. "POD minus FAR" could comprehensively explain both POD and FAR. This "PI formulation" is explained in detail in Section 3.3.  Table A1.   Table A1. Utilizing these climatological values during the training period, a threshold test of the real-time dual satellite observations was carried out and provided the PI class (or no fog) in a spatial grid. The AP as a function of PI class (i.e., Class obs ) was equivalent to the PI obs value, indicating the LSF probability from the real-time dual satellite observations. Fog detection methods from simultaneous observations of multiple satellites have been developed in two previous studies [9,36]. We analyzed various channel information obtained from advanced satellites and tried to present the probability index of fog occurrence more precisely in this study compared to previous studies. Table 5 shows a summary of the DSM differences between previous and the present studies. Table 4. Five LSF classes derived from 16 cases of the four variables (A: BTD 13.5-8.5 , B: ∆NDSI, C: R 0.65 , D: ∆BTD 3.7-11 ) based on the threshold test for detection of either LSF or clear-sky phenomena. For instance, Case 1, which passed the test for LSF detection by all four variables, corresponds to the "very high" probability of Class 1 in terms of the probability set "P(A ∩ B ∩ C ∩ D)". The superscript "c" in A C , which consists of the elements that are not in A in set theory, means the complement of a set A.  The maximum number of PI class (i.e., 16) will be available based on the probability of the four variables if the long-term dataset of LSF and clear-sky occurrences is given in the future (Table 4).

Results
The processes for deriving the optimal thresholds and their verification will be discussed in Sections 3.1 and 3.2, based on two-year observations and the RTM simulation. These sections explain the selection of optimal parameters for fog detection. The process of developing PI and its application to case studies to address the PI merits in LSF detection are presented in Section 3.3.

OptimalLSF Thresholds from the 2018 Control Data
To derive optimal LSF thresholds, the average values with standard deviation (±1σ) of satellite observations were investigated for LSF and clear-sky cases at dawn during the control (or training) period of 2018 ( Figure 4). The seven satellite-observed values (BTD 13.5-8.5 , R 0.65 , ∆NDSI, ∆BTD 3.7-11 , NDSI, ∆R 0.65 , and BTD 3.7-11 ) and the values of LST-BT 11 with ground observations of land surface temperature were analyzed to determine the range of each variable in LSF and clear-sky conditions. With the separation of LSF from the clear sky, the LSF threshold can be relatively reliable for LSF detection. Separation without overlap of the weather phenomena was noted for R 0.65 , which was traditionally utilized for fog detection. In addition to the two variables (∆NDSI and ∆BTD 3.7-11 ) derived from the dual-satellites of Himawari-8 and FY-4A, LSF and clear sky were distinct in the BTD 13.5-8.5 available from Himawari-8 only. This result suggests that the threshold values derived from dual satellites (∆NDSI and ∆BTD 3.7-11 ) can be worthy for LSF detection in addition to single satellite observations (BTD 13.5-8.5 and R 0.65 ). However, other thresholds (LST-BT 11 , BTD 3.7-11, ∆R 0.65 , and NDSI) have uncertainties for LSF detection due to mixed information that cannot discriminate between LSF and clear sky.

Results
The processes for deriving the optimal thresholds and their verification will be discussed in Sections 3.1 and 3.2, based on two-year observations and the RTM simulation. These sections explain the selection of optimal parameters for fog detection. The process of developing PI and its application to case studies to address the PI merits in LSF detection are presented in Section 3.3.

OptimalLSF Thresholds from the 2018 Control Data
To derive optimal LSF thresholds, the average values with standard deviation (±1σ) of satellite observations were investigated for LSF and clear-sky cases at dawn during the control (or training) period of 2018 ( Figure 4). The seven satellite-observed values (BTD13.5-8.5, R0.65, ΔNDSI, ΔBTD3.7-11, NDSI, ΔR0.65, and BTD3. [7][8][9][10][11] and the values of LST-BT11 with ground observations of land surface temperature were analyzed to determine the range of each variable in LSF and clear-sky conditions. With the separation of LSF from the clear sky, the LSF threshold can be relatively reliable for LSF detection. Separation without overlap of the weather phenomena was noted for R0.65, which was traditionally utilized for fog detection. In addition to the two variables (ΔNDSI and ΔBTD3.7-11) derived from the dual-satellites of Himawari-8 and FY-4A, LSF and clear sky were distinct in the BTD13.5-8.5 available from Himawari-8 only. This result suggests that the threshold values derived from dual satellites (ΔNDSI and ΔBTD3.7-11) can be worthy for LSF detection in addition to single satellite observations (BTD13.5-8.5 and R0.65). However, other thresholds (LST-BT11, BTD3.7-11, ΔR0.65, and NDSI) have uncertainties for LSF detection due to mixed information that cannot discriminate between LSF and clear sky. Using the data for the control period of 2018, the optimal thresholds were derived through the iteration process. The lower boundary value (i.e., Thlower in Table 6) that distinguished LSF and clear sky was obtained through iteration, while the upper boundary value (i.e., Thupper) was adjusted to the +1 σ range of the average value to avoid middle or high cloud contamination. The results of statistical verification for each threshold variable are shown in the contingency table and the index definition (Table A3), based on the five skill score indices (HSS, CSI, POD, PC, and FAR) shown in Table 3. HSS can be an overall skill indicator of detection accuracy [52], although the difference between POD and FAR (i.e., POD minus FAR) was analyzed in this study as an important indicator of detection accuracy. The difference values were higher in magnitude order for BTD13.5-8.5, ΔNDSI, R0.65, and ΔBTD3.7-11. In these four variables, the skill scores of POD minus FAR Using the data for the control period of 2018, the optimal thresholds were derived through the iteration process. The lower boundary value (i.e., Th lower in Table 6) that distinguished LSF and clear sky was obtained through iteration, while the upper boundary value (i.e., Th upper ) was adjusted to the +1 σ range of the average value to avoid middle or high cloud contamination. The results of statistical verification for each threshold variable are shown in the contingency table and the index definition (Table A3), based on the five skill score indices (HSS, CSI, POD, PC, and FAR) shown in Table 3. HSS can be an overall skill indicator of detection accuracy [52], although the difference between POD and FAR (i.e., POD minus FAR) was analyzed in this study as an important indicator of detection accuracy. The difference values were higher in magnitude order for BTD 13.5-8.5 , ∆NDSI, R 0.65 , and ∆BTD 3.7-11 . In these four variables, the skill scores of POD minus FAR (0.79-0.82) were roughly the same, within a 3% difference. Meanwhile, the POD minus FAR values for the other four variables did not exceed 0.47. Therefore, the threshold values of the top four variables were more useful for LSF detection than were those of the other four thresholds.  (0.79-0.82) were roughly the same, within a 3% difference. Meanwhile, the POD minus FAR values for the other four variables did not exceed 0.47. Therefore, the threshold values of the top four variables were more useful for LSF detection than were those of the other four thresholds.  To theoretically support these observational results, RTM simulation was performed considering various optical characteristics (e.g., altitude and optical thickness of the fog layer, altitude and optical thickness of upper clouds, effective radii of fog/cloud particles) in Table A2, and the results were shown in Figure 6.  Table A2. To theoretically support these observational results, RTM simulation was performed considering various optical characteristics (e.g., altitude and optical thickness of the fog layer, altitude and optical thickness of upper clouds, effective radii of fog/cloud particles) in Table A2, and the results were shown in Figure 6. The figure shows the simulation results in the domains of a) BTD13.5-8.5 vs. LST-BT11, b) ΔBTD3.7-11 vs. LST-BT11, c) ΔNDSI vs. LST-BT11, and d) R0.65 vs. LST-BT11. The average values of the Himawari-8 angles (i.e., SZA = 79.5°, RAA = 74.5°, and VZA = 41.5°) for 176 LSF cases during 2018 were used as RTM input. The satellite-observed mean values are included in the figure for LSF (black asterisk) and clear sky (gray-cross). The clear-sky values from the observations and simulations (pink asterisk) are depicted using black circles to show the approximate agreement between them. Five colors (yellow, blue, red, navy, and green) denote the heights (FH or CH: 1, 2, 4-6, 8-10 km) of the fog and upper cloud layers with respect to effective radii (ER: 2, 4, 8, 16 μm) and water/ice phase. The RTM details are given in Table  A2.
The DSM variables (i.e., ΔNDSI and ΔBTD3.7-11), as well as BTD13.5-8.5 and R0.65 from the single satellite observations (i.e., Himawari-8 only), had a significant tendency to separate LSF and clear sky at dawn (Figure 6a-d) due to the unique satellite orbit in space (i.e., VZA difference between the two satellites). The distributions of observed mean values for the two meteorological phenomena (LSF and clear sky) were consistent within the ordinate range of the simulation. The difference in "LST-BT11" between observations (5.49 K) and simulations (0.82 K) for clear-sky conditions could be explained in terms of the effect of water vapor and optically thin clouds during the warm season [53]. Considering that the LST-BT11 value is 19.7 K for LSF and 5.5 K for clear sky, the average LSF top height is 2-2.5 km, consistent with that of [11].  However, the simulated values of ΔBTD3.7-11 on the abscissa were somewhat out of range with the average of its observations (Figure 6b). According to the simulation, the upper clouds above the fog layer could cause errors in the LSF detection of satellite observations. The upper clouds generally show relatively large values of LST-BT11, compared to the lower clouds. Overall, the simulations were in good agreement with fog detection in DSM. The DSM variables (i.e., ∆NDSI and ∆BTD 3.7-11 ), as well as BTD 13.5-8.5 and R 0.65 from the single satellite observations (i.e., Himawari-8 only), had a significant tendency to separate LSF and clear sky at dawn (Figure 6a-d) due to the unique satellite orbit in space (i.e., VZA difference between the two satellites). The distributions of observed mean values for the two meteorological phenomena (LSF and clear sky) were consistent within the ordinate range of the simulation. The difference in "LST-BT 11 " between observations (5.49 K) and simulations (0.82 K) for clear-sky conditions could be explained in terms of the effect of water vapor and optically thin clouds during the warm season [53]. Considering that the LST-BT 11 value is 19.7 K for LSF and 5.5 K for clear sky, the average LSF top height is 2-2.5 km, consistent with that of [11].
However, the simulated values of ∆BTD 3.7-11 on the abscissa were somewhat out of range with the average of its observations (Figure 6b). According to the simulation, the upper clouds above the fog layer could cause errors in the LSF detection of satellite observations. The upper clouds generally show relatively large values of LST-BT 11 , compared to the lower clouds. Overall, the simulations were in good agreement with fog detection in DSM.

Verification of Satellite-Observed Thresholds during the Experimental Period of 2019
The eight LSF threshold values derived from the control period of 2018 for LSF detection at dawn were applied to the experimental data during 2019. Figure 7 shows statistical verification for the threshold values of eight satellite-observed variables in terms of a) POD minus FAR, b) POD, c) FAR, and d) CSI. The verification was performed based on the total ground-observations of 171 fog and 110 clear-sky occurrences at dawn near Japan. The thresholds of the top four variables (BTD 13.5-8.5 , R 0.65 , ∆NDSI, and ∆BTD 3.7-11 ) in the verification skill-scores were substantially better than those of the other four variables. The "POD minus FAR" scores of the top four variables during the experimental period were lower by~5% than those during the training period (Table 3). In particular, the "∆NDSI" threshold was most sensitive to the interannual dependence of threshold accuracy. However, more long-term data seemed to be required to estimate interannual variability in view of the two years of data in this study. Figure 8 presents the "±1σ" range of LSF and clear-sky cases for each variable in the same way as Figure 4 so that the biennial changes in satellite-observed variables under the weather phenomena were obtained from a comparison of the results shown in the two Figures. Since the threshold values were optimized for the training period data, statistical verification results decreased slightly (Figure 7 and Table 3). However, there were no significant differences in skill scores for the top four variables (BTD 13.5-8.5 , ∆BTD 3.7-11 , ∆NDSI, and R 0.65 ) between the two years. In addition, the interannual change of NDSI among the entire eight variables was noticeable. Although the skill scores of NDSI itself were not good, ∆NDSI showed excellent results due to the DSM advantage. This result was due to not only the large VZA difference between the two satellites but also to the relatively long optical path of FY-4A, which was viewed from outside of the nadir. Therefore, the DSM was likely to cancel the noise error and amplify the signal for LSF detection. This tendency was similar to that of the ∆BTD 3.7-11 case in a previous study [36].

PI development and Case Study
The probability index of LSF was developed to present fog occurrence as a probability concept rather than an alternative one (either no fog or fog), using the top four variables that showed excellent results (POD~0.8; POD minus FAR~0.7) during the experimental period. In Table 4, five LSF classes were derived from 16 probability sets of the variables based on the threshold test for detection of either LSF or clear-sky phenomenon. For instance, Case 1, which passed the test for LSF detection by all four variables, corresponded to Class 1 in terms of the set "P(A ∩ B ∩ C ∩ D)" with an assigned probability of 1.0. For Class 1, in which all four variables passed their own thresholds, the probability of LSF occurrence was classified as "very high". In the same way, three threshold values among the four variables were achieved for those in Class 2. In addition, any combination of the three variables was classified as a "high" possibility of LSF occurrence. In this regard, the number of cases that passed the threshold test for the four satellite-observed variables (i.e., BTD 13.5-8.5 in gray rectangle, ∆NDSI in orange, ∆BTD 3.7-11 in blue, and R 0.65 in green) was shown in the detection of a) LSF and b) clear sky in the so-called Venn diagram shown in Figure 9. In other words, the zones of red (Class 1 for the occurrence possibility of the two weather phenomena), yellow (Class 2), blue (Class 3), and white (Class 4) were defined and were also explained in the Table. Sens. 2021, 13, 1042 14 of 22 ΔNDSI, and R0.65) between the two years. In addition, the interannual change of NDSI among the entire eight variables was noticeable. Although the skill scores of NDSI itself were not good, ΔNDSI showed excellent results due to the DSM advantage. This result was due to not only the large VZA difference between the two satellites but also to the relatively long optical path of FY-4A, which was viewed from outside of the nadir. Therefore, the DSM was likely to cancel the noise error and amplify the signal for LSF detection. This tendency was similar to that of the ΔBTD3.7-11 case in a previous study [36].   ΔNDSI, and R0.65) between the two years. In addition, the interannual change of NDSI among the entire eight variables was noticeable. Although the skill scores of NDSI itself were not good, ΔNDSI showed excellent results due to the DSM advantage. This result was due to not only the large VZA difference between the two satellites but also to the relatively long optical path of FY-4A, which was viewed from outside of the nadir. Therefore, the DSM was likely to cancel the noise error and amplify the signal for LSF detection. This tendency was similar to that of the ΔBTD3.7-11 case in a previous study [36].    (Table 6). This indicated excellent detection of realistic PI values in the spatial distribution, rather than either "fog" or "no fog" from the traditional method. In other words, when LSF detection was performed using a single threshold, the existence of LSF was generally classified with a probability of either 0 or 1. Meanwhile, the PI in this study showed five classified probability values based on the four satellite-observed variables in DSM, for useful and reliable LSF detection.  Table 4.  Table 4.
When the LSF class in Table 4 was determined for a total of 171 LSF cases during the experimental period of 2019, 99 passed all four thresholds (i.e., Class 1). In addition, 43 cases passed three thresholds (i.e., Class 2). These results are described in detail in Table 6. Satellite-derived PI distributions of the five classes in the Table were weighted at 0.25 probability intervals for detection of LSF and clear-sky cases in 2019, and detection was quantitatively presented in five steps rather than two (i.e., yes or no). The probability values were also shown in terms of "assigned probability" in Table 4. If the PI results were evaluated for LSF cases during 2019, the 171 fog occurrences represented a detection score of 140.25. In clear-sky conditions, the score was 93.75 of 110 (Table 6). This indicated excellent detection of realistic PI values in the spatial distribution, rather than either "fog" or "no fog" from the traditional method. In other words, when LSF detection was performed using a single threshold, the existence of LSF was generally classified with a probability of either 0 or 1. Meanwhile, the PI in this study showed five classified probability values based on the four satellite-observed variables in DSM, for useful and reliable LSF detection.
The PI in the case study was applied to two foggy scenes at dawn near Japan (128-145 E, 30-45 N) to compare the results with those of the conventional method (i.e., using a single threshold value of either R 0.65 or BTD 3.7-11 ) in the LSF spatial distribution. Figure 10 shows the fog probability distribution at dawn (06:00 LST) on July 9, 2019, near Japan using the detection methods of a) R 0.65 threshold, b) BTD 3.7-11 threshold, and c) LSF PI in this study. The BT 11 distribution for approximate cloud-top temperature and the SYNOP map for station-observed fog occurrences were shown in Figure 10d,e, respectively. Fog occurrences (pink triangles in Figure 10a-d) were reported (station numbers 20, 26, and 31 in Table 2) at three stations by ASOS (automated surface observing system).
The R 0.65 threshold successfully detected two of the three fog occurrences at a detection rate of 0.66 (Figure 10a); the rate was 0.33 with the BTD 3.7-11 threshold (Figure 10b). The PI produced by combining the four variables using dual-satellite observations showed a probability of 0.66, higher than that in the BTD 3.7-11 case but similar to that in the R 0.65 case. Although PI accuracy was the same as that of the R 0.65 , it provided LSF spatial distribution in terms of probability in five classes. To examine the LSF-related information over the study area, the BT 11 distribution and the weather map chart were analyzed simultaneously (Figure 10d,e). The BT 11 was useful to estimate the approximate cloud top temperature. In addition, the map allowed the identification of the distributions of foggy and humid areas and lower or upper clouds. The BT 11 values were either high (≥~290 K) or low (≤~230 K) in areas where PI indicated no or low probability of fog due to almost clear-sky or optically thin cloud conditions. In Figure 10c, the LSF PI values indicated that the LSF occurrences were "very high" or "high" in areas where fog actually occurred on the weather map (pink triangle). For areas at station numbers 20 and 31 where fog detection failed using R 0.65 or BTD 3.7-11 , PI indicated a relatively low probability of fog but was the only method to generally succeed in the detection (PI = 0.75). distribution in terms of probability in five classes. To examine the LSF-related infor over the study area, the BT11 distribution and the weather map chart were an simultaneously (Figure 10d,e). The BT11 was useful to estimate the approximate clo temperature. In addition, the map allowed the identification of the distributions o and humid areas and lower or upper clouds. The BT11 values were either high (≥ ~ or low (≤ ~230 K) in areas where PI indicated no or low probability of fog due to clear-sky or optically thin cloud conditions. In Figure 10c, the LSF PI values indicat the LSF occurrences were "very high" or "high" in areas where fog actually occur the weather map (pink triangle). For areas at station numbers 20 and 31 whe detection failed using R0.65 or BTD3.7-11, PI indicated a relatively low probability of f was the only method to generally succeed in the detection (PI = 0.75).  Figure 11 is the same as Figure 10 except for 05:30 LST on 4 July 2019, and wi SYNOP map (not available at the time). According to the ASOS report, there wer fog occurrences at the stations (in Table 2; 6, 16, and 27). The time difference betwe Figure 10. Spatial distributions of fog probability at dawn (06:00 LST) on 9 July 2019, near Japan from detection methods of (a) R 0.65 threshold, (b) BTD 3.7-11 threshold, and (c) LSF PI in this study. (d) BT 11 distribution for cloud-top temperatures, and (e) SYNOP map for station-observed fog occurrences. Fog occurrences at the stations (pink triangles in Figure 9a,b) were reported by ASOS. Figure 11 is the same as Figure 10 except for 05:30 LST on 4 July 2019, and without a SYNOP map (not available at the time). According to the ASOS report, there were three fog occurrences at the stations (in Table 2; 6, 16, and 27). The time difference between the fog occurrence on the report and that on the weather map was greater than 30 minutes, so the low cloud and fog areas were analyzed with the distribution of BT 11 without the map. The PI and R 0.65 successfully detected all three fog occurrences, showing a detection probability of one. However, the BTD 3.7-11 threshold detected only one of the three, a probability of 0.33, only successful at station 6. These results show that PI was able to detect fog more accurately than BTD 3.7-11 and was more realistic (or natural) and excellent than R 0.65 in the LSF spatial distribution. In particular, detailed information on the LSF PI in five classes was useful for aviation and navigation. So far, it has been difficult to identify the occurrence of sea fog in the open ocean despite its frequent occurrence due to the limited number of ground observations. The spatial PI distribution of LSF over a large area could be useful to monitor fog phenomena over an extensive region including the ocean in almost real-time.
in five classes was useful for aviation and navigation. So far, it has been difficult to identify the occurrence of sea fog in the open ocean despite its frequent occurrence due to the limited number of ground observations. The spatial PI distribution of LSF over a large area could be useful to monitor fog phenomena over an extensive region including the ocean in almost real-time.

Discussion
The DSM used in this study is applicable to LSF (or fog) detection based on advanced satellite observations, including GEO-KOMPSAT-2A (GK2A; [54]) in the future, for realtime monitoring. Further improvement in nowcasting and forecasting (e.g., advective fog) can be expected. The Advanced Meteorological Imager (AMI) data of GK2A will be useful for the LSF analysis over the Korean Peninsula because of its nadir point. The DSM from FY-4A and Himawari-8 were applied to the LSF detection near Japan at dawn due to the Himawari-8 nadir. Similarly, the method can be utilized for the detection over the inland

Discussion
The DSM used in this study is applicable to LSF (or fog) detection based on advanced satellite observations, including GEO-KOMPSAT-2A (GK2A; [54]) in the future, for realtime monitoring. Further improvement in nowcasting and forecasting (e.g., advective fog) can be expected. The Advanced Meteorological Imager (AMI) data of GK2A will be useful for the LSF analysis over the Korean Peninsula because of its nadir point. The DSM from FY-4A and Himawari-8 were applied to the LSF detection near Japan at dawn due to the Himawari-8 nadir. Similarly, the method can be utilized for the detection over the inland region of China under the FY-4A nadir at dusk, although LSF generally occurs at dusk less than at dawn due to daytime solar heating.
Various applications of the DSM in the meteorological and environmental fields are possible to observe fog targets in a three-dimensional structure from a stereoscopic view, similar to observations with two eyes [37]. Thus, this method can be applied to any worldwide region where simultaneous observations by two geostationary-orbit satellites are available without additional expense. However, the difference in VZA between the two satellites is required to be between 40 • and 50 • to achieve optimal LSF detection. For this reason, the detection accuracy of ∆R 0.65 , which had been excellent in the COMS and FY-2D (VZA difference~46.5 • ) [9], was significantly reduced in this study. It seems that ∆R 0.65 is more sensitive to both VZA and seasonality than to the top-four variables in this study. Here the variables of ∆NDSI and ∆BTD 3.7-11 were newly utilized for LSF detection in terms of PI, despite the decrease in VZA, with the BTD 13.5-8.5 from the spectral difference of Himawari-8 only.
The RTM simulation of this study may have some limitations in the diverse weather conditions of multilayer, broken, or mixed-phase clouds [9]. Also, errors in the LSF detection due to the higher (or upper) clouds above the LSF layer have been discussed in detail, based on the simulation [9,36]. However, the adverse issues may not be so critical for the LSF (under~2.5km height) in relatively stable atmospheric conditions, considering the ratio (13-20%) of satellite-observed multi-layer cloudy pixels to all cloudy pixels [55] and the vertical properties of ice-over-water clouds [56]. In addition, for a comparison of the two GEO data in this study, the coarse-grid averaging of Himawari-8 observations from 2 km-by-2 km to 4 km-by-4 km may miss the detection of optically thin LSF phenomena in a small horizontal scale less than 4 km, particularly in the tests of BTD 13.5-8.5 and R 0.65 .

Conclusions
In this study, the dual-satellite method was applied to advanced satellites (Himawari-8 and FY-4A), and newly developed thresholds with additional channels and LSF spatial distribution in terms of PI were provided to improve its detection in Japan at dawn. In addition, the optical characteristics of LSF were estimated using the LUT from the RTM simulation under various conditions and were consistent with satellite-derived observations. The PI of LSF in the DSM was derived by combining the major threshold components (i.e., BTD 13.5-8.5 , ∆NDSI, ∆BTD 3.7-11 , and R 0.65 ) based on the top-four "POD minus FAR" scores (≥ 79%) among the eight variables.
The case study was carried out to address the PI merits of the LSF spatial distribution. The PI values in this study were similar to R 0.65 simply based on the success rate of fog detection due to the use of R 0.65 as one of the four equally-weighted variables for constructing the PI. Thus, the fog detection rate calculated from PI was expected to be similar to that from R 0.65 . However, the spatial distributions revealed the advantages of PI. Since multi-stage fog detection was possible in PI over a wide area, the LSF distribution in the DSM could be presented in more detail than that of the conventional methods of R 0.65 and/or BTD 3.7-11 . More specifically, R 0.65 only reported two classes of either presence or absence of fog, whereas PI presented fog occurrence probability as five classes. The classes can be increased to a maximum of 16 if a long-term dataset is available in the future, as shown in the statistical combination (Table 4). In this case, different probability weights would be assigned with respect to the variables. Thus, the LSF PI can be used as more practical and useful information for navigation and aviation. In particular, we believe that LSF PI is valuable for open oceanic areas where ground observations are limited in view of both accuracy and spatial distribution.

Acknowledgments:
We thank the employees at the CMA (for FY-4A) and the NMSC/KMA (for Himawari-8) for providing the satellite data. We are grateful to the anonymous reviewers for the constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. List of acronyms used in this study.

Acronyms
Original