{'fraction_non_alphanumeric': 0.033306412534324016, 'fraction_numerical': 0.01586173477628816, 'mean_word_length': 4.680429397192403, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Improved Cloud Gap-Filling Method for Longwave Infrared Land Surface Temperatures through Introducing Passive Microwave Techniques

: Satellite-derived land surface temperature (LST) data are most commonly observed in the longwave infrared (LWIR) spectral region. However, such data suffer frequent gaps in coverage caused by cloud cover. Filling these ‘cloud gaps’ usually relies on statistical re-constructions using proximal clear sky LST pixels, whilst this is often a poor surrogate for shadowed LSTs insulated under cloud. Another solution is to rely on passive microwave (PM) LST data that are largely unimpeded by cloud cover impacts, the quality of which, however, is limited by the very coarse spatial resolution typical of PM signals. Here, we combine aspects of these two approaches to ﬁll cloud gaps in the LWIR-derived LST record, using Kenya (East Africa) as our study area. The proposed “cloud gap-ﬁlling” approach increases the coverage of daily Aqua MODIS LST data over Kenya from <50% to >90%. Evaluations were made against the in situ and SEVIRI-derived LST data respectively, revealing root mean square errors (RMSEs) of 2.6 K and 3.6 K for the proposed method by mid-day, compared with RMSEs of 4.3 K and 6.7 K for the conventional proximal-pixel-based statistical re-construction method. We also ﬁnd that such accuracy improvements become increasingly apparent when the total cloud cover residence time increases in the morning-to-noon time frame. At mid-night, cloud gap-ﬁlling performance is also better for the proposed method, though the RMSE improvement is far smaller (<0.3 K) than in the mid-day period. The results indicate that our proposed two-step cloud gap-ﬁlling method can improve upon performances achieved by conventional methods for cloud gap-ﬁlling and has the potential to be scaled up to provide data at continental or global scales as it does not rely on locality-speciﬁc knowledge or datasets.


Introduction
Land surface temperature (LST) is classified as a Global Climate Observing System Essential Climate Variable (GCOS-ECV) [1]. It is widely used in both research and commercial applications, with its key domains of relevance including agriculture [2], urban landscape management [3], disaster risk analysis [4], and investigations on heat flux and hydrological features across the globe [5,6]. However, satellite LST data derived from brightness temperatures (BT) recorded in the longwave infrared (LWIR) spectral region have the key limitation that widespread cloud can completely obscure land surface from the sensors' view [7][8][9]. The resulting spatiotemporal 'cloud gaps' pose a challenge for applications reliant on regular and routine LST observations, both in terms of environmental models and derived data products [10,11]. The cloud gap problem may be further exacerbated when the sensor providing BT data is mounted on a polar orbiting satellite. In this scenario only a few (usually 1-2) views of an area in the tropical low latitudes and the temperate middle latitudes can be provided in every diurnal cycle, and these views can often be compromised by high (>50%) levels of cloud cover. This is a critical problem for both the monitoring of LST and its downstream users in applications, such as drought monitoring [12], urban island heat characterization [13], and fire detection [14,15].

Prior Cloud Gap-Filling Research
The filling of 'cloud gaps' in LWIR-derived LST data has been the subject of much research, with methods focusing primarily on data interpolation approaches based on temporal techniques such as the Savitzky-Golay filter [16,17], single spectrum analysis [18], and other similar methodologies [19]. The key assumption in such studies is a high continuity between daily LST dynamics at the scale of a remotely sensed pixel. However, daily continuity breaks down when a region suffers significantly from cloud shadowing, storm cells, differential urban heating, short-term precipitation events, and other environmental phenomena [20,21]. Therefore, there is an urgent need to enable a cloud gap-filling method that can account for these sub-daily changes in LST to minimize their negative effect on the intra-day continuity that temporal interpolation relies on. Geo-statistical interpolation is a variant of the above interpolation approaches for LST cloud gaps, in which areas of no-data are filled using functions built from the spatial distribution of the available clear-sky observations obtained at, or close to, the same local time at which the cloud gaps are observed [22][23][24]. However, this approach has been proven unreliable when the study area is affected by larger portions of contiguous data gaps [25] and is thus not suitable for use in areas suffering frequent cloud.
Another type of approach exploited to fill cloud gaps is based on LST variation and distribution patterns in the spatio-temporal space, sometimes referred to as "spatio-temporal data fusion" (STDF) methods. An STDF framework is usually constructed using clear-sky LST observations of spatially neighbouring pixels made at proximal or near-proximal dates. The methods deployed include linear models [26,27], Markov-based models [28], and those based on ratios between daily LST observations and the corresponding eight-day averages [29]. Improved performance of these methods, and thus more stable cloud gap LST estimates, have been obtained through the introduction of additional data beyond LST alone [25,30,31]. However, a flaw of this approach is highlighted by Song et al. [30] in that gap-filled pixels of MODIS LST data based on STDF methods represent an estimate of the LST that would be found under clear sky conditions, even if the gaps are caused by clouds. This is because the STDF approaches rely on relationships built from 'clear sky' LST observations. Therefore, a 'clear sky bias' within the STDF cloud gap-filled product can be introduced because the amount of outgoing longwave and incoming shortwave radiation in clouded areas is affected by the cloud cover in a manner not accounted for by the clear sky STDF relationship.
The use of physically based land surface energy balance (LSEB) modelling [31][32][33] and numerical weather prediction (NWP) modelling [34] are theoretically not affected by a clear-sky bias. However, since they require accurate local soil properties to work effectively, the LSEB models have large uncertainty on the model parameterization process. A recent investigation [35] has attempted to overcome this flaw by fusing LSEB theory with machine-learning-based methods. Despite the good accuracies (RMSE < 2.5 K) achieved in terms of results, their inputs still require spatially and temporally contiguous shortwave radiation products that are difficult to obtain globally with a high accuracy and at a high resolution from geo-stationary satellites [36]. In addition, validation of energy balance models of surface temperature has shown that there is a tendency for such models to suffer Remote Sens. 2021, 13, 3522 3 of 26 a cold bias [37,38]. The NWP modelling technique is much more sophisticated in describing the land surface and atmospheric conditions but has been shown to struggle in accuracy with regards to surface temperatures over a variety of landcovers [39,40]. Overall, physical models as they currently exist should be further improved to act as a suitable answer on their own to the filling of LST cloud gaps, especially on large (global) scales.
The final methodological approach available to tackle the cloud gap problem is to exploit passive microwave (PM) data. The advantage of PM-based approaches is the ability of microwave radiation to penetrate clouds, which allows the spaceborne sensor to record surface-emitted PM radiation and support the production of 'all-weather' LST data [41]. However, the ground pixel footprints of satellite-based PM instruments suitable for LST estimation are typically coarser than 10 × 10 km 2 . This is a far worse spatial resolution than is provided by LWIR sensors, such as the Moderate Spatial Resolution Imaging Spectroradiometer (MODIS) and the Visible infrared Imaging Radiometer (VIIRS), and insufficient for many applications of LST data given the strong spatial heterogeneity typical of land surface temperatures. Furthermore, PM-derived LST estimates represent the temperature at some depth in the soil, typically 2 mm to 1 cm, which is different to the surface skin temperature estimates typically provided by LWIR-based LST methods [42][43][44].
Due to the known limitations of PM-derived LST data, numerous investigations have been performed with respect to how PM-and LWIR-derived LST data records might be combined to produce an 'all-weather' daily LST record at a spatial resolution the same or similar to that of the original LWIR-derived LST record. Table 1 summarises key studies in this field from the last five years, almost all of which rely on the environmental characteristics of the particular region studied. Most of the approaches did not address the potential universality of their method beyond the particular study region, which is likely a limitation to producing the all-weather daily LST at continental/global scales. Table 1. Summary of LWIR and Passive Microwave data fusion studies related to the production of 'all-weather' daily LST information at a higher spatial resolution than PM data alone can produce.

Study
Main Method Description Spatial Scale Limitations

Kou et al. [45]
A Bayesian Maximum Entropy (BME) blending approach to merge PM and LWIR data, achieving a Root Mean Square Error (RMSE) in LST of between 2.3-4.5 K.
A relatively small 100 × 100 km 2 region Only used on night-time data over a small region, with its universality requiring more validation.
2. Duan et al. [43] An empirical model based on a digital elevation model (DEM) and clear sky LWIR-derived LST at neighboring pixels to downscale PM-derived data for achieving LST of cloudy LWIR pixels.

China
Downscaling of PM LST only relies on DEM, may be theoretically less effective in areas where the topographical variation is not important (e.g., low altitude plains).

Sun et al. [46]
A downscaling method for PM LST using NDVI and DEM data applied for gap-filling of LWIR LST at a finer resolution.

China
(1) Penetration depth difference between PM and LWIR not considered.
(2) The method was applied over a large area but only validated at limited sites.
(3) PM data were downscaled to 5 km but the feasibility of the method at a finer resolution (e.g., 1 km) requires further investigation. (1) This time-interpolation-like fusion method has strict requirement on the availability of its input datasets at neighboring dates [49].
(2) The method only addresses situations of temporally discontinuous pixel loss across relatively small study regions.
South Korea, about 10,000 km 2 ; Cauvery river basin in India, about 80,000 km 2 The physics behind machine learning models remains unclear, and it is difficult to justify the global universality of these models when the relationships they derive cannot be explicitly formalized.

Research Objective
Based on the research described above, both conventional STDF approaches and PM-LWIR data fusion approaches have a mixture of merits and limitations. Therefore, a combination of these has potential for compensating for the 'clear sky bias' present in the STDF approaches and for removing the regional restriction present in many existing PM-LWIR fusion methods. In this study, we build on past efforts by developing a novel two-step framework for better estimation of LST in cloud covered areas at a relatively high spatial resolution, based on a combination of LWIR-derived LSTs and PM LST data. We applied the approach at a national scale in tropical Africa, and through our efforts aim to further enhance the universality of an approach to generate all-weather, fine spatial resolution (~1 km) LST data. We also aim to investigate the effect of cloud insulation on the gap filling of LST, answering whether and how PM observations might mitigate biases in the cloud gap filled LST data generated by the conventional STDF methodology. Our work may benefit future efforts to produce globally quasi-full coverage daily LST data at resolutions similar to those provided by current LWIR sensors.

Study Area and In Situ LST Ground Observations
The study area for this work covers Kenya (East Africa) with a total area of around 582,646 km 2 . As illustrated in Figure 1, the far southwestern part of the country is characterized by highly vegetated cropland and forest, with complex topography at higher elevations (>1000 m). In contrast, the rest of the region is primarily dominated by flat plains and covered with sparser vegetation types (grass or shrub savannas). Most of Kenya is dominated by Arid (BSh, BWh) and Tropical (Af, Am, Aw) Köppen-Geiger climate zones with some Temperate (Cfb, Csb, Cwb) climate zones found in the south-west [53]. Snow and ice rarely occur except at high altitude. The local climate experiences a bi-modal rainfall pattern, with a season of 'short rains' typically between October and December and 'long rains' from April to late May or June. Other months are identified as dry seasons. Given the wide range of precipitation, topography and climate types, methodologies for improving LST cloud gap-filling developed in Kenya are likely applicable to many other areas globally, especially in the (sub) Tropics and the Temperate zones.
The ILRI Kapiti Research Station is located in south-west Kenya (−1.6083 • , 37.1327 • ), as shown in Figure 1a. This station is run by the International Livestock Research Institute (ILRI) and provides observations of LST and land surface energy balance, as well as CO 2 and H 2 O fluxes. The observation sites are instrumented using a series of mast mounted (Heitronics) LWIR radiometers to deliver upward (53 • view zenith angle) and downward (15 • view zenith angle) looking LWIR brightness temperatures. They are then processed to deliver logged LST estimates every 5-min according to the LST in situ measurement principles described in Gottsche et al. [54]. There are four sites in all (Masts 1-4 in Figure 1a) which can provide LWIR radiometer measurements suitable for the evaluation of satellitebased LST products since September 2018. The research station has a fenced 15 km × 20 km area of homogenous savannah, which can further be broken down into 94% grassland-soil complexes and 6% tree canopy cover. seen from the NDVI time series of five typical LWIR pixels in different locations a ferent land covers across Kenya in Figure 1b. Notice that the locations of the five i.e., Pixels A (at 3.63° N, 40.35° E), B (4.16° N, 35.23° E), C (0.78° S, 35.22° E), D (0.76° S, 4 and E (0.93° N, 39.37° E), have been marked in Figure 1c. They are respectively cove different land use types of shrubland, grassland, cropland, forest, and urban area. We employ the in situ LST data on Masts 1-4 of ILRI Kapiti Research Station as our primary benchmark (validation) dataset, focusing on the period from October 2018 to February 2019. A brief description of the derivation and accuracy evaluation of our benchmark LST data is provided in the Supplementary Material ( Figure S1). The certainty (range of possible error) of the ground observations ranges between 0.79 and 0.95 K per observation throughout the whole time series, across the four masts. The range of possible error consists of radiometric measurement uncertainty, surface emissivity uncertainty, and the uncertainty contribution of the logging equipment as driven by environmental temperature. This time period, encompassing both the short-rains and subsequent dry period, is sufficient in capturing the major climatic and vegetative cycles of an entire year (Figure 1b), by which land surface emissivity and the thermodynamic temperature are strongly affected [55]. In particular, the selected period includes the very rapid green-up (typically only a few days) when the rains start, the rainy sub-period accompanied by increasing vegetation canopy cover, and the sub-period of landscape drying when the rain ends. This phenological pattern is representative of the study area as a whole, as can be seen from the NDVI time series of five typical LWIR pixels in different locations and different land covers across Kenya in Figure 1b Figure 1c. They are respectively covered by different land use types of shrubland, grassland, cropland, forest, and urban area.

Satellite Datasets
Several satellite datasets were utilized for our study. The enhanced Aqua MODIS 1-km daily LST product (MYD21A1D/N.v061) was selected as the cloud gap-filling target. In addition, we also employed the MODIS MCD43A4 (v061) 500-m daily "Bidirectional Reflectance Distribution Function (BRDF)" adjusted reflectance dataset and the 90- In addition to the~5-month in situ data record coming from the ILRI Kapiti Research Station, we also employed the LST data from the geostationary Meteosat platform's Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) [56] as a second benchmark validation dataset. SEVIRI-derived LST data covering an entire climatic year (from October 2018 to September 2019) were used across the whole Kenya. The SEVIRI LST product is provided every 15 min by EUMETSAT's Land Surface Analysis Satellite Application Facility (LSA-SAF). Trigo et al. [56] describe the adapted split window algorithm for producing LST data from SEVIRI-measured LWIR brightness temperatures (BTs). Gottsche et al. [54] report that the product has an overall RMSE of 3.2 K compared to LWIR-radiometer derived LSTs collected at a West African validation station located in a tiger bush savanna biome. With a pixel size of around 4 km in Kenya, the primary role of Meteosat LST data was to aid evaluation of both the standard MODIS LST data and the cloud gap-filled MODIS LST data. The latter is possible because instances exist where MODIS pixels classed as cloudy (and thus with no MODIS LST data) are matched by Meteosat LST data that do have an LST observation classed as cloud free within a time lag of up to around 8 min, due to either cloud movements, differences in the cloud masking procedures, and/or due to the different viewing geometries of the two sensors (SEVIRI =~40-42 • view zenith angle (VZA) over Kenya, MODIS VZA = 0-65 • depending on overpass). This can be observed in the contrast show in Figure 4c,e in Section 4.

Methodology
The flowchart of our cloud gap-filling approach for MODIS LST data is schematically shown in Figure 2. It generally comprises two steps, namely (1) a conventional realization of the STDF methodology and (2) a PM-driven bias adjustment process to calibrate the data coming out of step (1), accounting for the fact that these data should represent cloud-covered rather than clear-sky LSTs. Therefore, there are three datasets tested in the validation step: (i) the clear sky original MODIS product (MODIS Clear ), (ii) MODIS with gaps filled by step (1)

Methodology Description
The 1 km MODIS LST data were subset to Kenya and pixels classed as cloud contaminated were removed using the "quality control (QC)" field contained within the MODIS LST product. In addition, all pixels flagged as 'cloud-edge' or low quality were also removed. A 500-m daily NDVI dataset for Kenya was calculated using the red-band and infrared-band reflectance data from the MODIS MCD43A4 product, which has been corrected for angular reflectance effects based on the BRDF. The resulting estimates were resampled using a spatial averaging procedure to the same 1-km resolution as the MODIS LST product. The spatial resampling process was also applied to the SRTM DEM dataset subset over Kenya. For each date (represented as t1) in the study period (from October  The 1 km MODIS LST data were subset to Kenya and pixels classed as cloud contaminated were removed using the "quality control (QC)" field contained within the MODIS LST product. In addition, all pixels flagged as 'cloud-edge' or low quality were also removed. A 500-m daily NDVI dataset for Kenya was calculated using the red-band and infrared-band reflectance data from the MODIS MCD43A4 product, which has been Remote Sens. 2021, 13, 3522 8 of 26 corrected for angular reflectance effects based on the BRDF. The resulting estimates were resampled using a spatial averaging procedure to the same 1-km resolution as the MODIS LST product. The spatial resampling process was also applied to the SRTM DEM dataset subset over Kenya. For each date (represented as t 1 ) in the study period (from October 2018 to September 2019), we used the STDF methodology version detailed in Song et al. [30] to fill cloudy pixels that have no LST data record within the MODIS daily LST product. The Song et al. [30] method was chosen as it has been applied to a large area in South China and obtained relatively good validation performance (with RMSEs between 1.5 and 3.5 K). A summary of the STDF method is given here, whilst the full details can be found in the work of Song et al. [30].
The STDF method first builds a transfer function based on the above-described datasets (LST, NDVI and DEM) as follows: where the superscript "*" indicates that this variable has been normalized to the range 0 to 1.0, based on the maximum and minimum values of that variable found across Kenya during our study period. Parameters a, b, c, and d are coefficients fitted using all pixels that have a clear-sky LST estimate on a specific date t 1 (LST * t 1 ) and at least one closely preceding date, t 0 (LST * t 0 ) (see below for further detail on how this was calculated). NDVI * t 1 indicates the corresponding (normalized) NDVI on the t 1 date. The NDVI data are available for almost every pixel continuously because the MCD43A4 daily product is generated using 16-day composited reflectance data, and therefore has good resistance to data loss from cloud contamination. Notice that day-time and night-time LST data must be separately treated in this transfer function. After deriving the coefficients of a, b, c, and d, the Equation (1) transfer function was used to fill all cloudy MODIS LST pixels on the t 1 date. For any t 1 date included in the study period, the t 0 date was iterated among all neighbouring dates of t 1 meeting the condition |t 0 − t 1 | <= 15 days (from the nearest date to the furthest date). An average of the estimated LST values for t 0 was then taken where a cloud gap pixel was filled more than once (based on the iterative t 0 dates), and the iteration was stopped when the fraction of pixels with effective LST values on t 1 was equal to or exceeded 0.9. In the original STDF method of Song et al. [30], the remaining small fraction (<0.1) of cloud gaps was then filled with LST data generated by a conventional inverse distance weighted (IDW) interpolation method [57]. However, the IDW step is not done here in order to avoid bringing in extra interpolation uncertainty that may negatively influence the validation step. The final modification to the original Song et al. [30] method is the use of iteration over all possible dates within a 30-day window (rather than only selecting one t 0 date that has a high fraction of clear-sky pixels) to increase the amount of data available to the procedure and therefore improve the applicability to the entire area of Kenya.

Method Sensitivity Analysis
In this section, we conducted sensitivity analysis for the STDF method expressed through Equation (1), in order to explore the existing uncertainty sources of the method which may influence performance of the subsequent bias adjustment process. Considering the linearity of Equation (1) and the fact that there are no interactive components between different input variables, the sensitivity of the output variable LST t 1 to each input variable X (denoting LST t 0 , NDVI t 1 , or DEM) can be expressed as the 1st order partial derivative to . It should be noted that these input variables are not normalized between 0-1. The sensitivity was calculated for all iterations conducted during Section 3.1.1. The statistical summary of this sensitivity analysis of the different iterations are reported in Table 2. From Table 2, it is apparent that the difference between daytime and night-time for sensitivity of output LST t 1 against LST t 0 and DEM is low. An LST t 0 difference of 1 K can cause a difference on average of 0.5-0.6 K for the (output) gap filled LST. In comparison an average decrease of about 0.3 K or 0.4 K was observed when elevation rises by 100 m. For NDVI, sensitivity is more important in the daytime compared with night-time, since an NDVI increase of 0.1 can result a decrease >1 K in the gap-filled outcome, while it has an almost negligible effect on the night-time outcomes.

Step 2: Cloud Shadowing Bias Adjustment with PM-Derived LST Data
As the STDF methodology is based on cloud-free LST observation data from MODIS but is attempting to estimate LST under cloud, biases may be introduced. Identifying and removing such bias is the focus of Step 2 and is based on the AMSR-2 PM-derived LST data. AMSR-2 provides both horizontally (H-) and vertically (V-) polarized BT observations in seven different frequencies [58]. LST data have been estimated by different combinations of these data [59][60][61][62][63], including at 6.9 GHz, 10.7 GHz, 18.7 GHz, 23.8 GHz, 36.5 GHz, and 89 GHz. The linear LST retrieval algorithm based on the 36.5 GHz V-polarized BT [61,64] is among the most widely used, though Song et al. [65] found it inferior to that based on the 18.7 GHz and 23.8 GHz bands [63,66], especially for pixels with higher fractions of water surface (including standing water and dynamic water signals like high levels of soil moisture). To minimize the influence from increased soil moisture on LST retrieval in the rainy season of our study area, we employ the AMSR-2 V-and H-polarized 18.7 GHz and 23.8 GHz BTs to retrieve LST using the methodology of Jones et al. [63]. The resulting PM-derived LST estimates have a spatial resolution of 25 km after the AMSR-2 Level1R BT data were resampled according to the global Equal Area Scalable Earth (EASE) Grid [67] projection system.
The PM BTs at 18.7 GHz and 23.8 GHz have a theoretical temperature sensing depth of around 1 cm [68], as compared to the MODIS LWIR-derived LST skin surface depth of ≤1 mm [69]. To deal with this we first averaged the MODIS 1-km LST data of Kenya to the 25-km EASE Grid and compared those cells that have more than 600 samples (i.e., >95%) of clear-sky MODIS pixels to the AMSR-2 retrieved LST for the same region. As seen in Figure 3, both daytime (ascending) and night-time (descending) MODIS and AMSR-2 LSTs show strong correlations but depart from the 1:1 line. This result agrees with the observations of Song et al. [65], who demonstrated that, apart from in desert and snow/ice covered areas, the LST derived from these two different sensors have a strong linear correlation at the global scale. The AMSR-2 observations have a negative bias of around −0.7 K in the daytime and a positive bias of around 2.5 K in the night-time. This is consistent with the logic that compared to temperature of sub-surface soils, that of skin surface is usually higher at mid-noon but lower at midnight. Therefore, the bias in Figure 3 can be largely ascribed to the different sensing depths. To enable the AMSR-2 LST data to provide the best estimate of what MODIS would have measured, we fit an independent linear equation to the data of day-time and night-time modes respectively, as shown by the orange text in Figure 3. These achieved high coefficients of determination (R 2 ≥ 0.9), with RMSEs between the AMSR-2 LST data adjusted via these equations and the true MODIS LST (RMSE unbias ) of less than 1.5 K. The small value of RMSE unbias proves a high degree of consistency between PM (AMSR-2)-derived LST and MODIS LWIR LST, thereby confirming the feasibility of the penetration depth bias adjustment process. Next, we inverted the linear relations in Figure 3 and thus transformed the AMSR-2 based LST data (LSTAMSR-2) into the corresponding best estimate of MODIS-derived LST at the 25-km pixel scale (LSTMODIS_25km): where and are the linear relationship parameters found in Figure 3. This removes the penetration depth bias. Data at 1-km resolution can then be further calibrated to remove the cloud shadowing bias based on the relationship: where the number of clear-sky MODIS LST pixels (LSTclear_sky) within a 25-km pixel is N1, and the number of cloud gap-filled pixels (LSTcloud_gap) is N2. We then assume that all cloud gap-filled MODIS pixels within the 25-km AMSR-2 pixel have the same shadowing bias (ΔLST). Therefore, ΔLST for all cloud gap-filled pixels within a certain 25-km grid can be calculated via rearrangement of Equation (3) to: . RMSEunbias denotes the fitting errors illustrated in Figure 3 and can be identified as the threshold value for random (non-systematic) error between LST retrievals made at different spatial resolutions. When ∑ △ lies below this threshold the LST error is mainly driven by this Next, we inverted the linear relations in Figure 3 and thus transformed the AMSR-2 based LST data (LST AMSR-2 ) into the corresponding best estimate of MODIS-derived LST at the 25-km pixel scale (LST MODIS_25km ): where k 0 and m 0 are the linear relationship parameters found in Figure 3. This removes the penetration depth bias. Data at 1-km resolution can then be further calibrated to remove the cloud shadowing bias based on the relationship: where the number of clear-sky MODIS LST pixels (LST clear_sky ) within a 25-km pixel is N1, and the number of cloud gap-filled pixels (LST cloud_gap ) is N2. We then assume that all cloud gap-filled MODIS pixels within the 25-km AMSR-2 pixel have the same shadowing bias (∆LST). Therefore, ∆LST for all cloud gap-filled pixels within a certain 25-km grid can be calculated via rearrangement of Equation (3) to:  Figure 3 and can be identified as the threshold value for random (non-systematic) error between LST retrievals made at different spatial resolutions. When ∑ N2 j=1 ∆LST N1+N2 lies below this threshold the LST error is mainly driven by this random error. The influence of this random error on ∆LST can be exaggerated if (N1 + N2) >> N2. To mitigate such influences, we use (N1 + N2) to substitute N2 as the new denominator in Equation (5) where N1+N2 ≤ RMSE unbias . This means that when random error dominates the overall shadowing bias of the 25-km grid, the bias will be equally allocated to all (N1 + N2) MODIS pixels within the grid regardless of whether they are clear sky or cloudy, as in Equation (5): Finally, the STDF cloud gap filled LST values are calibrated by the application of ∆LST, as found by Equation (4) or (5) depending on the amount of cloud cover within an AMSR-2 25-km grid cell, to give MODIS PMBC .

Overview of Cloud Gap-Filling
Our two-step cloud gap-filling processes improved the spatial coverage of the daily daytime and night-time MODIS Clear data of Kenya from <50% to >90% for the cloud gapfilled images, as can be seen in Figure 4a. However, it is important to note that the re-visit cycle of AMSR-2 at low latitudes such as those of Kenya is normally greater than a day, so not all cloud gap-filled MODIS LST pixels were able to be bias-adjusted based on a near-simultaneous AMSR-2 observation. Figure 4b shows the proportion distribution of bias adjustment fraction (BAF) for gap-filled pixels of each date. It indicates that BAF on most (≥49%) of dates is between 0.25 and 0.5. Although not high, this rate is sufficient to provide a basis for comparison between the bias adjusted and non-bias adjusted cloud gap-filled LST data. From Figure 4c-h, we employed the day-time data of 22 October 2018 as an example to demonstrate the spatial inputs, outcomes, and corresponding reference data of our methodology. Spatial distribution maps of LST are shown respectively for the non-bias adjusted result in Figure 4f and for the bias adjusted result in Figure 4g. The map of their difference is demonstrated in Figure 4h, revealing that some gap-filled LST values in the middle and south of the country are lower when bias adjustment is included. A series of quantitative validation steps against the reference datasets are carried out in the following sections.

Validation against In Situ LST Ground Observations
We first compared MODIS Clear observations from clear-sky pixels, MODIS STDF from Step 1, and MODIS PMBC from Step 2, to the contemporaneously collocated in situ LST data recorded at the four masts of the ILRI Kapiti Research Station site ( Figure 5 and Table 3). As the local observation times of MODIS and AMSR-2 are reported to have a variation between ±15 min from 13:30 for different dates and different observation swaths, the corresponding in situ data were determined by averaging all in situ LST samples from 13:15 to 13:45 local time (UTC+3) in the day-time, and all in situ LST samples from 01:15 to 01:45 local time at night-time. Moreover, in situ LSTs from Mast 2 and Mast 3 were averaged rather than being used independently, because these masts are in the same 1 km MODIS pixel. The use of multiple masts across and within different pixels (along with the geometric illumination model upscaling method applied, see Section Supplementary Figure S1) means that surface heterogeneity is well accounted for in the landcover under observation. We also examined the gap-filling results against the in situ LST data from Mast 1, Masts 2 and 3 (averaged), and Mast 4 as a time series (Figures 6 and 7). For the time series demonstration, each blue point (MODIS PMBC ) is accompanied by a corresponding green point (MODIS STDF ). Bias adjustments herein are made only on dates where LST estimates are available from both MODIS and AMSR-2. Bias-adjusted outcomes on dates when AMSR-2 data were not available are not shown. Table 3. Statistics resulting from the LST comparison shown in Figure 5. ** indicates significance at p < 0.05.   Table 3. The clearest benefit of the bias-adjusted gap filling is seen in the coolest temperatures of the day-time record.  Table 3. The clearest benefit of the bias-adjusted gap filling is seen in the coolest temperatures of the day-time record.  Table 3. The clearest benefit of the bias-adjusted gap filling is seen in the coolest temperatures of the day-time record.

Validation against SEVIRI Geostationary LST
Our cloud gap-filled MODIS LST data record was also validated against the LSA-SAF Meteosat SEVIRI LST dataset, based on the mean LST of all MODIS pixels within a SEVIRI pixel. This was carried out for the full twelve months of the year. Similar to the At night, all three datasets show a high degree of agreements with the in situ LST data, with absolute values of mean bias less than 1 K, and RMSEs below 1.5 K. The day-time performance of our cloud gap-filling approach appears poorer than that of night-time, this is likely due to the stronger spatial and temporal variations of LST at noon than at mid-night. However, MODIS PMBC agrees significantly better with the in situ LST data record in the daytime than does MODIS STDF . The former has an RMSE of 2.6 K compared to 4.3 K for the non-bias adjusted data, and this is similar to the 2.7 K RMSE shown by MODIS Clear . Based on all analyses above, the MODIS STDF values tend to overestimate the in-situ data in general, whilst MODIS PMBC values are in closer agreement in the daytime. At night, most of the cloud gap-filled satellite estimates are close to the in-situ values, regardless of whether they have been bias-adjusted or not.

Validation against SEVIRI Geostationary LST
Our cloud gap-filled MODIS LST data record was also validated against the LSA-SAF Meteosat SEVIRI LST dataset, based on the mean LST of all MODIS pixels within a SEVIRI pixel. This was carried out for the full twelve months of the year. Similar to the evaluation employed against the in situ LST data, SEVIRI data were compared in three independent scenarios. These scenarios are: (1) Figure 8, it is clear that daytime results differ more between clear-sky observations and cloud gap-filled estimates than do the night-time results. In the daytime, the best correspondence between MODIS and SEVIRI LSTs is achieved under clear-sky conditions, with a negligible bias and an RMSE of 3.2 K. The performance of the bias-adjusted cloud gap-filled result is slightly worse, with an RMSE of 3.6 K. However, this is an improvement upon the performance of the non-bias adjusted cloud gap-filled data, where an RMSE of 6.7 K and a positive LST bias of about 2.6 K were observed. At night, the performance of MODIS STDF and MODIS PMBC correction is similar, with RMSE no higher than 3.0 K.
Based on the yearly Kenya-wide comparison in Figure 8, further detailed analyses were conducted with respect to independent land cover types and daily timescales, as shown in Figures 9 and 10 respectively. The ESA CCI 300-m land cover map in 2018 (see Figure 1c) was exploited to provide land cover type information for the entirety of Kenya. The areal fraction of each land cover type in Figure 1c (except for water bodies) was calculated for all SEVIRI-viewed pixels (~4 km) within the country. A SEVIRI 4-km pixel containing any given land cover type with area fraction lower than 0.7 is identified as a 'mixed pixel' and thus eliminated, while the remaining 'pure' pixels were employed to evaluate the potential influence of land cover variation on the gap-filled results, as shown in Figure 9. The similar performance among different land covers was consistent with the overall RMSE results reflected in Figure 8. This suggests that the bias-adjusted cloud gap-filling method is effective over a wide range of vegetation covers in obtaining LST estimates with competitive accuracy against MODIS clear-sky LST observations. The poorest performance of MODIS PMBC is found over bare soil/rock during the day (RMSE ≈ 4.2 K) and over sparse vegetation at night (RMSE ≈ 4.0 K). However, this is acceptable considering that the difference between MODIS pmbc and MODIS clear RMSE is insignificant (<1.5 K). Based on the yearly Kenya-wide comparison in Figure 8, further detailed analyses were conducted with respect to independent land cover types and daily timescales, as shown in Figures 9 and 10 respectively. The ESA CCI 300-m land cover map in 2018 (see Figure 1c) was exploited to provide land cover type information for the entirety of Kenya. The areal fraction of each land cover type in Figure 1c (except for water bodies) was calculated for all SEVIRI-viewed pixels (~4 km) within the country. A SEVIRI 4-km pixel containing any given land cover type with area fraction lower than 0.7 is identified as a 'mixed pixel' and thus eliminated, while the remaining 'pure' pixels were employed to evaluate the potential influence of land cover variation on the gap-filled results, as shown in Figure 9. The similar performance among different land covers was consistent with the overall RMSE results reflected in Figure 8. This suggests that the bias-adjusted cloud gapfilling method is effective over a wide range of vegetation covers in obtaining LST estimates with competitive accuracy against MODIS clear-sky LST observations. The poorest performance of MODISPMBC is found over bare soil/rock during the day (RMSE ≈ 4.2 K) and over sparse vegetation at night (RMSE ≈ 4.0 K). However, this is acceptable considering that the difference between MODISpmbc and MODISclear RMSE is insignificant (<1.5 k).
The daily time series of difference between RMSE for gap-filled LST (MODISSTDF and MODISPMBC) and RMSE for clear sky LST (MODISClear) are demonstrated in Figure 10. The fraction of clear sky LST pixels over the study area for each day are also shown using the red bar on the right-hand y-axis. From the contrast between blue (MODISPMBC) and orange (MODISSTDF) lines, we can see that the relative performance between the bias adjusted gapfilling algorithm against the non-bias adjusted algorithm is stable over different seasons of the year. The performances of the non-bias adjusted vs. bias adjusted algorithms are not significantly influenced by the fraction of clear sky LST pixels that are needed to build the gap-filling model.     The daily time series of difference between RMSE for gap-filled LST (MODIS STDF and MODIS PMBC ) and RMSE for clear sky LST (MODIS Clear ) are demonstrated in Figure 10. The fraction of clear sky LST pixels over the study area for each day are also shown using the red bar on the right-hand y-axis. From the contrast between blue (MODIS PMBC ) and orange (MODIS STDF ) lines, we can see that the relative performance between the bias adjusted gapfilling algorithm against the non-bias adjusted algorithm is stable over different seasons of the year. The performances of the non-bias adjusted vs. bias adjusted algorithms are not significantly influenced by the fraction of clear sky LST pixels that are needed to build the gap-filling model.

Influence of Cloud Duration on PM-Based Calibration
As the PM-based bias adjustment was found to have more of an effect in the daytime than at night-time, we further investigated the impact of cloud duration during the day and night on bias correction performance. Firstly, data from the four selected dates highlighted in Figure 6 were employed to calculate an evaluation metric, ∆bias. This metric is defined as the difference between the absolute bias of the MODIS STDF (|∆LST STDF |) and that of MODIS PMBC (|∆LST PMBC |). The ∆bias metric (|∆LST STDF | − |∆LST PMBC |) was computed over the pixels in which ground validation masts are located. It can be used as an indicator of the accuracy improvement found through applying the PM-based LST bias-adjustment approach of Step 2. A time series indicating the presence or absence of morning cloud was then generated (1: cloud, 0: clear sky) from 7:30 a.m. to 1:30 p.m. local time for each of the four dates selected. Cloud state was provided by data from the sky-pointing radiometer deployed on the masts at ILRI Kapiti Research Station (see Supplementary Material Section S2, Figure S2). The dates selected for this were chosen at random from the days within four equal divisi)ons of the available time period. Based on this time series, a daily cloud duration fraction (CDF) was calculated, defined as the number of stable (occluded sky for a minimum of 15 min prior to the given minute) cloud minutes from 7:00 a.m. to 1:30 p.m., divided by the total number of minutes in that period. The results are displayed in Figure 11 and suggest that there is a relationship to be explored between ∆bias and CDF, especially considering the low ∆bias of 1.3 K and a small CDF of 0.14 on 2nd November 2018 as compared to the far higher ∆bias (4.8 K) and correspondingly high CDF of 0.85 of the morning of 15 January 2019. Results in Figure 11 suggest that the performance of our bias adjustment methodology in the daytime is influenced by the length of the cloud period experienced during the morning prior to the early afternoon MODIS overpass whose cloudy observations were filled by our proposed method. To corroborate the relationship found in Figure 11, a further analysis was carried out using the Meteosat SEVIRI LST data, as this data contains cloud masking information that is available over the entirety of Kenya at 15-minute intervals. This provides a more extensive test than with the in-situ data alone. A CDF based on the SEVIRI cloud mask was calculated for both daytime and night-time periods, following the same method as detailed above. For the day-time data, the CDF period was slightly extended to between 7:00 a.m. and 2:00 p.m. in order to account for differences in the observation time of the satellite To corroborate the relationship found in Figure 11, a further analysis was carried out using the Meteosat SEVIRI LST data, as this data contains cloud masking information that is available over the entirety of Kenya at 15-min intervals. This provides a more extensive test than with the in-situ data alone. A CDF based on the SEVIRI cloud mask was calculated for both daytime and night-time periods, following the same method as detailed above. For the day-time data, the CDF period was slightly extended to between 7:00 a.m. and 2:00 p.m. in order to account for differences in the observation time of the satellite sensors. Similarly, for the night-time data the period between 7:00 p.m. of the previous date and 2:00 a.m.of the observation date was used. Due to the temporal resolution of SEVIRI (15 min) as compared to that of the in-situ data (5 min), only 28 unique CDF values can be obtained in each seven-hour period.
The outcome (Figure 12) is that in the daytime, the mean ∆bias rises as the CDF increases from 0 to 1, while the standard deviations are close for all CDF-based groups. This phenomenon, together with the results in Figure 11, indicate that better day-time performance of our bias adjustment strategy is achieved under conditions of higher CDF ranges in the morning, especially when CDF exceeds 0.5. Essentially, longer periods of cloud cover prior to the afternoon Aqua MODIS overpass tend to provide improved performance of the bias adjustment approach. This is logical because, up to a certain limit, the longer a pixel is covered by cloud, the longer it is shadowed from direct sunlight and therefore the larger the temperature difference will be relative to clear-sky conditions. By contrast, at night there is only a very slight increase in ∆bias with increasing CDF, probably indicating that the change in LST is more influenced by the presence/absence of solar radiation in the day time rather than by the presence/absence of downwelling atmospheric radiation from cloud layers at night. phenomenon, together with the results in Figure 11, indicate that better day-time performance of our bias adjustment strategy is achieved under conditions of higher CDF ranges in the morning, especially when CDF exceeds 0.5. Essentially, longer periods of cloud cover prior to the afternoon Aqua MODIS overpass tend to provide improved performance of the bias adjustment approach. This is logical because, up to a certain limit, the longer a pixel is covered by cloud, the longer it is shadowed from direct sunlight and therefore the larger the temperature difference will be relative to clear-sky conditions. By contrast, at night there is only a very slight increase in Δbias with increasing CDF, probably indicating that the change in LST is more influenced by the presence/absence of solar radiation in the day time rather than by the presence/absence of downwelling atmospheric radiation from cloud layers at night. . This analysis highlights that during the day, the STDF+PMBC approach results in improved performance over the STDF-only approach, and this improvement increases with cloud residence time, as indicated by CDF. At night, bias correction offers little additional performance improvement.

Improvement and Universality of the Gap-Filling Methodology
Overall, our evaluation of the cloud gap-filled MODIS LST data against the benchmark data, based on both the in situ and Meteosat SEVIRI LST data records, provides a consistent narrative. For day-time Aqua MODIS overpasses gap-filled LST data using the conventional STDF methodology shows a positive bias. This supports our supposition that, if used alone, the STDF approach overestimates LST due to simulating what the LST would be under clear-sky type conditions, whereas in fact the location is or has recently been cloudy. Our results also reveal the effectiveness of using PM observations for calibration of the gap-filled LWIR-derived LSTs, especially during the daytime.
The implementation of this improved two-step cloud gap-filling framework on LST . This analysis highlights that during the day, the STDF+PMBC approach results in improved performance over the STDF-only approach, and this improvement increases with cloud residence time, as indicated by CDF. At night, bias correction offers little additional performance improvement.

Improvement and Universality of the Gap-Filling Methodology
Overall, our evaluation of the cloud gap-filled MODIS LST data against the benchmark data, based on both the in situ and Meteosat SEVIRI LST data records, provides a consistent narrative. For day-time Aqua MODIS overpasses gap-filled LST data using the conventional STDF methodology shows a positive bias. This supports our supposition that, if used alone, the STDF approach overestimates LST due to simulating what the LST would be under clear-sky type conditions, whereas in fact the location is or has recently been cloudy. Our results also reveal the effectiveness of using PM observations for calibration of the gap-filled LWIR-derived LSTs, especially during the daytime.
The implementation of this improved two-step cloud gap-filling framework on LST only depends on satellite remote sensing datasets. Compared to the existing methods (Table 1), it is easier to be implemented and is independent from all auxiliary datasets (e.g., reanalysis data) other than satellite remote sensing inputs. Our approach requires much less dependence on the particular environmental characteristics of the study area than some prior works (Table 1), especially when we take into account Figure 10 which shows that the accuracy of this method is not strongly influenced by the varied daily area fraction of cloudy pixels. The performance is stable under periods of continuous cloudy/rainy weather, as can be seen from the time series (e.g., between late November and early December 2018) in Figure 7. Taken together, these tests indicate the greater universality of the proposed method beyond its current study area, as compared to other existing methods.
It is noticeable that we find a relatively inconsequential cloud insulation effect on LST bias correction at night (Figure 12b), as opposed to the effect seen in the morning-tonoon window. The difference between daytime and night-time observations has rarely been discussed in previous studies related to cloud gap-filling of an LWIR-derived LST dataset, with some authors ignoring the problem by only tackling night-time cloud gapfilling [45]. Other works may have failed to find significant day/night differences, because they calibrated against air temperatures rather than the more physically accurate surface temperatures [30,70]. Our findings suggest that the difference between bias-adjusted and non-bias adjusted gap-filled LST, driven by cloud insulation, is less sensitive to ground longwave outgoing radiation and downwelling atmospheric radiation at night, as compared with the shadowing of solar short-wave incoming radiation during the day.
From Section 2.1 we can see that Kenya is representative of a variety of different geographical settings due to the varied land covers, climate zones and topographical conditions found across the country. Based on the stable evaluation of results across the different land covers in Figure 9, we suggest that the improved two-step cloud gap-filling framework should be effective in the vast majority of vegetation covers ranging from sparsely vegetated soils and grassland to forest found in low-and middle-latitude regions, where most of the world's human population resides. However, our current method is likely to be unsuitable for use in several land cover types, namely desert and snow/ice covered areas. This is because the temperature difference between microwave observations and LWIR skin surface observations is very high in such landcovers, whilst a conventional linear calibration (as applied in Equation (2)) insufficiently accounts for the depth-induced temperature differences encountered in these areas [65].

Uncertainty and Limitations of the Current Study
It is important to stress that neither the benchmark datasets (field radiometer derived LSTs and the Meteosat LST product) or the MODIS Clear data employed in this study represent the 'absolute truth' in terms of surface temperature observations, and we do not consider them as such. Each sensor is subject to its own inherent measurement errors. Moreover, measurements observed by the different sensors may differ on account of the heterogenous spatial scales at which observations are made and aggregated over sensor footprint geometry, as well as changes in environmental conditions occurring between non-temporally contemporaneous measurements. As such, we focus here upon the relative difference in performance metrics between the gap filled data (MODIS STDF and MODIS PMBC ) and MODIS Clear , rather than the absolute performance of these gap filled products vs. the SEVIRI and field radiometer benchmark LST datasets themselves. For example, using the Meteosat LST product as a benchmark, a relatively large absolute daytime RMSE (3.6 K) was obtained for pixels gap-filled using the MODIS PMBC dataset. This value is however much closer to the absolute daytime RMSE value (3.2 K) obtained from the 'raw' MODIS Clear dataset than the RMSE value (6.7 K) obtained using the MODIS STDF dataset, highlighting the effectiveness of our novel cloud gap-filling two-step framework. While the absolute RMSE found between the MODIS Clear and the benchmark datasets during the daytime was observed to be high (>3 K), exploration of the causes of this discrepancy is beyond the scope of the current study. Furthermore, the validation effort of Gottsche et al. [54] found an average RMSE of 3.2 K for Meteosat SEVIRI data over a similar biome, indicating RMSE values of this magnitude are currently inherent in LWIR LST products in semi-arid savannah regions.
A limitation of the current study is the relatively small number of field-based LST observations available to validate our cloud gap-filling methodology. Radiometer measurements made at the ILRI Kapiti Research Station are temporally limited, spanning a~5-month period. Fortunately, however, the period is sufficient to capture one entire climatic and vegetation cycle of local land surface (Figure 1b). In the spatial domain, the validation sites only coincide with ground footprints of three distinct 1-km MODIS pixels that are in close geographic proximity to one another relative to the size of the study area. Nevertheless, the establishment and maintenance of even a single site represents a substantial achievement. Indeed, only one other satellite validation site capable of producing similar LST validation data exists anywhere on the African continent and it is not currently operational [54]. In addition, the limited availability of ground validation data has been compensated for to some extent by the intercomparison of the MODIS-derived products against the independent Meteosat SEVIRI LST product [56] at the national scale through an entire climatic year. The intercomparison between the MODIS derived products and the ground/Meteosat-SEVIRI validation data provided consistent results. Prior studies on cloud gap-filling have primarily relied on validation strategies that artificially remove data to create "vacant sub-regions" of the imagery to replicate the presence of cloud gaps but for which the 'true values' of the data are actually known [26,71]. However, this strategy is unsuited to evaluating cloud gap-filling of LST data records, since the pixels for which the "true values" of LST are known have actually been observed under clear sky conditions and have thus not experienced the radiation interception effects of clouds that would be the case for 'true' cloud-covered LST pixels.
The sensitivity analysis results reported in Table 2 indicate that a theoretical error as large as 100 m in the DEM only leads to a difference of 0.3-0.4 K in LST on average. As the vertical error of the SRTM DEM is reported to be <16 m [72], the uncertainty of elevation data can be effectively neglected. On the other hand, the influence generated by uncertainty in the NDVI data can be much more significant. This is because the generation of daily-scale NDVI data from the 16-day composite MCD43A4 product relies on an assumption that vegetation change is not significant during the 16-day window. Unfortunately, vegetation changes rapidly during the onset of the rainy seasons in Kenya, which may result in a not-insignificant estimation error in NDVI estimates at certain dates. The resultant gap filled LST data is sensitive to such errors, especially during the daytime (as described in Section 3.1.2). Consequently, the differing performance between MODIS STDF and MODIS PMBC in the daytime cannot only be attributed to the effects of cloud as performance can also be impacted by NDVI uncertainty, particularly during the onset of the rainy season when rapid vegetation changes occur. However, the impact of cloud effects can still be identified as the primary driver of the performance differences between MODIS STDF and MODIS PMBC , based on the positive correlation observed between the CDF and ∆bias of the day-time data in Figures 11 and 12. Furthermore, the impact of the bias correction did not get significantly smaller at the rainy periods (typically between April-June and between October-December) in Figure 10.
Finally, despite the potential bias due to clouds as well as the uncertainty of NDVI input not being very significant for the night-time observations in this study, more validation experiments are required in other study areas to test the wider applicability of this finding. Currently, we recommend applying the PM bias correction process to both daytime and night-time data in future studies, in the event that larger night-time errors can arise from factors that were not present in the current study.

Future Work
Several points remain where further improvements to our approach might be investigated: (1) Biases calculated at the coarse (25-km) resolution of AMSR-2 are linearly allocated to each cloud gap pixel at the fine (1-km) resolution in the current study, based on an assumption that cloud geometry (e.g., cloud top height, cloud thickness, etc.) is homogenous within the AMSR-2 pixel. However, a bias-adjustment approach that considers the influence of cloud geometry differences at the AMSR-2 sub-pixel scale may offer further performance improvements to the current two-step LST cloud gap-filling framework.
(2) If applying the methodology at continental to global scales, the derived relations between AMSR-2 and MODIS LSTs (Equation (2)) should be re-explored because they are likely dependent on seasonal features and climate types that will vary at the continental/global scales.
(3) As the re-visit cycle of AMSR-2 is greater than one day at low latitudes, better spatial coverage may be achieved by the fusion of AMSR-2 data and data from other PM radiometers.
(4) Cloud gap-filling of Terra MODIS LST data (MOD11A1) could be similarly conducted using data from a PM radiometer similar to AMSR-2, but with local time of overpasses around 10:30 a.m./p.m. One candidate is the microwave radiation imager onboard the Fengyun-3C satellite [73].

Conclusions
In this study, we have presented an effective two-step framework for improving the cloud gap-filling performance of the LWIR-based daily land surface temperature (LST) data delivered by Aqua MODIS. The primary innovation of the framework lies in the coupling of a conventional "spatial-temporal data fusion (STDF)" methodology to a cloud shadowing bias-adjustment process, based on passive microwave (PM) LST data derived from the AMSR-2 instrument. We have evaluated the resulting STDF-based cloud gapfilled and PM-bias adjusted LST data against both in situ LST data and the geostationary SEVIRI-derived LST data record. In the daytime, the STDF-based outcomes show an RMSE of 4.3 K against the in-situ data and 6.7 K against the SEVIRI LST data. The RMSEs of the PM bias-adjusted outcomes are 2.6 K and 3.6 K respectively. At night, the performance of bias-adjusted outcome and that of STDF-based outcome are very similar.
Overall, we conclude that, compared to the STDF approach alone, our two-step framework provides a means to improve the performance of cloud gap-filling on LWIRderived LST through fusion with PM data, although the degree of improvement is far more significant for daytime compared to night-time observations. Finally, we have shown that the day-time accuracy improvements of our two-step approach to cloud gap-filling are increasingly apparent (relative to conventional STDF methods) under the condition of increased cloud cover residence time in the morning-to-noon timeframe. The proposed method has been shown to be stable in terms of RMSE change through time, even under spatially and temporally continuous pixel loss conditions, across a wide range of land cover types. In the future, this method could be applied to the next generation of sub-1 km all-weather LST products at the continental or global scale.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/rs13173522/s1, Figure S1: Summary of the available in situ LST record (Kelvin) for Kapiti site configuration 1. Data is available when at least one ground observing radiometer is active and returning admissible data. The upscaled mean is only available when at least two radiometer observations are available across all masts. Site time series are the downwelling and emissivity corrected LST values for each mast upscaled to the SEVIRI pixel scale. The final (E) time series is the Kapiti mean LST derived from all four sites and upscaled to the SEVIRI pixel scale. Figure S2: Stable cloud cover analysis for Mast 3. 1 is stable cloud present. 0 is stable cloud not present.