Development of Land Surface Albedo Algorithm for the GK-2A/AMI Instrument

The Korea Meteorological Administration successfully launched Korea’s next-generation meteorological satellite, Geo-KOMPSAT-2A (GK-2A), on 5 December 2018. It belongs to the new generation of GEO (Geostationary Elevation Orbit) satellite which offers capabilities to disseminate high spatial- (0.5–2 km) and high temporal-resolution (10 min) observations over a broad area, herein a geographic disk encompassing the Asia–Oceania region. The targeted objective is to enhance our understanding of climate change, owing to a bulk of coherent observations. For such, we developed an algorithm to map the land surface albedo (LSA), which is a major Essential Climate Variable (ECV). The retrieval algorithm devoted to GK-2A/Advanced Meteorological Imager (AMI) data considered Japan’s Himawari-8/Advanced Himawari Imager (AHI) data for prototyping, as this latter owns similar specifications to AMI. Our proposed algorithm is decomposed in three major steps: atmospheric correction, bidirectional reflectance distribution function (BRDF) modeling and angular integration, and narrow-to-broadband conversion. To perform BRDF modeling, the optimization method using normalized reflectance was applied, which improved the quality of BRDF modeling results, particularly when the number of observations was less than 15. A quality assessment was performed to compare our results to those of Moderate Resolution Imaging Spectroradiometer (MODIS) LSA products and ground measurement from Aerosol Robotic Network (AERONET) sites, Australian and New Zealand flux tower network (OzFlux) site and the Korea Flux Network (KoFlux) site from throughout 2017. Our results show dependable spatial and temporal consistency with MODIS broadband LSA data, and rapid changes in LSA due to snowfall and snow melting were well expressed in the temporal profile of our results. Our outcomes also show good agreement with the ground measurements from AERONET, OzFlux and KoFlux ground-based network with root mean square errors (RMSE) of 0.0223 and 0.0306, respectively, which is close to the accuracy of MODIS broadband LSA. Moreover, our results reveal still more reliable LSA products even when clouds are frequently present, such as during the summer monsoon season. It shows that our results are useful for continuous LSA monitoring.


Introduction
LSA (abbreviations not defined in the text and can be found in Table 1) represents the ratio of solar radiance reflected from land surface across the entire spectrum of solar radiation. Therefore, it is found useful to quantify the amount of solar radiance absorbed by the Earth's land surface [1][2][3]. As one of the ECVs, it functions as key forcing data for radiation budget, land surface climate, and biosphere models [4]. Changes in surface albedo due to artificial (e.g., deforestation, urbanization and grazing of rangelands) and natural (e.g., changing seasons, snowfall, wild fire, and desertification) factors are important indicators of land-cover and land-use changes in close relationship with regional climate systems as well as global climate change [5][6][7]. Therefore, accurate assessment of LSA is mandatory as it contributes significantly to characterizing radiation budgets and improving weather and climate model simulations [8]. In practice, some studies have shown that the use of dynamic albedo products in climate and weather forecasting models can help improve predictions of various weather variables, such as surface temperature, latent and sensible heat flux, and net ecosystem exchange [9].
Because the quantification of LSA is closely associated with the appraisal of several climate effect phenomena, LSA has been observed in various ways in recent decades [10][11][12][13][14]. Satellite remote sensing, among others, is a unique and effective tool that provides globally consistent LSA characteristics at high spatial and temporal resolution [15]. Early LSA calculations were as simple as finding the minimum value of TOA reflectance as observed from the satellite [2]. However, advances in numerical computational capabilities and radiation transfer theory have led to the development of the LSA retrieval algorithm using more sophisticated methods, such as the detection of cloud cover in satellite images and a correction of the effects of scattering and absorption by atmospheric components in the path of sunlight [16].
For polar orbit satellite, most LSA retrieval algorithm are divided into two types. One is a traditional algorithm that is composed of the most intuitive methods (atmospheric correction, BRDF modeling, and N2B conversion) for estimating LSA from satellite measurement. It is used in operating system for LSA retrieval from MODIS [17] and PROBA-V [18]. The other is a direct estimation algorithm. This algorithm produces directly LSA from single TOA reflectance without performing atmospheric correction and BRDF modeling to reduce the impact of errors occurring in each process of the traditional algorithm [15]. Nevertheless, atmospheric effect and surface BRDF are considered in the preparation of input data to establish the relationship between TOA reflectance and LSA [19]. This method was applied to the operational LSA retrieval algorithm for VIIRS [20] GLASS [5].
The initial method of calculating LSA from geostationary satellite observations was quite different from that for polar orbit [15]. The sensor onboard first-generation geostationary satellites (e.g., MVIRI) have only one broad visible channel; thus, there are problems to perform the atmospheric correction, especially, when the atmospheric impact is large (e.g., high aerosol loading). To overcome these shortcomings, Pinty et al. [21] designed a new approach for MVIRI LSA. This approach employs RPV model and 6S RTM to solve coupled surface-atmosphere radiation transfer problem and can estimate surface anisotropy and aerosol load simultaneously. This approach is currently utilized to geostationary surface albedo algorithm [22]. However, due to technical limitations on the near-real-time operation, the traditional algorithm as a simplified approach is used to estimate daily basis LSA from SEVIRI measurements [23]. In this method, the additional atmospheric variables were used to perform the atmospheric correction. The TPW and TCO were obtained from ECMWF and TOMS climatology, respectively, and AOD was determined as a function of latitude. In addition, several studies were conducted to improve each process of constructing the traditional algorithm applied to geostationary satellite. Proud et al. [24] present the method to improve the performance of atmospheric correction through modification of SMAC. Pokrovsky and Roujean [1] describes optimal scheme for angular Thanks to the aforementioned efforts, in recent years, many LSA products have been generated using polar orbit and tools that include the MODIS [27], AVHRR [11], MISR [28], VEGETATION [29] as well as geostationary (e.g., MVIRI [21], SEVIRI [23], GOES [26]) satellite-mounted sensors [19].
In particular, geostationary satellites are suitable for LSA estimations because they periodically observe the same wide area at high temporal resolution [30]. Unlike polar orbit satellites, this has the asset of collecting observations on the same land surface multiple times during a day and obtaining a multi-angle sampling dataset every day. As the number of surface observations increases from sunrise to sunset, it is expected to improve and stabilize the performance BRDF modeling, which is mandatory Remote Sens. 2020, 12, 2500 4 of 29 for accurate LSA retrieval [31]. The frequent observations are also possible to improve the completeness of LSA retrievals by facilitating observations in areas with frequent clouds, such as the tropical regions and East Asia during the summer monsoon season [32]. Moreover, the high temporal resolution of GEO satellites shortens the sampling period for BRDF modeling. This allows the satellite-based LSA product to better reflect the surface albedo variability that ensues when hazardous events within a short timeframe due to events such as wildfires, snowfall, and snow melting [20,33]. In light of these advantages, the SCOPE-CM, which is a network of agencies and operators of environmental satellite systems and interfaces, makes an effort to estimate global LSA products as ECVs based on a constellation of GEO satellites through one of the pilot projects [22]. In this respect, GEO satellites are the best tool with which to produce continuous, long-term, and high-quality LSA over a large region for the purpose of studying and monitoring the Earth's surface.
Several agencies worldwide have launched or are planning to launch a new generation of geostationary satellites, including GOES-R in the US, MTG in Europe, Himawari-8 in Japan, and FengYun-4A in China [34][35][36][37]. GK-2A satellite, which is a next-generation geostationary meteorological satellite that replaces South Korea's COMS, was successfully launched on 4 December 2018 [38]. The GK-2A AMI is located at 36,000 km above the equator (central longitude of 128.2E • ) and observes the eastern hemisphere every 10 min to obtain geostationary meteorological data for continuous monitoring of meteorological phenomena in the Asia-Oceania region [39]. Lee et al. [32] proposed an LSA-retrieval algorithm using one broadband visible channel of COMS/MI, but it is limited to the East Asia region. Moreover, COMS/MI provided its observations until March 2020, after which it was terminated.
In this study, we propose a LSA-retrieval algorithm for GK-2A/AMI aimed at achieving sustained and consistent LSA monitoring in the Asia-Oceania region. The spatio-temporal resolution and spectral information of AMIs is superior to that obtained using MI and can improve the quality of LSA products. Because AMI data are available from July 2019 following in-orbit testing, we developed the LSA algorithm using Himawari-8/AHI, which has very similar channel specification and observation cycle to GK-2A/AMI [40], with the benefit to analyze seasonal variation in LSA. To develop the LSA-retrieval algorithm for GK-2A/AMI, we used a traditional algorithm approach suitable for operational purposes. The AOD determined as a function of latitude in the SEVIRI LSA algorithm cannot consider its spatial and temporal variation, thus, our algorithm used real-time AOD data. In addition, optimization process was used to improve the quality of BRDF modeling.
The remainder of this paper is organized as follows. Section 2 describes the data used in the development of the LSA algorithm. The methods used to estimate LSA and for validation are explained in Section 3. In Section 4, we present a quality assessment of our results compared to other satellite-based LSA products and in situ data. Finally, in Section 5, we summarize our findings.

Satellite and Reanalysis Data
Himawari-8, operated by the MSC of the JMA, began operation on 7 July 2015. The AHI on board Himawari-8, as mentioned in Section 1, has a similar channel configuration to that of AMI and the same spatial (0.5-2 km) and temporal (10 min) resolution for the full-disk domain. In addition, the area observed by AHI overlaps most of that observed by the AMI in addition to their similar geometrical solar/surface/satellite relationships [41]. Thus, AHI is well suited to replace the simulated AMI data. Of AHI's 16 channels, we used TOA radiance for five channels in the shortwave region with central wavelengths of 0.47, 0.51, 0.63, 0.86, and 1.61 µm. The 2.3 µm channel data were excluded because the channel was not equipped in the AMI. Because TOA radiance data have a different spatial resolution for each channel, all TOA radiance data were used after resampling at a 2 km resolution. AHI produces full-disk observations every 10 min, representing one of the highest temporal resolutions.
As AHI does not provide TOC reflectance data, we used various auxiliary data to produce TOC reflectance for each shortwave channel. These data are summarized in Table 2. To eliminate cloud-contaminated pixels, we used CLP data from AHI. The CLP product includes information for cloud effective radius, cloud optical thickness, cloud top height, cloud top temperature, cloud type, and QA. Each cloud type pixel in CLP data has 1 of 12 categorical values, including "clear," "unknown," "fill value," and nine cloud types according to the ISCCP classification [42]. TOC reflectance of each channel was calculated only for pixels classified as "clear." The ARP L2 product contains AOD, angstrom exponent, etc. It is produced using the Lambertian assumption-based retrieval algorithm [43]. The AOD over the land was validated by Zhang et al. [44] and shows good agreement with 58 AERONET sites. The CLP and ARP used in this study were obtained through the JAXA Himawari monitor (P-Tree System) website (https://www.eorc.jaxa.jp/ptree/index.html). CLP and APR only cover a limited area with latitude 60 • S-60 • N and longitude 80 • E-160 • W, which is shaded light blue in Figure 1. TOC reflectance was estimated for the pixels corresponding to land. Land pixels were classified using land/water mask data provided by the NMSC. The snow cover product generated by applying the GK-2A SC algorithm from AHI data was used to distinguish snow-covered pixels from snow-free land. The GK-2A SC algorithm is based on the threshold method and dynamic wavelength warping method proposed by Lee et al. [45] and was developed through optimization for AHI [46]. Although TOA radiance, CLP, ARP, and SC data were provided every 10 min, we used hourly data for computation efficiency. Level 2.0 (quality-assured). The AERONET flux data at the surface were validated by Garcia et al. [51]. We calculated LSA using Level 2.0 surface upward and downward flux data, which are inversion products from the Version 3 algorithm, at the 27 anchor stations (Table 3). These data are derived using aerosol properties measurements from AERONET sites and are available at AERONET website (https://aeronet.gsfc.nasa.gov/). The inversion products have uncertainties (random error plus biases) that are associated with uncertainty in input data and inversion process such as AOD measurement and absolute sky radiance calibration; thus, various input and post inversion screening criteria were applied for its quality assurance [52]. Only measurements observed at local noon were used for comparison with satellite-based LSA to minimize uncertainty due to differences in observation time between the two observation methods.

OzFlux
OzFlux is an established research network aimed at providing a consistent view of energy, carbon and water exchange to the Australian, New Zealand, and global ecosystem modeling communities [26]. One site (Boyagin) was used in this study due to data availability. This site is located within an area of wandoo woodland (Table 3) and surrounded by broadacre farming practices. The consistent observations such as fluxes, vegetation dynamics, and hydrology were provided from September 2017. The measurements of downward and upward shortwave radiations by a net radiometer (Model CNR4, Kipp & Zonen) at height of 4m has been subject to quality control and post-processing and were used for validation of our results.

KoFlux
The CRK is part of the KoFlux network. This site is located in the Korean Peninsula (Table 3) and is surrounded by rice paddy. The downward and upward shortwave radiations for the entire sky were observed at a height of 9.6 m with net radiometer (Model CNR4, Kipp & Zonen B.V., Delft, The Netherlands) throughout 2017. More detailed information about the CRK site is provided in Hwang et al. [53]. The shortwave radiation observed on this site is averaged every 30 min. For example, the recoded shortwave radiation at 14:00 is the average of the observations from 13:30 to 14:00. To match the reference time between the CRK site observation data and the satellite data, we used the noon observation data. We also used the clearness index presented by Lee et al. [54] to exclude cloudcontaminated shortwave radiation measurements from further analyses.  The atmospheric condition data (TPW and TCO) were provided by the ECMWF CAMS. The daily data were used because these atmospheric parameters have insufficient sensitivity to TOC reflectance estimation [47].
The MCD43A3 provides black-sky LSAs at local noon and white-sky LSAs for seven spectral band (bands 1-7) and three broad band domains (0.3-0.7 µm for visible, 0.7-5.0 µm for near-infrared, and 0.3-5.0 µm for shortwave) [17,27]. We used the black-and white-sky shortwave LSA with a training of broadband conversion coefficients and a comparison of our results' purpose. The MODIS shortwave LSA was estimated through traditional algorithm [33]. Although uncertainty may occur in the MODIS shortwave LSA due to errors occurring in each processing step, the MODIS shortwave LSA was validated in various studies and results show good agreement with tower measurements for both snow-free and snow-covered conditions [17,48]. The MCD43A2 was also used. These data were delivered in a sinusoidal projection at a 500 m spatial resolution. For the purpose of a straight comparison of MODIS data with our AHI results, these were re-projected to geostationary projection and resampled to 2 km using the MRT. Only pixels with a quality flag of 0 (i.e., good quality) were used for further analyses, with the exception of time-series analyses.

AERONET Sites
AERONET is a ground-based observation network that provides atmospheric aerosol properties from the direct and diffuse solar radiation measured using Cimel sun/sky-radiometers [49]. This network has provided globally distributed, long-term, and continuous data for several decades. The retrieval algorithm has developed gradually from Version 1 to the current Version 3, which performs quality control and cloud screening entirely automatically [50]. The AERONET inversion code provides spectral and broadband fluxes at TOA and surface as well as aerosol parameters for three levels of data quality: Level 1.0 (Unscreened), Level 1.5 (Cloud-screened and quality-controlled), and Level 2.0 (quality-assured). The AERONET flux data at the surface were validated by Garcia et al. [51]. We calculated LSA using Level 2.0 surface upward and downward flux data, which are inversion products from the Version 3 algorithm, at the 27 anchor stations (Table 3). These data are derived using aerosol properties measurements from AERONET sites and are available at AERONET website (https://aeronet.gsfc.nasa.gov/). The inversion products have uncertainties (random error plus biases) that are associated with uncertainty in input data and inversion process such as AOD measurement and absolute sky radiance calibration; thus, various input and post inversion screening criteria were applied for its quality assurance [52]. Only measurements observed at local noon were used for comparison with satellite-based LSA to minimize uncertainty due to differences in observation time between the two observation methods.

OzFlux
OzFlux is an established research network aimed at providing a consistent view of energy, carbon and water exchange to the Australian, New Zealand, and global ecosystem modeling communities [26]. One site (Boyagin) was used in this study due to data availability. This site is located within an area of wandoo woodland (Table 3) and surrounded by broadacre farming practices. The consistent observations such as fluxes, vegetation dynamics, and hydrology were provided from September 2017. The measurements of downward and upward shortwave radiations by a net radiometer (Model CNR4, Kipp & Zonen) at height of 4m has been subject to quality control and post-processing and were used for validation of our results.

KoFlux
The CRK is part of the KoFlux network. This site is located in the Korean Peninsula (Table 3) and is surrounded by rice paddy. The downward and upward shortwave radiations for the entire Remote Sens. 2020, 12, 2500 7 of 29 sky were observed at a height of 9.6 m with net radiometer (Model CNR4, Kipp & Zonen B.V., Delft, The Netherlands) throughout 2017. More detailed information about the CRK site is provided in Hwang et al. [53]. The shortwave radiation observed on this site is averaged every 30 min. For example, the recoded shortwave radiation at 14:00 is the average of the observations from 13:30 to 14:00. To match the reference time between the CRK site observation data and the satellite data, we used the noon observation data. We also used the clearness index presented by Lee et al. [54] to exclude cloud-contaminated shortwave radiation measurements from further analyses.

Algorithm Development
The overall flow of GK-2A algorithm for LSA retrieval is illustrated in Figure 2. The LSA is estimated through three major steps: atmospheric correction, BRDF modeling, and N2B conversion. These processing steps are intuitive methods of solving the problems in LSA estimation using multi-spectral TOA radiance observed by a satellite sensor [19]. This method is also used for LSA retrieval from the observations of operational satellites, such as MODIS [55] and SEVIRI [30]. The proposed LSA-retrieval algorithm can produce TOC reflectance, surface anisotropy (BRDF parameters), and normalized reflectance for each channel as intermediate products in addition to LSA. These products are applicable in many fields, such as cloud detection [56], aerosol effective height retrieval [57], and monitoring vegetation cover [58]. Details for each processing step are described in the following sections.

Step 1: Atmospheric Correction
The radiance measured by satellites encounters perturbation as it passes through the atmosphere. Eliminating these atmospheric effects, such as absorption and scattering by molecules and gasses, is an essential process for estimating LSA. We used the 6S vector code to correct for atmospheric effects in TOA radiance measured by AHI. The 6S code is a radiative transfer model that is accurate for simulating atmospheric effect, compared to other methods, such as RT3, MODTRAN, and SHARM [59]. Further details for 6S are provided by Vermote et al. [60]. The 6S model provides three atmospheric correction coefficients: xa (of inverse of the transmittance), xb (of the scattering term of the atmosphere), and xc (of spherical albedo) for given input (atmospheric, aerosol, spectral, and geometric) conditions [61]. To overcome the disadvantage of 6S, which is the long calculation time for application to satellite data, we previously calculated atmospheric correction coefficients for various input conditions as shown in Table 4 and stored them in LUT form. Atmospheric correction coefficients were calculated with Lambertian surface assumption. This assumption is conflicts with

Step 1: Atmospheric Correction
The radiance measured by satellites encounters perturbation as it passes through the atmosphere. Eliminating these atmospheric effects, such as absorption and scattering by molecules and gasses, is an essential process for estimating LSA. We used the 6S vector code to correct for atmospheric effects in TOA radiance measured by AHI. The 6S code is a radiative transfer model that is accurate for simulating atmospheric effect, compared to other methods, such as RT3, MODTRAN, and SHARM [59]. Further details for 6S are provided by Vermote et al. [60]. The 6S model provides three atmospheric correction coefficients: xa (of inverse of the transmittance), xb (of the scattering term of the atmosphere), and xc (of spherical albedo) for given input (atmospheric, aerosol, spectral, Remote Sens. 2020, 12, 2500 9 of 29 and geometric) conditions [61]. To overcome the disadvantage of 6S, which is the long calculation time for application to satellite data, we previously calculated atmospheric correction coefficients for various input conditions as shown in Table 4 and stored them in LUT form. Atmospheric correction coefficients were calculated with Lambertian surface assumption. This assumption is conflicts with BRDF modeling and can lead the uncertainty in LSA retrieval. However, Lambertian assumption can reduce calculations of 6S and increase the process efficiency of atmospheric correction by reducing the volume of the 6S LUT. The interval of input parameters that make up the 6S LUT for atmospheric correction can lead to errors in the surface reflectance retrieval, which can cause significant positive or negative errors due to the rapid change in atmospheric effects, especially at high SZA above 70 • [40]. Thus, we refine the constructed 6S LUT using the method proposed by Lee et al. [40], which can reduce the error that occurs due to the interval of input parameters, especially high SZA.
We calculated TOC reflectance for each channel using the Equation (1), which is suggested by Vermote et al. [62], for the pixel corresponding to cloud-free, land, and daytime (SZA < 80 • ) where ρ i and L i refer to TOC reflectance (unitless) and TOA radiance (Wm -2 sr -1 µm -1 ) at channel i, respectively; xa i , xb i , and xc i refer to the atmospheric correction coefficient for channel i selected in the 6S LUT for the pixel's observation conditions.

Step 2: BRDF Model Inversion and Angular Integration
The purpose of this step is the estimation of spectral (narrowband) LSA. The surface albedo means the ratio of downward and upward shortwave radiance for all sun and view geometries. The land surface has anisotropic reflection properties, and it depends on illumination direction, viewing direction, canopy structure such as vegetation density and background characteristic, and wavelengths [63]. Thus, in order to estimate narrowband LSA, the process of estimating the anisotropic reflection characteristics of the surface through the BRDF model and integrating it into the range of the sun/view angle for the hemisphere is necessary [19].
The angular integration was performed to consider surface anisotropy for all angular ranges. The spectral (narrowband) LSA for black-sky and white-sky can be expressed as a linear combination of kernel integration and BRDF parameters (Equations (2) and (3)). where In the above equations, α bs,β and α ws,β refer to black-sky and white-sky spectral LSA of channel β. h b,n and h w,n are kernel integration for black-sky and white-sky respectively. θ s , θ v , and ∅ are SZA, VZA, and RAA, respectively. f n and K n,β are the physical kernels and their weights, respectively, which characterize the anisotropic reflection properties of the land surface.
To describe the land surface anisotropy, we used the kernel-driven BRDF approach. In the semi-empirical model, surface bidirectional reflectance was expressed through a linear combination of three (isotropic, geometric, and volumetric) scattering terms [64]. Each term consists of a kernel that expresses the scattering mechanism of the land surface according to the sun-surface-satellite geometry and the parameter that is the weight of the kernel (Equation (6)).
In the above equation, ρ model means modeled reflectance by kernel-driven BRDF model at certain angular condition. K 0 , K 1 , and K 2 refer to BRDF parameters for each scattering term, and f 1 and f 2 are the physically based kernels for geometric and volumetric scattering. Various kernel-driven BRDF models with different kernels have been proposed, but this study used the Roujean model, which was developed by Roujean et al. [65] and is currently used to obtain LSA from SEVIRI observations [30]. This BRDF model consists of a Ross-Thick (volumetric) kernel and a Roujean (geometric) kernel, and each kernel is calculated as follows: The empirical BRDF parameters can be calculated through BRDF inversion. For this, linear square regression was applied and the least three observations with various angular conditions in the clear sky are required. Unfortunately, AHI is not a multi-angle sensor; thus, we used time-series data accumulated within a short timeframe to retrieve BRDF parameters. The scheme for BRDF modeling in this algorithm is illustrated in Figure 3. We set the synthesis period to five days, considering the accuracy, processing efficiency, and ratio of the pixel with filling value due to a lack of observations. This period may be shortened, as we use observations measured every 10 min. BRDF modeling is performed with a daily sliding window, and BRDF parameters are retrieved every day. These results are assigned to the final date of the synthesis period.
BRDF modeling has several advantages in analyzing the land surface properties from satellite data, but it occasionally presents significant errors depending on the degree of variability in the angular distribution (e.g., when directional signatures are insufficient) [66,67]. To reduce the errors in BRDF modeling caused by these factors, Roujean [68] proposed using prior knowledge of BRDF parameters. In practice, the current BRDF modeling method used in the MCD43A3 algorithm (version 6) uses a typical BRDF database that is updated with up-to-the-minute retrieval with full inversion for each pixel as prior information [69].
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 28 The optimization process scheme is detailed in Figure 4. First, was calculated using BRDF parameters derived from BRDF inversion using the linear regression method. Then K0 was modulated with . Lastly, BRDF inversion was performed again for recalculation of K1 and K2. In the second BRDF inversion, K0 was not calculated and was assumed to be a known variable. This process is repeated until the BRDF parameter value stabilizes with little change after that process. The results of the application of the optimization process using BRDF modeling are presented in Section 4.1. The proposed optimization process in this study can performed without any additional information such as prior knowledge.

Step 3: Broadband LSA Calculation
Limited narrowband wavelengths can be observed by multi-spectral sensors, such as AMI and AHI [72]. This spectral information can be combined via N2B conversion to provide broadband information [73], and this method is applied to several satellite sensors to provide broadband LSA [74]. In this study, broadband LSA for the shortwave domain was calculated using linear combinations of the spectral albedos of AHI and their weights (N2B coefficients), as shown in equation (11): where Λ is black-sky (white-sky) broadband LSA, and is black-sky (white-sky) narrowband LSA of channel .
is the N2B coefficient of the corresponding spectral LSA. N2B coefficients for snowfree and snow-covered are determined through linear regression analyses using MCD43A3 as reference data. Each dataset contains 18,304 and 803 data entries for snow-free and snow-covered conditions, respectively, with various land types (forest, grasslands, barren, urban, and shrublands). The N2B coefficients for the snow-free and snow-covered conditions of each AHI spectral LSA and However, this method can constrain the BRDF inversion [70] and is not suitable when the surface condition changes rapidly or the prior knowledge is old. Therefore, in this study, we used normalized reflectance to minimize the errors in BRDF modeling. The normalized reflectance was estimated through the normalization method proposed by Yeom and Kim [71] (Equation 10), which considers the observation characteristics of the geostationary orbit (fixed viewing angles): where ρ norm is the normalized reflectance. ρ model and ρ measured denote reproduced TOC reflectance using BRDF parameters and measured TOC reflectance, respectively. θ s, mean and ∅ mean are the mean SZA and RAA during the BRDF synthesis period. θ s, i and ∅ i are the SZA and RAA of the observation number i. The n means total number of valid observations for synthesis period. The optimization process scheme is detailed in Figure 4. First, ρ norm was calculated using BRDF parameters derived from BRDF inversion using the linear regression method. Then K 0 was modulated with ρ norm . Lastly, BRDF inversion was performed again for recalculation of K 1 and K 2 . In the second BRDF inversion, K 0 was not calculated and was assumed to be a known variable. This process is repeated until the BRDF parameter value stabilizes with little change after that process. The results of the application of the optimization process using BRDF modeling are presented in Section 4.1. The proposed optimization process in this study can performed without any additional information such as prior knowledge.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 28 The optimization process scheme is detailed in Figure 4. First, was calculated using BRDF parameters derived from BRDF inversion using the linear regression method. Then K0 was modulated with . Lastly, BRDF inversion was performed again for recalculation of K1 and K2. In the second BRDF inversion, K0 was not calculated and was assumed to be a known variable. This process is repeated until the BRDF parameter value stabilizes with little change after that process. The results of the application of the optimization process using BRDF modeling are presented in Section 4.1. The proposed optimization process in this study can performed without any additional information such as prior knowledge.

Step 3: Broadband LSA Calculation
Limited narrowband wavelengths can be observed by multi-spectral sensors, such as AMI and AHI [72]. This spectral information can be combined via N2B conversion to provide broadband information [73], and this method is applied to several satellite sensors to provide broadband LSA

Step 3: Broadband LSA Calculation
Limited narrowband wavelengths can be observed by multi-spectral sensors, such as AMI and AHI [72]. This spectral information can be combined via N2B conversion to provide broadband information [73], and this method is applied to several satellite sensors to provide broadband LSA [74]. In this study, broadband LSA for the shortwave domain was calculated using linear combinations of the spectral albedos of AHI and their weights (N2B coefficients), as shown in Equation (11): where Λ is black-sky (white-sky) broadband LSA, and α i is black-sky (white-sky) narrowband LSA of channel i. ω i is the N2B coefficient of the corresponding spectral LSA. N2B coefficients for snow-free and snow-covered are determined through linear regression analyses using MCD43A3 as reference data. Each dataset contains 18,304 and 803 data entries for snow-free and snow-covered conditions, respectively, with various land types (forest, grasslands, barren, urban, and shrublands). The N2B coefficients for the snow-free and snow-covered conditions of each AHI spectral LSA and the accuracy of linear regression are shown in Table 5. The results show that all regression results were reasonably accurate with a high R over 0.9 and low standard error under 0.05. In our algorithm, when the TOC reflectance with snow flag was more than half of the total observations during the BRDF modeling synthesis period, snow broadband LSA was calculated using N2B coefficients for the snow-covered conditions. Table 5. Results of multiple linear regression. ω 0 and ω i are constant and coefficients for each AHI channel i for converting narrowband albedo to broadband albedo. The quality of broadband LSA data was directly affected by the accuracy of the BRDF modeling [75]. Thus, the quality flag of broadband LSA was designated "good" when the number of observations was higher than seven and the RMSE of BRDF modeling was 0.07 or less; otherwise it was designated "bad". This classification method is similar to the quality flag classification method used in the MODIS LSA algorithm.

Quality Assessment
The quality assessment in this study focused on black-sky and white-sky broadband LSA only, although the LSA algorithm we proposed generates several products, including TOC reflectance, BRDF parameters, normalized reflectance, narrowband LSA, and broadband LSA. First, the inter-comparison of AHI and MODIS broadband LSA was carried out to evaluate the overall statistical discrepancy in relation to other satellite-based LSA, spatial consistency, and temporal consistency. Direct comparison with ground measurements from 27 AERONET stations was also performed to evaluate absolute accuracy, and then the AHI results were compared to those of MODIS. The three validation metrics in the form of scalars (RMSE, bias, and R) were used as quantitative indicators of quality. These validation metrics are the most common metrics used to measure accuracy for continuous variables and are expressed as follows: In the above equations, x and y denote reference and estimation, respectively; x and y denote mean of x and y, respectively; and n is the number of collocated data entries used for quality assessment.
The validation results of satellite-based LSA compared to ground measurements using the direct comparison method are not always valid due to differences in spatial representation between measurement systems [76]. TC analyses has been suggested as a means of reducing this effect on validation results. This was introduced by Stoffelen et al. [77] and is a method used to estimate the random error variances from three given collocated datasets that have the same geophysical variables [78]. This method makes it possible to estimate the errors of satellite-based products with wide spatial resolution by correcting the ground measurements' spatial representation errors. For this reason, the TC method has been used to compare various satellite-based geophysical parameters, including LSA, to ground observations [79]. In this analysis, an independent three-measurement system was required, and the RMSE and R of each triplet's measurement system were calculated [80] as follows: where Q ij is the covariance between the LSA from the measurement system i and measurement system j. Each element of the above matrix corresponds to each measurement system of the triplet. In this study, validation metrics calculated using the TC method and validation metrics calculated using direct comparison were compared to determine the effectiveness of the TC method.
In the further analyses, with the exception of time-series analyses, only MODIS and AHI broadband LSA with "good" quality flags were used. The different LSA estimation strategies for snow-covered pixels between AHI and MODIS can also present a difference in broadband LSA [18]. To minimize the error caused by this, when the pixel was snow-covered throughout the entire BRDF modeling synthesis period, it was used to validate snow-covered broadband LSA. However, if the pixel showed no snow cover throughout the entire BRDF modeling synthesis period, it was used for snow-free broadband LSA validation.

Optimization of BRDF Modeling
To determine the iteration number for optimization, we analyzed changes in BRDF parameters, ρ norm , and the reproduced TOC reflectance using BRDF parameters according to the number of optimization processes. These results for three pixels with different cases for channel 3 are shown in Figure 5. The same results were obtained for the other channels and have been omitted for brevity. located at high latitudes in the northern hemisphere. As shown in Figure 5d, all observations of the pixel were made at a similar local time. Therefore, the angle information sampled from this pixel was limited to a fixed view and a narrow solar geometry range (77.93-79.84°). In such cases, an error can occur when BRDF inversion is due to insufficient angular information [54]. In practice, as shown in Figure 5c,d, all three BRDF parameters calculated prior to optimization showed abnormal values, and the reproduced TOC reflectance was less than 0, indicating a value outside the valid range. The BRDF inversion result improved after the optimization process and showed similar improvements to the previous case. The last row (Figure 5e,f) shows the case when BRDF inversion RMSE was slightly increased after optimization from 0.0285 to 0.0369. Even before the optimization was performed, the TOC reflectance reproduced from the BRDF parameters appeared to simulate the surface anisotropy reflectance of the observed TOC reflectance. Even before the optimization was performed, the reproduced TOC reflectance closely simulates the surface anisotropy characteristics of the observed TOC reflectance (Figure 5f). Nevertheless, K0 has a negative value outside the valid range, as shown in Figure 5e. After optimization, K0 showed a positive value, and the reproduced TOC reflectance closely simulated the observation TOC reflectance prior to optimization. These results show that the optimization process using improves errors that occur for various reasons, such as cloud contamination and limited angular range in BRDF inversion. This can be useful when no a prior information is available. The changes in BRDF parameters and normalized reflectivity according to the number of repetitions of optimization change up to three times, but little change occurs subsequently. Moreover, the reproduced TOC reflectance does not show a significant difference depending on the number of optimization iterations. Therefore, considering the processing efficiency and accuracy, it was deemed appropriate to perform the optimization process three times. The first row (Figure 5a,b) shows when the RMSE of BRDF modeling was reduced from 0.1208 to 0.0704 after optimization. As Figure 5a shows, before the optimization process, the value of K 0 was abnormally high (almost 0.4), although this pixel's land type is an evergreen broadleaf forest located in Northern Australia (tropical north). For this reason, the value of K 2 was also abnormally high at over 3. Therefore, the TOC reflectance reproduced using the BRDF parameter without an optimization process shows a higher value than the typical reflectance for that land type, and the variation is also considerable. This error was caused by a cloud classification error that occurred on December 30 at 01:00 UTC (Figure 5b). After the optimization process, K 0 and K 2 decreased. The reproduced TOC reflectance showed better agreement with the TOC reflectance observed in cloudless conditions, and the variability was also reduced. This result showed that optimization processing can help to improve the quality of BRDF inversion when some observations are contaminated by clouds. The second row (Figure 5c,d) shows the results for snow-covered pixels located at high latitudes in the northern hemisphere. As shown in Figure 5d, all observations of the pixel were made at a similar local time. Therefore, the angle information sampled from this pixel was limited to a fixed view and a narrow solar geometry range (77.93-79.84 • ). In such cases, an error can occur when BRDF inversion is due to insufficient angular information [54]. In practice, as shown in Figure 5c,d, all three BRDF parameters calculated prior to optimization showed abnormal values, and the reproduced TOC reflectance was less than 0, indicating a value outside the valid range. The BRDF inversion result improved after the optimization process and showed similar improvements to the previous case. The last row (Figure 5e,f) shows the case when BRDF inversion RMSE was slightly increased after optimization from 0.0285 to 0.0369. Even before the optimization was performed, the TOC reflectance reproduced from the BRDF parameters appeared to simulate the surface anisotropy reflectance of the observed TOC reflectance. Even before the optimization was performed, the reproduced TOC reflectance closely simulates the surface anisotropy characteristics of the observed TOC reflectance (Figure 5f). Nevertheless, K 0 has a negative value outside the valid range, as shown in Figure 5e. After optimization, K 0 showed a positive value, and the reproduced TOC reflectance closely simulated the observation TOC reflectance prior to optimization.
These results show that the optimization process using ρ norm improves errors that occur for various reasons, such as cloud contamination and limited angular range in BRDF inversion. This can be useful when no a prior information is available. The changes in BRDF parameters and normalized reflectivity according to the number of repetitions of optimization change up to three times, but little change occurs subsequently. Moreover, the reproduced TOC reflectance does not show a significant difference depending on the number of optimization iterations. Therefore, considering the processing efficiency and accuracy, it was deemed appropriate to perform the optimization process three times. Figure 6a shows the mean and standard deviations of the RMSE of the reproduced TOC reflectance by BRDF modeling based on the number of observations for each channel before and after the optimization process over one year. Mean and standard deviation of RMSE tend to decrease as the number of observations increases in all channels except channel 5. The mean RMSE of BRDF modeling was lowered in all channels after the optimization process when the number of observations was lower than nine. However, when the number of observations was 15 or more, RMSE increased after the optimization process. However, the increase in RMSE is small, up to 0.0053 in channel 4. This is also a reasonable result when RMSE increases despite the improvement of BRDF inversion, as shown in Figure 5e,f. The standard deviation of RMSE was improved in all cases after the optimization process, however. In particular, the fewer observations, the more remarkable the improvement, and when fewer than six observations were made, the standard deviation decreased by about half. In Figure 6f, pixels with 24 observations or fewer account for more than 50% of the total valid pixels every day. Therefore, the optimization process using ρ norm improved the accuracy of BRDF inversion without a prior knowledge accumulation for most pixels.

Inter-Comparison with MODIS LSA Product
The broadband LSA derived from AHI was compared to that from MODIS. Figure 7a,c show the difference in broadband LSA estimated from AHI and MODIS for black-sky and white-sky under snow-free conditions. Most pixels were located close to the 1:1 line. Although the AHI broadband LSA values are generally higher than those of MODIS, with the positive bias of 0.0017 for black-sky and 0.0034 for white-sky, they are negligible and generally show good agreement overall. The validation metrics of black-sky (RMSE of 0.0176, bias of 0.0017, and R of 0.90) showed slightly improved results compared to those of white-sky (RMSE of 0.0236, bias of 0.0034, and R of 0.84). Moreover, each point of white-sky is more widely distributed from the diagonal than those of black-sky. Figure 7b,d show comparison results of broadband LSA in snow-covered pixels. Most of pixels are located close to the 1:1 line and show good agreement overall. The higher RMSE that is 0.0523 for black-sky and 0.0529 for white-sky in snow-free conditions is reasonable when the high broadband LSA of the snow-covered pixels is considered [3,81]. Unlike the comparison results for snow-free conditions, AHI broadband LSA was lower than that of MODIS with a negative bias of -0.0052 for black-sky and −0.0051 for white-sky. The accuracies of snow-covered pixels for black-sky and white-sky were similar with close validation metrics (RMSE, bias, and R). Consequently, these results indicate that black-sky and white-sky AHI broadband LSA agree well overall with MODIS broadband LSA, which was well validated in both snow-free and snow-covered conditions, although some differences are evident between these sensors, such as cloud detection, snow cover detection, the synthesis period of BRDF modeling, spectral information, and angular sampling.
the optimization process over one year. Mean and standard deviation of RMSE tend to decrease as the number of observations increases in all channels except channel 5. The mean RMSE of BRDF modeling was lowered in all channels after the optimization process when the number of observations was lower than nine. However, when the number of observations was 15 or more, RMSE increased after the optimization process. However, the increase in RMSE is small, up to 0.0053 in channel 4. This is also a reasonable result when RMSE increases despite the improvement of BRDF inversion, as shown in Figure 5e, f. The standard deviation of RMSE was improved in all cases after the optimization process, however. In particular, the fewer observations, the more remarkable the improvement, and when fewer than six observations were made, the standard deviation decreased by about half. In Figure 6f, pixels with 24 observations or fewer account for more than 50% of the total valid pixels every day. Therefore, the optimization process using improved the accuracy of BRDF inversion without a prior knowledge accumulation for most pixels.   Figure 8 shows the statistical difference for one year between AHI and MODIS broadband LSA in geographic distribution. In the vicinity of the equator and in southern China, there were no high-quality MODIS pixels for a year, so it was marked in white to indicate the fill value. As shown in Figure 8a,b, the salient characteristic is that most pixels have a low bias within ±0.025 for both black-sky and white-sky. Relatively large positive or negative bias was shown in the regions of western China. The region is characterized by frequent cloud cover [82]. Thus, these features in bias may be related to the difference between satellite-based cloud retrievals from the two sensors. The regions with high VZAs, such as eastern India, also show a relatively high positive bias. This may also be related to cloud masking of AHI because geostationary satellites perform relatively weakly in cloud masking for regions with large VZA [83]. In the higher latitudes of the northern hemisphere, such as southern Russia and the Kamchatka Peninsula, a slight negative bias is present due to snow cover in the winter season. Generally, the white-sky bias was larger than that of black-sky. In practice, white-sky bias showed a slightly positive bias, unlike black-sky in Australia. This result is in accordance with the results presented in Figure 7. The geographic distribution of RMSE for black-sky and white-sky (Figure 8c,d) shows similar results to those related to bias. Of all valid black-sky and white-sky pixels, 96.41% and 95.73%, respectively, show accuracy within the 0.02-0.05 range that is required as absolute retrieval accuracy by climate models across a range of spatial and temporal scales [84] compared to MODIS. Therefore, AHI broadband LSA has good spatial consistency with MODIS broadband LSA. and 0.0034 for white-sky, they are negligible and generally show good agreement overall. The validation metrics of black-sky (RMSE of 0.0176, bias of 0.0017, and R of 0.90) showed slightly improved results compared to those of white-sky (RMSE of 0.0236, bias of 0.0034, and R of 0.84). Moreover, each point of white-sky is more widely distributed from the diagonal than those of blacksky. Figure 7b,d show comparison results of broadband LSA in snow-covered pixels. Most of pixels are located close to the 1:1 line and show good agreement overall. The higher RMSE that is 0.0523 for black-sky and 0.0529 for white-sky in snow-free conditions is reasonable when the high broadband LSA of the snow-covered pixels is considered [3,81]. Unlike the comparison results for snow-free conditions, AHI broadband LSA was lower than that of MODIS with a negative bias of -0.0052 for black-sky and −0.0051 for white-sky. The accuracies of snow-covered pixels for black-sky and whitesky were similar with close validation metrics (RMSE, bias, and R). Consequently, these results indicate that black-sky and white-sky AHI broadband LSA agree well overall with MODIS broadband LSA, which was well validated in both snow-free and snow-covered conditions, although some differences are evident between these sensors, such as cloud detection, snow cover detection, the synthesis period of BRDF modeling, spectral information, and angular sampling. Figure 7. Density scatter plots for AHI LSA versus MODIS LSA: (a) case with black-sky and snowfree (b) case with black-sky and snow-covered (c) white-sky and snow-free (d) case with white-sky and snow-covered. The color bar denotes the ratio of frequency in 0.005 X 0.005 bins. The terms "slope" and "offset" represent slope and interception of regression line, and "N_data" means the number of data entries used in this analysis. The regression line and 1:1 line are shown as dash and solid gray lines, respectively. Figure 8 shows the statistical difference for one year between AHI and MODIS broadband LSA in geographic distribution. In the vicinity of the equator and in southern China, there were no high- Figure 7. Density scatter plots for AHI LSA versus MODIS LSA: (a) case with black-sky and snow-free (b) case with black-sky and snow-covered (c) white-sky and snow-free (d) case with white-sky and snow-covered. The color bar denotes the ratio of frequency in 0.005 X 0.005 bins. The terms "slope" and "offset" represent slope and interception of regression line, and "N_data" means the number of data entries used in this analysis. The regression line and 1:1 line are shown as dash and solid gray lines, respectively.

Inter-Comparison with MODIS LSA Product
The temporal comparisons of the black-sky broadband LSA with the MODIS data for five land types (barren [41. Figure 9. The results for white-sky were similar to those for black-sky; thus, they have been omitted for brevity. For the barren land (Figure 9a), the LSA shows a stable value with around 0.22 throughout the year. However, the LSA of AHI shows little temporal variability compared to that of MODIS. The misclassification of clouds may occur in the barren region due to its spectral signature, which is similar to cloud cover [85]. This indicates that the error in TOC reflectance due to misclassification of clouds affects the accuracy of the LSA. The broadband LSA of AHI and MODIS showed a difference within ±0.05, indicating a good overall agreement. In the case of pixels corresponding to croplands (Figure 9b), seasonal variation in the LSA was well expressed in both AHI and MODIS. The broadband LSA tends to decrease over time until summer, which corresponds to the vegetation growing season, and in the later period, the broadband LSA increases in association with the aging of vegetation. During some summer periods (DOY 170-199 and 215-231) when the land surface was frequently covered by clouds, broadband LSA was not estimated from MODIS whereas AHI broadband LSA was retrieved. This demonstrates the advantages of AHI's high temporal resolution. This advantage can be maximized by taking observations every 10 min. As Figure 9c illustrates, broadband LSA has a large range from about 0.2 during the snow-free period to around 0.8 during the snow-covered period. The broadband LSA from both AHI and MODIS have slightly different values due to the different broadband LSA calculation strategies for snow-covered pixels but show a similar change during the snow melting period in the spring and snowfall periods in winter. However, at the end of the snow melting season (DOY 89-97), MODIS broadband LSA was high when covered with snow, although it was not classified as a "snow" flag in the MCD43A2. This is thought to be due to a cloud and snow cover detection error in the snow melting season and MODIS' relatively long BRDF synthesis period (16 days). AHI broadband LSA closely followed the changes in the ratio of snow covered observation to valid observation during the BRDF synthesis period. Even in the pixel corresponding to the mixed forest (Figure 9d), the change in AHI albedo in winter (DOY 1-85 and 315-365) corresponds well with the ratio of snow-covered observations. The MODIS broadband LSA was not retrieved during the same period. In the other periods, the AHI broadband LSA shows a slightly larger variation than MODIS in summer. The majority of retrieved LSAs show a difference of 0.05 or less. Figure 9e shows the time-series of LSA from AHI and MODIS for open shrublands. This case is hardly affected by clouds and snow. Thus, AHI broadband LSA closely follows that of MODIS all year round with a difference within ±0.03.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 28 quality MODIS pixels for a year, so it was marked in white to indicate the fill value. As shown in Figure 8a ,b, the salient characteristic is that most pixels have a low bias within ±0.025 for both blacksky and white-sky. Relatively large positive or negative bias was shown in the regions of western China. The region is characterized by frequent cloud cover [82]. Thus, these features in bias may be related to the difference between satellite-based cloud retrievals from the two sensors. The regions with high VZAs, such as eastern India, also show a relatively high positive bias. This may also be related to cloud masking of AHI because geostationary satellites perform relatively weakly in cloud masking for regions with large VZA [83]. In the higher latitudes of the northern hemisphere, such as southern Russia and the Kamchatka Peninsula, a slight negative bias is present due to snow cover in the winter season. Generally, the white-sky bias was larger than that of black-sky. In practice, whitesky bias showed a slightly positive bias, unlike black-sky in Australia. This result is in accordance with the results presented in Figure 7. The geographic distribution of RMSE for black-sky and whitesky (Figure 8c,d) shows similar results to those related to bias. Of all valid black-sky and white-sky pixels, 96.41% and 95.73%, respectively, show accuracy within the 0.02-0.05 range that is required as absolute retrieval accuracy by climate models across a range of spatial and temporal scales [84] compared to MODIS. Therefore, AHI broadband LSA has good spatial consistency with MODIS broadband LSA.   1-85 and 315-365) corresponds well with the ratio of snow-covered observations. The MODIS broadband LSA was not retrieved during the same period. In the other periods, the AHI broadband LSA shows a slightly larger variation than MODIS in summer. The majority of retrieved LSAs show a difference of 0.05 or less. Figure 9e shows the time-series of LSA from AHI and MODIS for open shrublands. This case is hardly affected by clouds and snow. Thus, AHI broadband LSA closely follows that of MODIS all year round with a difference within ±0.03. The above result shows that AHI broadband LSA was temporarily consistent with that of MODIS. Furthermore, AHI broadband LSA expresses seasonal variation well in addition to rapid changes due to snowfall and snowmelt, while the temporal variation in AHI was slightly higher than that of MODIS. The AHI can provide broadband LSA with improved temporal completeness thanks to its high temporal resolution in areas with frequent cloud or where snow and cloud are difficult to distinguish. When all observations taken every 10 min are used to estimate the LSA, it is expected that the temporal completeness will improve further and that the quality of the data produced will also improve.
The broadband LSA varies somewhat according to land type and time. The AHI broadband LSA was compared to that of MODIS for some land types and months ( Figure 10). The forest type indicates that all five land types (evergreen broadleaf, evergreen needle-leaf, deciduous broadleaf, deciduous needle-leaf, and mixed forest) are combined into one. Both black-sky and white-sky exhibit low RMSE (from 0.0157 to 0.0270 in black-sky and from 0.0216 to 0.0302 in white-sky) and bias (from −0.0042 to 0.0066 in black-sky and from −0.0117 to 0.0014 in white-sky) for all land types (Figure 10a,b). The variation in bias in different land types was slightly greater in white-sky than black-sky, but it was nonetheless negligible. In Figure 10c,d, the RMSE between AHI and MODIS broadband LSA was almost stable in black-sky (from 0.0160 to 0.0203), while the monthly RMSE in white-sky showed a slightly larger variation according to time (from 0.0180 to 0.0287). The reason for this error tendency may be that the error generated in the atmospheric correction and BRDF modeling step was amplified by the integration of the BRDF kernel for the SZA to calculate the white-sky spectral LSA [32]. Nevertheless, both black-sky and white-sky show reasonable differences with an RMSE smaller than 0.05. nonetheless negligible. In Figure 10c,d, the RMSE between AHI and MODIS broadband LSA was almost stable in black-sky (from 0.0160 to 0.0203), while the monthly RMSE in white-sky showed a slightly larger variation according to time (from 0.0180 to 0.0287). The reason for this error tendency may be that the error generated in the atmospheric correction and BRDF modeling step was amplified by the integration of the BRDF kernel for the SZA to calculate the white-sky spectral LSA [32]. Nevertheless, both black-sky and white-sky show reasonable differences with an RMSE smaller than 0.05. Figure 10. The RMSE (blue) and mean bias (red) of black-sky and white-sky AHI broadband LSA compared to MODIS broadband LSA according to land types and month: (a) RMSE according to land types for black-sky broadband LSA (b) RMSE according to land types for white-sky broadband LSA (c) RMSE according to month for black-sky broadband LSA (d) RMSE according to land types for white-sky broadband LSA.

Accuracy Evaluation Using Ground Measurements
The results of validation using all available spatial-temporal matchup data between satellite-based LSA and ground measurements from AERONET sites and Boyagin site are shown in Table 6. In general, AHI shows the tendency to overestimate with a bias of 0.0125. This leads to an RMSE that is somewhat slightly worse in AHI (0.0223) than that in MODIS (0.0195). At some sites, such as "Beijing-CAMS", "Beijing", "Dongsha_Island", "Lumbini", and "Yonsei_University", bias over 0.025 is indicated and shows a relatively higher magnitude of overestimation than other sites. These sites have over three times more retrievals in the case of AHI than in the case of MODIS in common. Therefore, this trend of overestimation may be caused by differences in the cloud mask between the two sensors. Nonetheless, AHI broadband LSA outlines good overall consistency with ground measurements, which is comparable to the validation results of [26]. Moreover, in the case of the Boyagin site, MODIS shows a greater positive bias than AHI. In most sites (except three: "Irkutsk," "Omkoi," and "Pusan_NU"), the number of AHI samples was higher than that of MODIS, and the total number of AHI samples was 541, far larger than 317. This indicates that AHI has better completeness and can provide a larger number of valid pixels in the Asia-Oceania region, thanks to its high temporal resolution.
The histogram of the difference between the ground measurements and the satellite-based LSA from AHI and MODIS ( Figure 11) has almost normal distribution for AHI and right-skewed distribution for MODIS. Although the AHI tends to overestimate compared to the ground measurement with a positive mean (0.0125) and median (0.0125) of difference, 97.4% of collocation data show an absolute difference within 0.05. The largest difference was 0.079 at the Beijing site, which may be attributable to an error in TOC reflectance by large aerosol loading [26]. In addition, AHI has a slightly better mode (i.e., the value where the histogram reaches its peak) than MODIS. As mentioned, when comparing a satellite-based LSA with a kilometer-scale spatial resolution to ground observations, some errors may occur due to differences in spatial representation between the data. We conducted an evaluation using the TC method to address the issue related to spatial representation, and then compared the results with the results achieved using a direct comparison method ( Table 7). The results of the direct comparison presented in Table 7 differ from the results shown in Table 6 and Figure 11 because these analyses were used only for comparison when all three data sets (ground measurement, AHI, and MODIS) were valid. For both AHI and MODIS, the TC method yielded slightly better RMSE and R than direct comparison. This indicates that the TC method can help to compare between two data sets with different spatial representations. As mentioned by Wu et al. [80], because the TC technique does not completely correct for the difference in spatial representation between the ground measurement and the satellite data, no significant improvement in RMSE or R 2 was observed for either AHI or MODIS. Another reason for this result is that the TC technique involves several assumptions, which are not always met [79]. As mentioned, when comparing a satellite-based LSA with a kilometer-scale spatial resolution to ground observations, some errors may occur due to differences in spatial representation between the data. We conducted an evaluation using the TC method to address the issue related to spatial representation, and then compared the results with the results achieved using a direct comparison method ( Table 7). The results of the direct comparison presented in Table 7 differ from the results shown in Table 6 and Figure 11 because these analyses were used only for comparison when all three data sets (ground measurement, AHI, and MODIS) were valid. For both AHI and MODIS, the TC method yielded slightly better RMSE and R than direct comparison. This indicates that the TC method can help to compare between two data sets with different spatial representations. As mentioned by Wu et al. [80], because the TC technique does not completely correct for the difference in spatial representation between the ground measurement and the satellite data, no significant improvement in RMSE or R 2 was observed for either AHI or MODIS. Another reason for this result is that the TC technique involves several assumptions, which are not always met [79]. Figure 12 shows the time-series of satellite-derived broadband LSA and ground measurements for the CRK site. The AHI and MODIS broadband LSA reasonably captured the seasonal variation in   Figure 12 shows the time-series of satellite-derived broadband LSA and ground measurements for the CRK site. The AHI and MODIS broadband LSA reasonably captured the seasonal variation in ground measurements. In addition, the RMSE of AHI and MODIS broadband LSA was 0.0309 and 0.0316, respectively, which is less than the required accuracy for climate models. Both AHI and MODIS broadband LSA showed negative bias of −0.0193 and −0.0229, respectively. This error may be due to following two reason: (1) the difference in the reference time between the ground measurement from the CRK site (average from local noon to the next 30 min) and the satellite-based broadband LSA (instantaneous at local noon), (2) difference in spatial representation between two measurement systems. The AHI broadband LSA value was greater than that of MODIS throughout 2017 due to its high temporal resolution. This advantage is particularly pronounced during the summer months when it is hot and humid and the rain can persist for long periods. During the DOY 170-236, which is usually a long rainy season, AHI can provide broadband LSA while MODIS has no LSA retrieval Remote Sens. 2020, 12, 2500 23 of 29 capabilities. This indicates that AHI is more useful for the retrieval of broadband LSA when clouds are frequently present.
capabilities. This indicates that AHI is more useful for the retrieval of broadband LSA when clouds are frequently present.

Conclusions
In this study, we developed an LSA-retrieval algorithm with the operational purpose of retrieving LSA in clear-sky conditions from GK-2A/AMI, South Korea's next-generation geostationary meteorological satellite. Due to GK-2A/AMI's short operation period, the proposed algorithm was developed using Himawari-8/AHI, which has similar channel specifications to GK-2A/AMI in the shortwave domain. The LSA was retrieved through three major steps: atmospheric correction, BRDF modeling, and N2B conversion. To perform each processing step, we adopted the 6S RTM, the Roujean BRDF model, and the multiple linear regression method and optimized the geostationary satellite observation system. In particular, the iterative optimization process used in BRDF modeling would work well to improve the accuracy of BRDF inversions, and the improvement was evident when the number of observations taken during the BRDF composite period was below nine. The multiple linear regression results for deriving N2B conversion coefficients were also reliable for both snow-free and snow-covered conditions with a low standard error of less than 0.05 and a high R of over 0.9.
The black-sky and white-sky broadband LSA retrieved from AHI was compared to MODIS LSA product to evaluate its overall quality over the period of a year. In general, black-sky (white-sky) broadband LSA from AHI shows good overall consistency with that of MODIS in snow-free conditions with RMSE of 0.0176 (0.0236), bias of 0.0017 (0.0034), and R of 0.90 (0.84). Although the RMSE was slightly higher than the required accuracy of broadband LSA in the climate model (0.02-0.05) for snow-covered conditions, it was reasonable when its high value was considered. In addition, AHI broadband LSA shows good spatial and temporal consistency with MODIS broadband LSA. In most AHI full-disk domains, except for some pixels that were frequently covered by cloud, AHI broadband LSA shows a difference lower than 0.05 for both black-sky and white-sky, and AHI broadband LSA expressed diurnal variation in broadband LSA well over one year for five land types. The results of validation with ground measurements from 27 AERONET sites show that our results showed slightly higher RMSE (0.0223) than the MODIS RMSE (0.0195), and a weak positive bias is indicated. Although AERONET downward and upward fluxes used in this study were well validated [51], these are not direct measurement and most sites are located in heterogeneous areas. Thus, validation with direct measurements from other flux networks is needed for a more reliable quality assessment. AHI provides more valid LSA data than MODIS due to its high temporal resolution and despite of its relatively short BRDF synthesis period. Moreover, the mode value of the difference was close to 0. This advantage will be increased when AHI data observed every 10 min is used. When the TC technique was used, the accuracy was slightly improved, but no significant improvement was observed. This is similar to the results presented by Wu et al. [80]. AHI broadband LSA reasonably captured diurnal variation in ground measurements from the CRK site and showed improved temporal completeness in particular when clouds were frequently observed.
Our retrieval algorithm performs well to provide high-quality LSA with high spatial and temporal completeness from geostationary satellite observations in the Asia-Oceania region. We believe that this algorithm will be applied to GK-2A, the Republic of Korea's next-generation meteorological satellite. LSA data with high spatial and temporal resolution from the proposed algorithm can be utilized in various fields, such as validation of model results [11], simulation of weather and climate [86], and radiation budget studies [87]. Nevertheless, further research is necessary to improve and ensure the quality of LSA derived from the proposed algorithm. Further study is required to assess the stability of the LSA using long-term data and to validate LSA for more extensive spatial and temporal coverage, including snow-covered conditions to meet the validation stage 3 presented by the Committee on Earth Observation Satellites Land Product Validation subgroup.