Evaluation of Six High-Spatial Resolution Clear-Sky Surface Upward Longwave Radiation Estimation Methods with MODIS

Surface upward longwave radiation (SULR) is a critical component in the calculation of the Earth’s surface radiation budget. Multiple clear-sky SULR estimation methods have been developed for high-spatial resolution satellite observations. Here, we comprehensively evaluated six SULR estimation methods, including the temperature-emissivity physical methods with the input of the MYD11/MYD21 (TE-MYD11/TE-MYD21), the hybrid methods with top-of-atmosphere (TOA) linear/nonlinear/artificial neural network regressions (TOA-LIN/TOA-NLIN/TOA-ANN), and the hybrid method with bottom-of-atmosphere (BOA) linear regression (BOA-LIN). The recently released MYD21 product and the BOA-LIN—a newly developed method that considers the spatial heterogeneity of the atmosphere—is used initially to estimate SULR. In addition, the four hybrid methods were compared with simulated datasets. All the six methods were evaluated using the Moderate Resolution Imaging Spectroradiometer (MODIS) products and the Surface Radiation Budget Network (SURFRAD) in situ measurements. Simulation analysis shows that the BOA-LIN is the best one among four hybrid methods with accurate atmospheric profiles as input. Comparison of all the six methods shows that the TE-MYD21 performed the best, with a root mean square error (RMSE) and mean bias error (MBE) of 14.0 and −0.2 W/m2, respectively. The RMSE of BOA-LIN, TOA-NLIN, TOA-LIN, TOA-ANN, and TE-MYD11 are equal to 15.2, 16.1, 17.2, 21.2, and 18.5 W/m2, respectively. TE-MYD21 has a much better accuracy than the TE-MYD11 over barren surfaces (NDVI < 0.3) and a similar accuracy over non-barren surfaces (NDVI ≥ 0.3). BOA-LIN is more stable over varying water vapor conditions, compared to other hybrid methods. We conclude that this study provides a valuable reference for choosing the suitable estimation method in the SULR product generation.


Introduction
The Earth's surface radiation budget (SRB) is a driving factor in the exchange of energy between the atmosphere and oceans/land [1]. The SRB has a great influence on climate and land cover changes and is critical in ecological, hydrological and biogeochemical land surface processes [2]. The surface upward longwave radiation (SULR; the total surface upward radiative flux in the 4 to 100 μm spectral using TOA and nonlinear regression (TOA-NLIN), (5) the hybrid method using TOA and ANN (TOA-ANN), and (6) the new hybrid method using BOA and linear regression (BOA-LIN). Both TE-MYD21 and BOA-LIN are newly developed SULR methods presented here. The data we used for the evaluation is presented in Section 2. The details of the two physical temperature-emissivity methods and four hybrid methods are introduced in Section 3. The evaluation results and discussions are demonstrated in Section 4. Finally, conclusions are given in Section 5.

Data and Processing
We conducted the evaluation with two datasets: a simulation dataset for the hybrid method comparison and a dataset containing two years of MODIS products and corresponding SURFRAD in situ measurements for the evaluation of all six methods. For the evaluation and coefficients calculation of hybrid methods, representative databases of typical atmospheric profiles and surface spectral emissivities are required to generate a simulation dataset of SULR and corresponding target variables (e.g., TOA radiances or BOA radiances in this study). The construction of the simulation dataset is introduced in Section 2.1. In the evaluation using MODIS and in situ measurements, the MODIS products (including TOA radiances, viewing zenith angle (VZA), geolocation, LST, narrowband emissivity, and cloud information), atmospheric profiles (from ERA5), and in situ SULR measurement (from SURFRAD) are used. The data and corresponding processing steps are introduced in Section 2.2.

The Simulation Dataset for Hybrid Method
Representative databases of atmospheric profiles, LSTs and spectral LSEs are important to make the hybrid method fit real situations as much as possible. In this study we adopted the widely used Thermodynamic Initial Guess Retrieval (TIGR) database [23,24] database, which is composed of 2311 profiles measured from global atmospheric conditions. Based on the criterion that in all layers the relative humidity should be less than 90% [22], we selected 946 clear-sky profiles. The number of selected clear-sky atmospheric profiles for tropical, mid-latitude, and polar atmospheres is 236, 258, and 452, respectively. Their water vapor (wv) varies from 0.06 to 6.3 g/cm 2 [25,26]. The LSTs were set with a (−10 K, 15 K) offset to the bottom temperature of the selected atmospheric profiles in a 5 K step. The LSEs (ε(λ)) were extracted from the MODIS UCSB (University of California, Santa Barbara) spectral library [27] which contains 158 emissivity spectra in the spectral range of 3.3 to 14.6μm. To balance computation efficiency and surface representativeness, 35 typical emissivities were selected from different surface types, including three spectra for water, one spectrum for ice, one spectrum for snow, 13 spectra for soils and minerals, and 17 spectra for vegetation.
Postprocessing for the spectral LSEs is required because the spectral emissivity above 14 μm is not available in the MODIS UCSB library. Two extrapolation methods have been proposed to fill the gap between 14 μm and the upper wavelength boundary (i.e., 100 μm) for SULR estimation. Tang et al. [12] set the gap emissivity as unity. Wang et al. [28] developed an emissivity extrapolation method which is a linear combination of MODIS band 29, 31, and 32 emissivities (Equation (1)) and has been adopted in this paper.
where , , and are the band-effective emissivity for MODIS band 29, 31, and 32, where is the emissivity value in spectral range 14 to 25 μm calculated with narrow-band MODIS emissivity, and is the emissivity value in spectral range 25 to 100 μm, which is set with the value of .
The MODerate resolution atmospheric TRANsmission (MODTRAN) model was selected to calculate the spectral atmospheric downwelling radiance (L↓(λ)), spectral atmospheric upwelling radiance (L↑(λ)), and spectral atmospheric transmittance (τ(λ)) with the inputs of 946 clear-sky TIGR atmospheric profile [29]. The surface spectral thermal emission (B(λ,T)) were obtained by the Plank function for the given wavelength (λ) and LST (T). The surface spectral outgoing radiances (Ioutgoing(λ)) include the self-emitted and reflected atmospheric downwelling radiance as shown in Equation. (2). The SULR can be integrated with the Ioutgoing(λ) by Equation.

The Evaluation Datasets with MODIS, ERA5 and in situ Measurements
The flowchart of data processing for the MODIS, ERA5, and in situ measurements is shown in Figure 1. First, the MODIS clear-sky pixels were screened with the MODIS cloud product. Then, the clear-sky LSTs, LSEs, and the in situ measured DLR in SURFRAD stations were used to drive the temperature-emissivity methods (TE-MYD21 and TE-MYD11) and the clear-sky TOA radiances, and the coefficients of TOA methods were used to drive the TOA hybrid methods (TOA-LIN, TOA-NLIN, and TOA-ANN). The data processing for the BOA-LIN method includes four steps: (1) Extracting the atmospheric profiles form ERA5 global reanalysis product based on geolocation (latitude and longitude) and the Shuttle Radar Topography Mission (SRTM) digital evaluation model (DEM) [30]; (2) calculating the atmospheric parameters (τ, L↓ and L↑) using ERA5 atmospheric profiles; (3) performing the atmospheric correction from clear-sky TOA radiances to BOA radiances using the atmospheric parameters; (4) calculating the SULR using the BOA radiances with BOA coefficients. Finally, the six methods were evaluated using the SURFRAD SULR in situ measurements.

MODIS Datasets
The MODIS instruments onboard satellite Terra (Aqua) have been observing the Earth continuously since 1999 (2002). MOD21 product has not been produced beyond 2005 because it requires bands 29, 31, and 32 as inputs and the infrared bands 27-30 of Terra-MODIS have been negatively impacted by an optical crosstalk issue since 2005. Observations from 2017 to 2018 were considered here; therefore, only the MYD products generated from Aqua-MODIS are used. The Aqua-MODIS is a whisk-broom sensor with a cross track and along track swath of 2330 km x 10 km. Its maximum VZA at the image edge is >60 degrees. The sensor has 36 spectral bands located in the 0.405 to 14.385 μm range, with spatial resolutions of 1000 m in the TIR bands. MODIS currently has 44 data products, which can be grouped into 5 classes: calibration, atmosphere, land, cryosphere and ocean. All these data can be found and downloaded from NASA's EARTHDATA system (https://earthdata.nasa.gov).
The MODIS datasets of MYD021KM, MYD03, MYD35_L2, MYD11B1, MYD11_L2, MYD21_L2, MOD13A2, and MYD13A2 from MODIS collection 6 were used for SULR estimation and analysis in this study. The list of the MODIS products used in this study is summarized in Table 1. The MYD021KM product provides the MODIS/Aqua 5-min L1B calibrated radiances at 1 km swath data and the radiances of bands 29, 31, and 32 were selected for the SULR estimation using hybrid methods. The MYD03 provided the geolocation and viewing geometry information for the swath data. The MYD35_L2 product provided the 1 km resolution cloud mask and also ancillary parameters. Four clear-sky confidence levels (confident clear, probable clear, probable cloudy, and cloudy) are given and thin cirrus information is also detected and flagged in the MYD35_L2 cloud product. The pixels with confident/probable clear and not cirrus contaminated flags were treated as non-cloud pixels. Next, the selected pixels are further filtered to ensure that the 3x3 neighborhood pixels are non-cloud pixels. The MYD11 products (MYD11B1 and MYD11_L2) were used to drive one temperature-emissivity method (TE-MYD11), and the MYD21_L2 product was used to drive another temperature-emissivity method (TE-MYD21). The MOD13A2 and MYD13A2 products contain the 16-day composed vegetation index information (Normalized Difference Vegetation Index (NDVI)) and are used to analyze the performance of TE-MYD21 and TE-MYD11 over barren and non-barren surfaces. Two years (2017 and 2018) of MODIS product data (~1 TB and ~6000 scenes) were rigorously processed in this study and resulting in 3197 available clear sky pixels at last. The BOA-LIN method uses MODTRAN for atmospheric correction requiring additional atmospheric profile information for MODIS data processing. Li et al. [31] evaluated National Center for Environmental Prediction (NCEP) operational global analysis data and MODIS atmospheric products (MOD07) for single channel LST retrieval with ground measurements and the result showed that the NCEP performed better than the MOD07 with an RMSE of 1. ERA5 hourly data on pressure levels is the latest version of global reanalysis data released by ECMWF on June, 2018 that aimed to replace the ERA-Interim with better spatial and temporal resolution. The ERA5 data has better spatial-temporal resolution (0.25° x 0.25° and 1 hour) than all the eight reanalysis products introduced above (no better than 0.5° x 0.625° and 3 hours). Therefore, the ERA5 data is selected to perform atmospheric correction in this study.
The ERA5 hourly pressure levels data contains 16 variables at 37 fixed pressure levels (from 1 mb to 1000 mb). The variables of geopotential height, pressure, temperature, and relative humidity are extracted and temporally interpolated to the overpass time of MODIS, vertically interpolated to the elevation of the station based on 30m SRTM DEM dataset. The interpolated atmospheric profiles are then used with MODTRAN 4.3 to calculate the atmospheric transmittance and upwelling radiance for the atmospheric correction of MODIS TOA radiances.

Ground Measurements
The SURFRAD is a U.S. national-scale network that measures the SRB in a continuous mode beginning with four stations in 1995, then two more stations were added in 1998 and the latest station was installed in 2003 [33]. The SURFRAD stations are maintained by the National Oceanic and Atmospheric Administration (NOAA), and the data can be downloaded anonymously from ftp://aftp.cmdl.noaa.gov/data/radiation/surfrad/. The information of the seven SURFRAD sites is as shown in Table 2. The in situ measured SULR and DLR in the seven SURFRAD sites over the two years of interest (2017 and 2018) were used to evaluate the six methods and to drive the temperature-emissivity methods, respectively. The SULR and DLR in SURFRAD were measured with pyrgeometers which were mounted on a tower ~8 m above the ground. The pyrgeometer has an accuracy of about 1% and is calibrated annually with carefully selected calibration devices [34]. The SURFRAD has been widely used to evaluate satellite SULR estimates [3,12,13,17]. The SULR data was aggregated to one-minute averages beginning in 2009. The data was interpolated to the exact overpass time of the satellite.

Temperature-Emissivity Method
LST and broadband emissivity products are required for the calculation of SULR using the temperature-emissivity method. The MODIS product MYD11 and MYD21 can both provide the LST and emissivity information, the equation of the temperature-emissivity method TE-MYD11 is shown in Equation (6): where the subscript MYD11 represents the TE-MYD11 method, εbb is the surface broadband emissivity that is calculated with the narrow-band emissivities, λ1 and λ2 are the spectral ranges (4-100 μm) for SULR, T is the surface temperature, and B stands for the Plank function. To calculate the SULR, three parameters are needed: εbb,T, and DLR.
For TE-MYD11 method, we use the level-2 1km swath product MYD11_L2 to provide the swath LST and the global level-3 daily 6km product MYD11B1 to provide the three narrow-band emissivities. Wan et al. [35][36][37] evaluated the accuracy of the MYD11_L2 LST product and the result shows that the LST accuracy is within 1K for relatively wide ranges of surface and atmospheric conditions. The emissivity in MYD11B1 is calculated using the day/night physical algorithm [38] that was developed by Li and Becker for the Advanced Very High Resolution Radiometer (AVHRR) data [39]. The accuracy of the narrow-band emissivity in MYD11B1 is reported to be 0.01 [35,36].
The equation of the temperature-emissivity method TE-MYD21 is shown in Equation (7), in which the subscript MYD21 represent the TE-MYD21 method, and other symbols remain the same as in Equation. (6). For the TE-MYD21 method, we use the level-2 1km swath product MYD21_L2 product to provide the LST and LSE. The water vapor scaling (WVS) method [40] is applied in the production of MYD21 to improve the atmospheric correction accuracy. The evaluation result of MYD21 product showed a similar accuracy (~1K) as MYD11 for LST over vegetated and ice/snow surfaces but a significant improvement in accuracy over dryland regions (from >3K of MYD11 to <1K of MYD21) and emissivity uncertainties below 0.015 using a LST and emissivity uncertainty simulator [41,42].
The SURFRAD in situ measured DLR values are used in Equations (6) and (7) to calculate the SULR. The broadband emissivity calculation method by Wang et al. [28] (Equation (8)) is adopted to calculate the broadband emissivity from three narrow-band emissivities of MODIS bands 29, 31, and 32 with a maximum error of 0.006.
where , , and are the MODIS narrow-band emissivity for MODIS bands 29, 31, and 32.

Hybrid Method
The hybrid method has long been used for the variable estimation in remote sensing with the assumption that the sensor recorded TOA radiance or brightness temperature contains the information to be retrieved. It mainly was used for the estimation of DLR [43] and was introduced into the estimation of SULR by Wang et al. [13] for its simplicity and good accuracy. The main steps of the hybrid method include the construction of the simulation dataset (see Section 2.1) and the regression based on the predefined model. The flowchart of the coefficient estimation for the hybrid methods is given in Figure 2. There are four hybrid models we discussed in this part, including three existing models (TOA-LIN, TOA-NLIN, and TOA-ANN) and one newly developed model (BOA-LIN). The main differences between the BOA-LIN and the TOA model are that the BOA model takes the surface outgoing radiances as SULR regression variables and a radiative transfer model (i.e., MODTRAN) is introduced in the calculation of the surface outgoing radiances from the satellite TOA radiances.

TOA-LIN
The TOA linear model used by Wang et al. [13] and Cheng et al. [3], a linear combination of TOA radiances for MODIS bands 29, 31, and 32 (Equation (9)), is directly adopted for consistency. It should be noted that the regression coefficients change with VZA, where is the VZA, , − , are the regression coefficients at , , , , , and , are MODIS TOA radiances for bands 29, 31, and 32, respectively. The coefficients were generated at different view zenith angles (VZA= 0° -60° with a 10° step) as shown in Table 3. The SULRVZA was calculated with a linear interpolation of two SULR values of adjacent VZAs.  [19] proposed a nonlinear TOA model method to estimate the SULR from the MODIS bands 31 and 32. The SULR is calculated with the introduction of an equivalent temperature ( ) that can be calculated from the TOA brightness temperature of the split-window bands and the VZA. The basic assumption of this model is that the TOA linear model does not eliminate the effect of atmospheric scattering and taking into account of the impact of atmospheric wv contents, thus may cause large errors in the estimation, and the split-window algorithm can minimize the effect of the wv. The model can be expressed as functions (10) and (11): where and are the TOA brightness temperatures of MODIS bands 31 and 32, respectively, M( ) is the radiation calculated by the Stefan-Boltzman law. , − , , , and are coefficients of the regression with the simulation dataset. The initial can be deduced by Planck law of the radiance in the wavelength corresponding to MODIS band 31 or 32, without loss of generality. Here, we use MODIS band 31.
The TOA-NLIN coefficients were generated at VZA from 0° to 60° with a 10° step as shown in Table 4. The ANN method has been widely used in the retrieval of environmental variables from remotely-sensed data. This is largely due to its ability to solve nonlinear problems such as estimation of near-surface air temperature [44], soil moisture [45], leaf area index (LAI) [46], and SULR [3,20]. With a suitable set of connecting weights and transfer functions, it has been shown that the multilayer perceptron can approximate any smooth, measurable function between the input and output vectors [47]. Back-propagation neural network (BPNN) is one of the most widely used ANN architectures. The basic BPNN structure is composed of an input layer, a hidden layer, and an output layer. Usually, more hidden layers and hidden neurons will lead to more accuracy fit results, but also may cause overfitting and an increase in computation time. After many iteration experiments of different hidden layers and hidden neurons, a network structure of 3-7-1 was chosen. The structure and training parameters of the BPNN model used for each VZA are as shown in Table 5. The surface outgoing radiances at the BOA are more closely related to the SULR than TOA radiances because the atmospheric attenuation and path radiation are potential error sources in the estimations process. Therefore, the BOA-LIN method is proposed here with the basic assumption which is similar to the TOA nonlinear model (the atmospheric correction is needed for SULR calculation with hybrid methods). The main difference between BOA model and TOA nonlinear model is that the BOA model chooses a physical atmospheric correction model (MODTRAN) to remove the effect of atmosphere using assimilated atmospheric profiles while the TOA-NLIN model uses a parametric approach (SW algorithm). The BOA-LIN method can be expressed as function (12): where , , , and are the regression coefficients, and , , , , and , are the surface outgoing radiances for bands 29, 31, and 32 after atmospheric correction. For a specific satellite TIR sensor, the surface outgoing radiance of band i can be calculated from TOA radiance as follows, where , is the surface outgoing radiance considering SRF of band , , is the TOA radiance of band , , and ↑, are the atmospheric transmittance and upwelling radiance considering the SRF of band , respectively.
The coefficients of BOA-LIN model are shown in Table 6. It should be noted that the VZA effect has been considered in the atmospheric correction using MODTRAN and thus the regression coefficients in the BOA linear model is VZA independent.

Results and Discussion
The estimation accuracy is described using three indicators: Root Mean Square Error (RMSE; Equation (14)), Mean Bias Error (MBE; Equation (15)), and R 2 .
where is the number of the samples, is the estimated SULR for the ith sample, and , is the simulated (in situ measured) truth value of SULR for the ith sample from the process of extensive radiative transfer simulation (ground measurement validation).

Results and Analysis Based on Simulated Datasets
The evaluation results for the four hybrid methods (TOA-LIN, TOA-NLIN, TOA-ANN, and BOA-LIN) based on the simulation datasets (see Section 2.1) are presented. The scatterplots between MODTRAN simulated and estimated SULR of four methods with VZA=0° are as shown in Figure 3. Details of the remaining VZA values for the TOA hybrid methods are given in Table 7.  Table 7. For the TOA hybrid methods, the fitting RMSE decreases with the increase of the model nonlinearity. The scatterplots of the TOA hybrid methods show that the bias increases with an increase of SULR values, while the BOA-LIN method has a relatively stable bias distribution at all SULR range.
The bias histograms and CDFs (Cumulative Distribution Functions) for the four hybrid methods are shown in Figure 4. It should be noted that the binsize for TOA hybrid methods is 1 W/m 2 and that for BOA-LIN is 0.5 W/m 2 . The bias histograms of all hybrid methods showed a normal distribution which centered at about 0 W/m 2 . Taking the bias values with CDF=1% and CDF=99% to indicate the bias level, we can see that the bias range for the BOA-LIN (−3.   Using accurate atmospheric parameters ( ↑, and ) the BOA linear hybrid method achieved a better theoretical accuracy than the other three hybrid methods in Figures 3,4. However, the atmospheric parameters are sensitive to the wv [48][49][50], which is one of the most temporally and spatially dynamic variables in the atmosphere. To better understand the influence of wv accuracy in the SULR estimation of the BOA-LIN method, additional simulation analysis with two levels of wv errors were taken. The wv in the 946 TIGR atmosphere profiles were scaled with the number of 0.95 and 1.05 (0.9 and 1.1) to calculate the biased atmospheric parameters for the evaluation of the performance of BOA-LIN with a 5% (10%) wv error. It should be noted that the estimation accuracy of the BOA-LIN method under certain wv error (e.g., 5% or 10%) is also VZA dependent as shown in Table 7. This is because with the increase of the VZA, the optical path increases and thus the accuracy of the MODTRAN calculated atmospheric parameters (τ, L↓ and L↑) decreases.
The results of BOA-LIN with accurate wv, 5% wv error (BOA-LIN_5%), and 10% wv error (BOA-LIN_10%), as well as the other three TOA hybrid methods from 0° to 60° are listed in Table 7. It can be seen that the TOA-LIN, BOA-LIN_5%, and BOA-LIN_10% have a larger RMSE and bias range with an increase of VZA. However, the results of VZA 10°-60° are only evident for TOA-NLIN and TOA-ANN because of their nonlinearity. The BOA-LIN with no added wv error has a very small RMSE (1.75 W/m 2 ) for all VZA values as expected. The RMSE of BOA-LIN_5% and BOA-LIN_10% increase sharply with the increase in VZA. All methods have a small MBE (<0.37 W/m 2 ) and a high R 2 (>0.993). Considering the mean RMSE values of all methods in Table 7, the performance in descending order is BOA-LIN, BOA-LIN_5%, TOA-ANN, TOA-NLIN, BOA-LIN_10%, and TOA-LIN.

Results and Analysis Based on in situ Measurements
The evaluation results for all the six methods (TE-MYD11, TE-MYD21, TOA-LIN, TOA-NLIN,  TOA-ANN, (Table 7) and proves the importance of taking the atmospheric correction into consideration. The TOA-NLIN method (with a RSME of 16.1 W/m 2 and MBE of −4.6 W/m 2 ) performs slightly better than the TOA-LIN method (with a RSME of 17.2 W/m 2 and MBE of −5.5 W/m 2 ). However, the TOA-ANN method (with a RSME of 21.2 W/m 2 and MBE of −2.5 W/m 2 ) performed the worst among all six methods which is contrary to the result shown in Table 6. Figure 5e shows that TOA-ANN has significant underestimation when SULR>600 W/m 2 , probably due to the overfitting of the ANN structure.

Six Methods Comparison for Daytime and Night-time Data
The validation results are separated for daytime and night-time as shown in Figure.6. The numbers of available in situ measurement clear-sky pixel for the daytime and night-time are 1422 and 1775, respectively. As we can see from the Figure 6a, the RMSE values in the daytime for all the six methods are greater than those in the night-time. The averaged RMSE (MBE) over the six methods in the daytime and night-time are 21.1 (−2.9) W/m 2 and 12.8 (−4.7) W/m 2 , respectively. The higher RMSE in the daytime can be explained by the larger SULR variation range in the daytime (with SULR standard derivation value of 103 W/m 2 in the daytime and 53 W/m 2 in the night-time) and significant thermal radiation directionality (TRD) in the daytime [22,51].
The RMSE values for the daytime range from 18.0 W/m 2 to 27.6 W/m 2 and from 9.5 W/m 2 to 14.9 W/m 2 during the night-time. As seen in Figure.

Six Methods Comparison for Each Site
The RMSE and MBE values for all six methods over seven SURFRAD sites are listed in Table 8. The accuracy varies for different observation sites. In general, the results over FPE and PSU are the best. The RMSE of SXF, BOU, GCR, and BON are in the middle, and the result of DRA is the worst. The DRA site has an obvious negative bias (from −8.8 W/m 2 to −25.3 W/m 2 ) and a relatively high RMSE (from 15.9 W/m 2 to 33.1 W/m 2 ) which is similar to the previous studies of Wang et al. [13] and Cheng et al. [3]. The performance at the DRA site for TE-MYD21 (RMSE=15.9 W/m 2 and MBE=-8.8 W/m 2 ) and BOA-LIN (RMSE = 17.5 W/m 2 and MBE = −11.3 W/m 2 ) is much better than the four remaining methods. This may be due to more accurate LST values from the MYD21 product over barren surfaces compared with the MYD11 product [52,53] and the accurate atmospheric correction in the BOA-LIN method, respectively. The TE-MYD11 and TE-MYD21 method use the same algorithm and same source of in situ measured DLR but different source of LST and LSE as described in section 3.1. The significant improvement of TE-MYD21 method over DRA barren site has been discussed in Section 4.2.3. The performances of TE-MYD21 and TE-MYD11 over barren and non-barren surfaces are further analyzed at each site, considering that the non-barren sites (sites except DRA) will perform like barren in some seasons because of the vegetation annual growth cycle or man's activities (e.g. tillage practices on cropland). The NDVI data at the seven SURFRAD sites is retrieved from the MODIS products (MOD13A2 and MYD13A2, 16-day composed with 1 km resolution) and used to describe the richness of vegetation. The retrieved NDVI values are shown in Figure 7;  The land surfaces with NDVI lower than 0.3 is recognized as barren surfaces in this study. The percentage and mean NDVI of barren and non-barren observations at the seven SURFRAD sites and corresponding MBEs and RMSEs for the TE-MYD21 and TE-MYD11 methods are summarized in Table 9. The available observations are all barren (NDVI<0.3) for DRA site and are all non-barren (NDVI>=0.3) for GCR and PSU sites. For barren surfaces, the TE-MYE21 performs better than the TE-MYD11 method at four sites (DRA, FPE, SXF, and BOU) among five available sites and the weighted mean RMSE (MBE) are 14.2 W/m 2 (-4.2 W/m 2 ) and 22.4 W/m 2 (-14.9 W/m 2 ) for TE-MYD21 and TE-MYD11 methods. For non-barren surfaces, the TE-MYE11 performs slightly better than the TE-MYD21 method at four sites (BON, PSU, SXF, and BOU) among six available sites and the weighted mean RMSE (MBE) are 13.8 W/m 2 (3.7 W/m 2 ) and 13.6 W/m 2 (-2.5 W/m 2 ) for TE-MYD21 and TE-MYD11 methods. The TE-MYD21 shows a much better accuracy than TE-MYD11 method at barren surfaces and a similar accuracy at non-barren surfaces. The TOA hybrid methods only require TOA radiances and coefficients when calculating the SULR, while the BOA-LIN method needs additional inputs of atmospheric profiles for atmospheric correction. The performance of the BOA-LIN method highly depends on the wv accuracy of the atmospheric profile as shown in Table 6. In this section, the RMSE values of four hybrid methods with different wv values (step=0. 5) are compared.
The wv histogram (Figure 8a) shows that 93.0% of the observations have a wv in the range of [0,4]g/cm 2 . The observations with wv>4.0 g/m 2 are not considered here because of the low frequency. The RMSE comparison of four hybrid methods at different wv are shown in Figure 8b. It can be seen that the BOA-LIN method has the lowest RMSE for wv<3 g/m 2 and a relatively small RMSE for 3 g/m 2 <wv<4 g/m 2. The TOA-ANN has the largest RMSE for wv>1 g/cm 2 . The TOA-LIN and TOA-NLIN show an intermediate accuracy and the later one is more accurate for wv<2 g/cm 2 . The MBE comparison of four hybrid methods at different wv are shown in Figure 8c. The BOA-LIN method has the most stable MBE likely due to the step of atmospheric correction. The TOA-ANN method has the most variable MBE maybe due to the over fitting of the ANN-structure. TOA-NLIN method has a MBE closer to zero than the TOA-LIN for wv<2 g/cm 2 similar to the RMSE. The TE-MYD21 and BOA-LIN are shown to be the best methods in their respective categories while the TE-MYD21 performed slightly better than the BOA-LIN method (see Table 8). The DLR used in TE-MYD21 was in situ measured data however in the global SULR generation, the DLR is usually obtained from remote sensing or reanalysis products. It is reported that the DLR estimation can achieve an accuracy of ~25 W/m 2 [4,48]. The RMSE of TE-MYD21 method with varying DLR errors were compared with RMSE of BOA-LIN in Figure 9. A 50 W/m 2 bias of DLR leads to a <0.074 W/m 2 bias to the RMSE of the TE-MYD21. Therefore, the TE-MYD21 method likely performs better than the BOA-LIN in applications. Nevertheless, the BOA-LIN method is also a good choice for the satellite sensors that cannot utilize the TES algorithm to estimate accurate LST and LSE.

Conclusion
This paper evaluated six clear-sky SULR estimation methods with simulation datasets and Aqua-MODIS measurements. The recently released MYD21 product was used to estimate the SULR. A new method (BOA-LIN) was developed under the general framework of the hybrid method but also includes atmospheric correction to eliminate the effect of atmosphere. Some of the new insights into the SULR estimation methods are as follows.
(1) For the theoretical analysis of TOA hybrid methods (TOA-LIN, TOA-NLIN, TOA-ANN), the fitting RMSE decreases with increasing model nonlinearity. The fitting RMSE of BOA-LIN (1.75 W/m 2 ) is much less than the RMSE of the TOA hybrid methods assuming an accurate atmospheric correction has been achieved. The performance of BOA-LIN decays with the increase of atmospheric profile wv error. The BOA hybrid method has great potential in application if accurate atmospheric profiles can be obtained as input.
(2) The TE-MYD21 performs the best among all the six methods with RMSE of 14.0 W/m 2 and MBE of −0.2 W/m 2 , and the BOA-LIN performs best among the four hybrid methods with RMSE of 15.2 W/m 2 and MBE of −2.3 W/m 2 based on the two-year satellite products. The performance of six methods in descending order is TE-MYD21, BOA-LIN, TOA-NLIN, TOA-LIN, TE-MYD11, and TOA-ANN. TE-MYD21 has a much better accuracy than the TE-MYD11 over barren surfaces (NDVI<0.3) and a similar accuracy over non-barren surfaces (NDVI>=0.3). The BOA-LIN is more accurate than other TOA hybrid methods due to the inclusion of atmospheric correction.
SULR is an important parameter in the estimation of the radiation budget that needs to be estimated accurately with remotely sensed data. More adequate evaluation is needed in the future with more spatial locations and longer temporal coverage. All six SULR estimation methods presented here are under the thermal isotropic assumption because MODIS can only supply one viewing angle observation. Further studies should focus on the development of SULR calculation methods that take into consideration the thermal directionality.