An Improved Parameterization for Retrieving Clear-Sky Downward Longwave Radiation from Satellite Thermal Infrared Data

Surface downward longwave radiation (DLR) is a crucial component in Earth’s surface energy balance. Yu et al. (2013) developed a parameterization for retrieving clear-sky DLR at high spatial resolution by combined use of satellite thermal infrared (TIR) data and column integrated water vapor (IWV). We extended the Yu2013 parameterization to Moderate Resolution Imaging Spectroradiometer (MODIS) data based on atmospheric radiative simulation, and we modified the parameterization to decrease the systematic negative biases at large IWVs. The new parameterization improved DLR accuracy by 1.9 to 3.1 W/m2 for IWV≥3 cm compared to the Yu2013 algorithm. We also compared the new parameterization with four algorithms, including two based on Top-of-Atmosphere (TOA) radiance and two using near-surface meteorological parameters and water vapor. The algorithms were first evaluated using simulated data and then applied to MODIS data and validated using surface measurements at 14 stations around the globe. The results suggest that the new parameterization outperforms the TOA-radiance based algorithms in the regions where ground temperature is substantially different (enough that the difference between them is as large as 20 K) from skin air temperature. The parameterization also works well at high elevations where atmospheric parameter-based algorithms often have large biases. Furthermore, comparing different sources of atmospheric input data, we found that using the parameters interpolated from atmospheric reanalysis data improved the DLR estimation by 7.8 W/m2 for the new parameterization and 19.1 W/m2 for other algorithms at high-altitude sites, as compared to MODIS atmospheric products.


Introduction
The surface downward longwave radiation (DLR) is radiative flux density reaching the surface in the thermal infrared (TIR) part of the electromagnetic spectrum (4-100 µm) and represents a key component in land-atmosphere interactions.DLR is related to the downward radiance at the Earth's surface as [1]: where I v (z = 0, −µ) is the spectral radiance at the surface at a given local zenith angle θ and µ = cos(θ) and ν is the wavenumber.Under cloud-free conditions, multiple scatter can be omitted, thus spectral radiance can be expressed as: where z t is the top of atmosphere, B v (z ) is the Planck function evaluated with the temperature at level z , and T v (0, z ; −µ) is the monochromatic transmittance from surface to height z .A global, long-term DLR database is needed for studies of surface energy budget and global climate change [2][3][4][5].Although surface radiation fluxes can be obtained by field measurement, their sparsity makes downward radiation estimation difficult over larger regions or regions without in-situ measurement stations, while remote sensing provides more complete spatial coverage.
Surface DLR results from atmospheric absorption, emission, and scattering within the entire atmospheric column.If detailed atmospheric information is available, including vertical profiles of temperature and humidity, ozone and other trace gases, aerosols, as well as cloud physical properties such as cloud height and cloud temperature, DLR can be accurately calculated by detailed radiative-transfer (RT) modeling.This method has been applied to data from the TIROS Operational Vertical Sounder (TOVS) of National Oceanic and Atmospheric Administration (NOAA) satellites [6], International Satellite Cloud Climatology Project (ISCCP, [7]), and Clouds and the Earth's Radiant Energy System (CERES, [8]).
Considering that RT methods require extensive computation and multiple input parameters that are unavailable for many satellites, researchers developed more practical algorithms, including empirical and parameterized algorithms.Because the major part of surface DLR originates from the first few hundred meters of the atmospheric layer above the earth surface under clear sky conditions [9], the empirical algorithms establish empirical relationships between clear-sky DLR and meteorological parameters such as near-surface air temperature and humidity [10][11][12][13].Such algorithms have been applied to the Meteosat series [14][15][16], Moderate Resolution Imaging Spectroradiometer (MODIS, [17,18]), Geostationary Operational Environmental Satellite (GOES, [19]), and CERES data [20].Empirical algorithms are usually developed using in-situ observations from certain surface sites and must be adjusted for other locations, which makes them unsuitable for use at larger scales.The parameterized algorithms express DLR as a function of key atmospheric parameters, and most of these algorithms determine the relationship between DLR and input parameters through extensive radiative-transfer model calculation and statistical regression.Gupta et al. [21,22] defined clear-sky DLR in terms of atmospheric effective temperature and column-integrated water vapor (IWV).Inamdar and Ramanathan [23] established a relationship between clear-sky DLR and satellite-retrieved parameters, including top-of-atmosphere (TOA) outgoing longwave radiation, IWV, land surface temperature, and near-surface air temperature.Zhou et al. [24,25] determined clear-sky DLR from surface upward longwave radiation and IWV.
Meteorological data used by the last two types of algorithms come from either atmospheric reanalysis datasets or satellite retrieval products, thus algorithm error is affected by the retrieval accuracy of atmospheric parameters [1].In addition, atmospheric data are often at coarse spatial resolution (≥1 • ), which may have large errors when they are directly used at fine scale.Current global radiation products such as ISCCP-FD, Global Energy and Water Cycle Experiment-Surface Radiation Budget (GEWEX-SRB), and CERES-FSW were found to have large errors in high-altitude regions such as the Qinghai-Tibetan Plateau, which were mainly caused by errors in surface air temperature and water vapor content [26][27][28].Zhang et al. [29] pointed out that an uncertainty of ~2-4 K in surface air temperature would induce ~15 W/m 2 of uncertainty in the DLR of ISCCP.
To avoid propagating retrieval errors of satellite-derived atmospheric parameters into final radiation estimates, some researchers have developed hybrid algorithms that use satellite-observed thermal infrared (TIR) data to calculate clear-sky DLR directly.Algorithms for the high-resolution infrared spectrometer (HIRS/2) of NOAA [1], MODIS [30][31][32][33][34], and GOES Sounder [35] have been developed during the past two decades.In these algorithms, the DLR is expressed by linear or nonlinear combinations of multispectral infrared radiances or brightness temperature (BT) at the TOA.Because this type of algorithm is not only free from real-time atmospheric and surface information but also can derive DLR at the same spatial resolution as satellite thermal data, it is regarded as a promising trend in radiative flux estimation [32].However, hybrid models that use TOA radiance or land surface temperature always greatly overestimate DLR during daytime over dry-arid regions where surface skin temperature is significantly higher than near-surface air temperature [27,30].
To overcome the DLR overestimation over dry-arid regions for TOA radiance-based algorithms, Yu et al. [36] (hereafter Yu2013) proposed a correction based on ground minus air temperature difference, denoted as δT s,a .They developed a parameterization that estimates clear-sky DLR from the BT of Huan Jing InfraRed Scanner (HJ-1B IRS) thermal channel (10.8 µm) and the synchronous MODIS water vapor product.They found that taking into account δT s,a improved DLR at dry-arid sites.The Yu2013 parameterization provides a new approach to obtaining DLR from satellite TOA data.This algorithm is not restricted to satellites with multispectral infrared channels, such as MODIS, but can also be used for satellites with a single thermal channel, such as Fengyun-2E (FY-2E) and HJ-1B IRS.However, Yu et al. [36] also found obvious negative errors under high IWV values upon validation of the parameterization by a testing dataset that was simulated using a radiative transfer model and global radio-soundings.
The purpose of this paper is to extend Yu2013 parameterization to MODIS data by coefficient reconstruction, taking advantage of the high-quality input data from MODIS to improve the parameterization and decrease the systematic negative biases at high IWV values and to compare the improved parameterization with several other state-of-the-art algorithms with good performance.The rest of the paper is structured as follows: the algorithm introduction, extension, and improvement are described in Section 2. Section 3 assess the improved algorithm based on a simulated dataset and actual MODIS Terra data, respectively.Section 4 further discusses the effect of different sources of atmospheric products on DLR estimation.Section 5 gives the conclusion.

Clear-Sky DLR Algorithms
As mentioned in the above section, two types of algorithms have been widely used in DLR estimation, i.e., parameterizations using atmospheric input data and hybrid algorithms using TOA radiance of TIR channels.Table 1 gives four widely used or newly developed algorithms.The algorithms proposed by Zhou et al. [25] and Gupta et al. [22], hereafter referred to as Zhou-Cess and Gupta2010, were used to produce global radiation flux by the CERES Surface Radiation Budget (SRB) team [37].The two TOA radiance algorithms developed by Tang and Li [30] and Liang and Wang [31], hereafter referred to as Tang-Li and Wang-Liang, define clear-sky DLR by MODIS TOA radiance or brightness temperature.Zhou-Cess, Gupta2010, and Tang-Li algorithms were developed using global distributed atmospheric profiles, whereas the Wang-Liang algorithm was developed based on profiles over the North American continent only.
The Yu2013 parameterization is a hybrid of the two types of algorithms above.This parameterization was developed at the hypothesis that the atmosphere is a grey body, and clear-sky DLR can be written as σε a T 4 e , where σ is Stefan-Boltzmann constant, ε a is effective emissivity of the entire atmospheric column, and T e is the effective temperature of the atmosphere.Considering that the TOA BT of TIR channels sensitive to near surface (without making atmospheric corrections) has strong correlation with near-surface air temperature T a [1,15,36], which is relevant to the atmospheric effective temperature under a clear sky, BT is used to represent the effective temperature [1,36].Then, the definition of atmospheric effective emissivity should be adjusted when TOA BT substitutes for atmospheric temperature: where T k is BT of the TIR channel, and subscript k denotes the channel number.Because atmospheric effective emissivity ε k can be empirically specified as a function of atmospheric IWV by regression analysis [1,36]: where IWV has units of cm, and a i (i = 0, 2) represent coefficients.Therefore, DLR is expressed as: For the Yu2013 algorithm, T k used the TIR channel of HJ1B/IRS with a center wavelength of 10.8 µm.The regression coefficients are the look up tables (LUTs) of elevation H, satellite view zenith angle VZA, and δT s,a , where δT s,a = T s -T a .For ordinary regions that T s is similar to T a , the regression coefficients are be written as a i = a i (H, VZA).For the condition that T s is very distinct from T a , such as in arid regions, a i = a i (H, VZA, δT s,a ).The coefficients were regressed using the dataset simulated by the MODerate resolution atmospheric TRANsmission (MODTRAN) radiation model [38] and thermodynamic Initial Guess Retrieval (TIGR) atmospheric profiles [39].

Authors
Abbreviation Algorithm 1   Zhou et al. [25] Zhou-Cess Tang and Li [30] Tang-Li DLR = a 0 + 1 L i is the top-of-atmosphere (TOA) irradiance of channel i; a i and b i are regression coefficients, SULW is the surface upwelling longwave flux computed using 2-m air temperature T a with an emissivity of unity, integrated water vapor (IWV) is in unit of cm, T e is effective temperature which computed from the surface skin temperature T s and the mean temperatures of the lowest two atmospheric layers (T 1 and T 2 ), and H is surface elevation.

Extension of Yu2013 Parameterization to MODIS
When extending the Yu2013 algorithm to MODIS satellite data, T k in ( 5) is replaced by BT of MODIS channel 31 (T 31 ) with a central wavelength of 10.8 µm, thus ε k is expressed as ε 31 .As shown in Figure 1, BT has a strong correlation with atmospheric effective temperature, thus the algorithm is feasible for MODIS.The flowchart of algorithm construction is shown in Figure 2. Firstly, DLR and TOA radiance of MODIS channel 31 were simulated by atmospheric radiative transfer model and global profiles.Then, the relationship between DLR and MODIS T 31 was derived by regression analysis-actually, the relationship between ε 31 and IWV was derived.Finally, the algorithm was validated using simulated data and actual data.

Training and Testing Datasets Construction
In order to adequately cover atmospheric variations, two global atmospheric profile databases were used to simulate the training dataset.The atmospheric profile database TIGR2002 of the previous study was used [36].TIGR2002 includes 2311 atmospheric profiles of temperature, relative humidity, and ozone at 40 fixed pressure levels from 1013 to 0.05 hPa.However, TIGR2002 lacks the profiles with high IWV [36].Therefore, a database of 13,495 profiles generated by Chevallier et al. [40] was also used.These profiles were sampled from European Centre for Medium-Range Weather Forecast (ECMWF) 40-year re-analyses, which are globally distributed and representative of a wide range of atmospheric conditions.The ECMWF profiles include temperature, relative humidity, and ozone in 60 fixed pressure levels from 1013 to 0.1 hPa.We converted the ECMWF profiles to the same pressure levels as TIGR2002 using the dataset self-contained program.Because only clear-sky conditions were considered, profiles with relative humidity <85% in each layer were selected.Finally, 875 TIGR2002 and 7730 ECMWF profiles were used in the simulation.
In order to adequately cover atmospheric variations, two global atmospheric profile databases were used to simulate the training dataset.The atmospheric profile database TIGR2002 of the previous study was used [36].TIGR2002 includes 2311 atmospheric profiles of temperature, relative humidity, and ozone at 40 fixed pressure levels from 1013 to 0.05 hPa.However, TIGR2002 lacks the profiles with high IWV [36].Therefore, a database of 13,495 profiles generated by Chevallier et al. [40] was also used.These profiles were sampled from European Centre for Medium-Range Weather Forecast (ECMWF) 40-year re-analyses, which are globally distributed and representative of a wide range of atmospheric conditions.The ECMWF profiles include temperature, relative humidity, and ozone in 60 fixed pressure levels from 1013 to 0.1 hPa.We converted the ECMWF profiles to the same pressure levels as TIGR2002 using the dataset self-contained program.Because only clear-sky conditions were considered, profiles with relative humidity < 85% in each layer were selected.Finally, 875 TIGR2002 and 7730 ECMWF profiles were used in the simulation.For each atmospheric profile, the DLR, TOA radiance, and BT of MODIS thermal channels were simulated for various surface types, viewing angles, elevations, and δTs,a.Version 4.0 of the MODTRAN radiation model was used for simulation.Downward surface radiances at nineteen different zenith angles (0 o to 85 o at 5 o interval, and 89 o ) were computed for wave numbers ranging between 100 and 3333 cm −1 at a resolution of 1 cm −1 .Such values were then integrated to obtain the total DLR according to Equation (1).The model was configured to use ozone and trace gases from a

Training and Testing Datasets Construction
In order to adequately cover atmospheric variations, two global atmospheric profile databases were used to simulate the training dataset.The atmospheric profile database TIGR2002 of the previous study was used [36].TIGR2002 includes 2311 atmospheric profiles of temperature, relative humidity, and ozone at 40 fixed pressure levels from 1013 to 0.05 hPa.However, TIGR2002 lacks the profiles with high IWV [36].Therefore, a database of 13,495 profiles generated by Chevallier et al. [40] was also used.These profiles were sampled from European Centre for Medium-Range Weather Forecast (ECMWF) 40-year re-analyses, which are globally distributed and representative of a wide range of atmospheric conditions.The ECMWF profiles include temperature, relative humidity, and ozone in 60 fixed pressure levels from 1013 to 0.1 hPa.We converted the ECMWF profiles to the same pressure levels as TIGR2002 using the dataset self-contained program.Because only clear-sky conditions were considered, profiles with relative humidity < 85% in each layer were selected.Finally, 875 TIGR2002 and 7730 ECMWF profiles were used in the simulation.For each atmospheric profile, the DLR, TOA radiance, and BT of MODIS thermal channels were simulated for various surface types, viewing angles, elevations, and δTs,a.Version 4.0 of the MODTRAN radiation model was used for simulation.Downward surface radiances at nineteen different zenith angles (0 o to 85 o at 5 o interval, and 89 o ) were computed for wave numbers ranging between 100 and 3333 cm −1 at a resolution of 1 cm −1 .Such values were then integrated to obtain the total DLR according to Equation (1).The model was configured to use ozone and trace gases from a For each atmospheric profile, the DLR, TOA radiance, and BT of MODIS thermal channels were simulated for various surface types, viewing angles, elevations, and δT s,a .Version 4.0 of the MODTRAN radiation model was used for simulation.Downward surface radiances at nineteen different zenith angles (0 • to 85 • at 5 • interval, and 89 • ) were computed for wave numbers ranging between 100 and 3333 cm −1 at a resolution of 1 cm −1 .Such values were then integrated to obtain the total DLR according to Equation (1).The model was configured to use ozone and trace gases from a climatological data set and the "Rural" aerosol with horizontal meteorological visible thickness 23 km, all available from MODTRAN.Seven surface types (conifer, dry grass, sea water, fine snow, wetland, concrete, and soil) with distinctly different emissivity spectra, seven elevations, eight VZA, and nine δT s,a values were used to simulate TOA radiance and brightness temperature.Surface emissivity spectra were from the ASTER spectral library, and their spectral emissivity are described in Figure 3.To obtain atmospheric profiles at various elevations, profiles from 40 to 34 layers were used, and average elevations of the lowest seven layers were 0.002, 0.454, 0.934, 1.410, 1.885, 2.670, and 3.511 km.
The VZA varied from 0 • to 70 • at a 10 • interval.T s was assumed to be different from T a , and the δT s,a was between −20 K and 20 K at 5 K intervals.Knowing DLR and BT, ε 31 defined by BT were calculated using Equation (3).IWV for each condition was calculated by MODTRAN directly.
km, all available from MODTRAN.Seven surface types (conifer, dry grass, sea water, fine snow, wetland, concrete, and soil) with distinctly different emissivity spectra, seven elevations, eight VZA, and nine δTs,a values were used to simulate TOA radiance and brightness temperature.Surface emissivity spectra were from the ASTER spectral library, and their spectral emissivity are described in Figure 3.To obtain atmospheric profiles at various elevations, profiles from 40 to 34 layers were used, and average elevations of the lowest seven layers were 0.002, 0.454, 0.934, 1.410, 1.885, 2.670, and 3.511 km.The VZA varied from 0° to 70° at a 10° interval.Ts was assumed to be different from Ta, and the δTs,a was between −20 K and 20 K at 5 K intervals.Knowing DLR and BT, ε31 defined by BT were calculated using Equation (3).IWV for each condition was calculated by MODTRAN directly.The testing dataset was used to evaluate whether the new algorithm was appropriate under various atmospheric conditions and to quantify algorithm accuracy under error-free conditions of the input parameters.A global, cloud-free atmospheric profile database, Cloudless Land Atmosphere Radiosounding (CLAR, [41]), was used.CLAR contains 382 radiosoundings taken at uniformly distributed meteorological stations and has a wide range of temperature, humidity, and elevation.Thus, the database can serve as an ideal testing dataset for accuracy evaluation.The testing dataset was generated using the same method as the training dataset except that only these atmospheric profiles at the surface elevation were used, and the ozone was from the climatological data set provided by MODTRAN.

Improvement of DLR Underestimation at High Water Vapor Content
Yu et al. [36] found that DLRs predicted by Equation ( 5) had increasing negative errors with IWV when IWV was large.The underestimation was caused by a significant underestimation of atmospheric emissivity at large IWV, as shown by the fit1 line of Figure 4a.Adding a large number of profiles with large IWV to derive the coefficients can obviously decrease the negative errors, but the underestimation persists (Figure 4b).The testing dataset was used to evaluate whether the new algorithm was appropriate under various atmospheric conditions and to quantify algorithm accuracy under error-free conditions of the input parameters.A global, cloud-free atmospheric profile database, Cloudless Land Atmosphere Radiosounding (CLAR, [41]), was used.CLAR contains 382 radiosoundings taken at uniformly distributed meteorological stations and has a wide range of temperature, humidity, and elevation.Thus, the database can serve as an ideal testing dataset for accuracy evaluation.The testing dataset was generated using the same method as the training dataset except that only these atmospheric profiles at the surface elevation were used, and the ozone was from the climatological data set provided by MODTRAN.

Improvement of DLR Underestimation at High Water Vapor Content
Yu et al. [36] found that DLRs predicted by Equation ( 5) had increasing negative errors with IWV when IWV was large.The underestimation was caused by a significant underestimation of atmospheric emissivity at large IWV, as shown by the fit1 line of Figure 4a.Adding a large number of profiles with large IWV to derive the coefficients can obviously decrease the negative errors, but the underestimation persists (Figure 4b).
Considering that the DLR underestimation was caused by underestimation of atmospheric emissivity ε 31 , we intended to develop a robust formulation of atmospheric emissivity ε 31 .Water vapor is the most important atmospheric gas contributing to clear-sky DLR, and the variation of other gases (such as CO 2 and O 3 ) have small effects on DLR [32,34].Therefore, we used IWV to predict atmospheric emissivity ε 31 .As shown in Figure 4, the atmospheric emissivity ε 31 increased with atmospheric IWV and then was nearly unchanged when IWV achieved a specific value.By data analysis, we found that using a logarithmic form of IWV predicted the relationship between atmospheric emissivity ε 31 and IWV well, which is consistent with the results in previous studies [22,25,42].Since the logarithmic function decreases very rapidly with decreasing water vapor below 1 cm, replacing ln(IWV) by ln(1+IWV) can avoid flux underestimation under low water vapor conditions [25].Finally, based on the regression result, the emissivity term of Equation ( 5) is modified to the following form: where 1 + IWV is used to ensure that V is greater than zero under low water vapor conditions.In Equation ( 6), V is expressed by ln(1 + IWV) rather than ln(1 + IWV) because DLR is overestimated at very large IWV when using ln(1 + IWV) (as shown in Figure 4a).Figure 4 compares atmospheric emissivity ε 31 fitted from three different formulations.Using √ IWV or ln(1 + IWV) will underestimate or overestimate the atmospheric emissivity ε 31 at very large IWV, especially when the number of training profiles is small (Figure 4a).On the contrary, using ln(1 + IWV) produces much better results than the two former formulations, and the result is not restricted by sample number.Considering that the DLR underestimation was caused by underestimation of atmospheric emissivity ε31, we intended to develop a robust formulation of atmospheric emissivity ε31.Water vapor is the most important atmospheric gas contributing to clear-sky DLR, and the variation of other gases (such as CO2and O3) have small effects on DLR [32,34].Therefore, we used IWV to predict atmospheric emissivity ε31.As shown in Figure 4, the atmospheric emissivity ε31 increased with atmospheric IWV and then was nearly unchanged when IWV achieved a specific value.By data analysis, we found that using a logarithmic form of IWV predicted the relationship between atmospheric emissivity ε31 and IWV well, which is consistent with the results in previous studies [22,25,42].Since the logarithmic function decreases very rapidly with decreasing water vapor below 1 cm, replacing ln(IWV) by ln(1+IWV) can avoid flux underestimation under low water vapor conditions [25].Finally, based on the regression result, the emissivity term of Equation ( 5) is modified to the following form: ( ) where 1+IWV is used to ensure that V is greater than zero under low water vapor conditions.In Equation ( 6), V is expressed by ln(1 ) IWV + rather than ln(1+IWV) because DLR is overestimated at very large IWV when using ln(1+IWV) (as shown in Figure 4a).Figure 4 compares atmospheric emissivity ε31 fitted from three different formulations.Using IWV or ln(1+IWV) will underestimate or overestimate the atmospheric emissivity ε31 at very large IWV, especially when the number of training profiles is small (Figure 4a).On the contrary, using ln(1 ) IWV + produces much better results than the two former formulations, and the result is not restricted by sample number.
The new formulation can greatly improve results compared with the Yu2013 algorithm at large IWV, especially when the profile number for regression is not large.Figure 5 shows the results of the The new formulation can greatly improve results compared with the Yu2013 algorithm at large IWV, especially when the profile number for regression is not large.Figure 5 shows the results of the two parameterizations when their coefficients were derived using only TIGR2002 profiles.The atmospheric emissivity was underestimated at large IWV for the Yu2013 parameterization, causing negative bias of DLR (Figure 5a,c).With the new parameterization, the negative bias of atmospheric emissivity was corrected, thereby improving the predicted DLRs (Figure 5b,d).As shown in the figure, the new parameterization reduced the root mean square error (RMSE) by 1.9 to 3.1 W/m 2 for δT s,a = 0 K.The DLR statistical error distribution for IWV greater than 3.0 cm is shown in Table 2.For small VZA (<60 • ), more than 99% of the DLRs had accuracy within 20 W/m 2 , an improvement of 0.5% compared with Yu2013 algorithm.At large VZA (=60 • ), the accuracy increased by 16.9%.
Subsequently, we derived a LUT of coefficients for the new parameterization at various elevation, VZA, and δT s,a values from the training dataset.We found the coefficients did not yield obvious differences between surface types except for dry grass, the same result as that in [36], because the emissivity spectrum of dry grass is obviously lower than other surface types in the thermal infrared range.Therefore, two sets of coefficients were derived, dry grass-only coefficients for arid surfaces and the all surfaces-average coefficients for other surfaces.Moreover, because δT s,a is difficult to retrieve accurately using MODIS data, we suggest that δT s,a is only used to surface with large δT s,a , such as bare ground, arid region, snow, ice, etc.Also, we suggest setting δT s,a within ±20 K in the DLR calculation to limit the error caused by inaccurate δT s,a .
Table 2.For small VZA (< 60°), more than 99% of the DLRs had accuracy within 20 W/m 2 , an improvement of 0.5% compared with Yu2013 algorithm.At large VZA (= 60°), the accuracy increased by 16.9%.Subsequently, we derived a LUT of coefficients for the new parameterization at various elevation, VZA, and δTs,a values from the training dataset.We found the coefficients did not yield obvious differences between surface types except for dry grass, the same result as that in [36], because the emissivity spectrum of dry grass is obviously lower than other surface types in the thermal infrared range.Therefore, two sets of coefficients were derived, dry grass-only coefficients for arid surfaces and the all surfaces-average coefficients for other surfaces.Moreover, because δTs,a

Evaluation of Algorithm Accuracy
DLRs were calculated from IWV and the simulated TOA BT of the testing dataset using the new parameterization and evaluated using MODTRAN-simulated DLRs.Table 3 indicates that most surfaces had similar retrieval accuracy (RMSE = 5.1-8.4W/m 2 ), except of dry grass, which has a large negative error (bias= −6.0 W/m 2 ).We further gave the DLR error statistics under different climatic conditions for the conifer surface.We found the RMSEs of daytime and nighttime, with and without temperature inversion, were 6.2, 4.8, 5.5, and 5.6 W/m 2 , respectively.Daytime DLRs had a larger positive error (~3.4 W/m 2 ) than nighttime, and DLRs under temperature inversion conditions had a bias about −3.1 W/m 2 greater than the conditions without temperature inversion.Based on detailed Remote Sens. 2019, 11, 425 9 of 27 analysis, we found that the temperature inversion had a very small effect on DLR estimation for this algorithm, except that the near-surface temperature profile had a very thick layer with a strong lapse rate.Furthermore, some polar winter profiles had large negative errors.Figure 6a,d gives RMSE and bias variation with VZA for the conifer surface.RMSE increased with VZA from 5.6 to 20.2 W/m 2 for δT s,a = 0 K, and the bias was 1.8 to −6.8 W/m 2 .RMSE were 5.4 to 20.0 W/m 2 and 6.7 to 20.8 W/m 2 for δT s,a = −20 and 20 K, respectively.Table 3.The mean value of actual and calculated DLRs (DLR act and DLR cal ), root mean square error (RMSE), bias, maximum error and minimum error between actual DLRs, and those calculated from training datasets for different surface types when δT s,a = 0 K and VZA = 0 • (unit: W/m 2 ).surfaces had similar retrieval accuracy (RMSE=5.1-8.4W/m 2 ), except of dry grass, which has a large negative error (bias= −6.0 W/m 2 ).We further gave the DLR error statistics under different climatic conditions for the conifer surface.We found the RMSEs of daytime and nighttime, with and without temperature inversion, were 6.2, 4.8, 5.5, and 5.6 W/m 2 , respectively.Daytime DLRs had a larger positive error (~ 3.4 W/m 2 ) than nighttime, and DLRs under temperature inversion conditions had a bias about −3.1 W/m 2 greater than the conditions without temperature inversion.Based on detailed analysis, we found that the temperature inversion had a very small effect on DLR estimation for this algorithm, except that the near-surface temperature profile had a very thick layer with a strong lapse rate.Furthermore, some polar winter profiles had large negative errors.Figure 6a and 6d gives RMSE and bias variation with VZA for the conifer surface.RMSE increased with VZA from 5.6 to 20.2 W/m 2 for δTs,a = 0 K , and the bias was 1.8 to −6.8 W/m 2 .RMSE were 5.4 to 20.0 W/m 2 and 6.7 to 20.8 W/m 2 for δTs,a = −20 and 20 K, respectively.The new algorithm was compared with the four algorithms in Section 2.1.Though these algorithms were developed under different conditions, it was worth investigating their advantages and drawbacks by applying them to the testing data.The results of algorithms based on testing data are shown in Figures 6 and 7.As shown in Figure 6, the TOA radiance-based algorithms, Tang-Li and Wang-Liang, performed poorer than the new algorithm, and their DLR accuracies were largely dependent on δT s,a and VZA.The Tang-Li algorithm had the best result at δT s,a = 0 K and performed worse when the absolute value of δT s,a became larger.This is because Ta was considered equal to Ts in the algorithm development [30].For the Wang-Liang algorithm, RMSE and bias were minimal at δT s,a = 10 K and increased at other δT s,a values.Though MODIS-derived Ta and Ts were used in the MODTRAN simulation during their algorithm development [31], the coefficients were not regressed for each δT s,a value separately.Therefore, they represented the average δT s,a of all training profiles.When VZA ≤ 40 • , DLR RMSE is 17.0 to 23.8 W/m 2 for the Tang-Li algorithm at δT s,a = 0 K and 18.7 to 31.4 W/m 2 for the Wang-Liang algorithm at δT s,a = 10 K.When VZA was larger, RMSEs of the two algorithms increased rapidly with VZA, especially for Wang-Liang algorithm, with RMSE reaching 151.5 W/m 2 .Figure 7a,b indicate that the large negative biases in DLR occurred at large VZA, which happened at high IWVs.Considering that the Wang-Liang algorithm was developed over the North American continent, 38 CLAR profiles in this region were analyzed.Results in the North American continent were similar to those worldwide except that the best result was at δT s,a = 5 K, and RMSE varied from 14.9 to 163.4 W/m 2 with VZA increase at δT s,a = 10 K.
dependent on δTs,a and VZA.The Tang-Li algorithm had the best result at δTs,a = 0 K and performed worse when the absolute value of δTs,a became larger.This is because Ta was considered equal to Ts in the algorithm development [30].For the Wang-Liang algorithm, RMSE and bias were minimal at δTs,a = 10 K and increased at other δTs,a values.Though MODIS-derived Ta and Ts were used in the MODTRAN simulation during their algorithm development [31], the coefficients were not regressed for each δTs,a value separately.Therefore, they represented the average δTs,a of all training profiles.When VZA ≤ 40°, DLR RMSE is 17.0 to 23.8 W/m 2 for the Tang-Li algorithm at δTs,a = 0 K and 18.7 to 31.4 W/m 2 for the Wang-Liang algorithm at δTs,a = 10 K.When VZA was larger, RMSEs of the two algorithms increased rapidly with VZA, especially for Wang-Liang algorithm, with RMSE reaching 151.5 W/m 2 .Figure 7a and 7b indicate that the large negative biases in DLR occurred at large VZA, which happened at high IWVs.Considering that the Wang-Liang algorithm was developed over the North American continent, 38 CLAR profiles in this region were analyzed.Results in the North American continent were similar to those worldwide except that the best result was at δTs,a = 5 K, and RMSE varied from 14.9 to 163.4 W/m 2 with VZA increase at δTs,a = 10 K.The Gupta2010 and Zhou-Cess algorithms performed well, with RMSEs of 10.3 and 11.2 W/m 2 (Figure 7).The new parameterization had a smaller RMSE (about 5 W/m 2 ) than the two algorithms for VZA ≤ 50 • (comparing Figures 6 and 7) but had a slightly larger RMSE (about 1 to 9 W/m 2 ) than the other two algorithms at larger VZAs.However, when the two algorithms were applied to actual satellite data, DLRs from them were greatly influenced by the uncertainty of T a , which was difficult to retrieve accurately using satellites, especially for those without the TIR sounding channels.

Sensitivity Analysis of Input Parameters
The sensitivity of three main inputs parameters, IWV, BT, and δT s,a , were analyzed for three various climate types.As shown in Table 4, Cases 1, 2, and 3 were at the conditions of low-BT and low-IWV, high-BT and low-IWV, and high-BT and high-IWV, respectively.The three cases represent a dry and cold climate such as a polar region, a hot and dry climate such as a desert, and a wet and hot climate such as a tropical ocean, respectively.The sensitivity was defined as the DLR error versus the error of an input parameter while keeping the others fixed.The definition of each case and the error range of input parameters are given in Table 4, and the results of sensitivity analysis are displayed in Figure 8.

Sensitivity to Error in Temperature Difference between Ground and Air
Figures 8d and 8e show the DLR sensitivity to the error of δTs,a when the ground was 10 K cooler or warmer than the surface air.The DLR error linearly increased with the absolute error δTs,a and reached its positive (or negative) maximum when the δTs,a error was −10 K (or 10 K) for δTs,a = −10 K (or 10 K) because δTs,a was set within ±20 K in the DLR calculation.DLR was more sensitive to δTs,a for high-BT and low-IWV conditions (Case 2) than other cases.Furthermore, the DLR error difference between 0° and 50° was only 1 to 2 W/m 2 for Cases 1 and 2, but it reached 17.3 W/m 2 for Case 3. The figure indicated that if δTs,a was ignored in DLR estimation, which meant the δTs,a error was 10 K (or −10 K) for δTs,a = −10 K (or 10 K), DLR would be underestimated by 27.3, 48.3, and 25.9 W/m 2 (or overestimated by 23.1, 40.9, and 25.2 W/m 2 ) for Cases 1, 2, and 3 for VZA = 0°.Therefore, we suggest that δTs,a should be considered for the bare soil, sand, and dry grass surfaces, especially for deserts where higher Ts and lower Ta and high BT and low IWV frequently exist in the daytime.Also, the δTs,a is non-negligible for an ice surface where Ts is much lower than Ta.
The above sensitivity analysis indicates that DLR accuracy in actual application is greatly affected by uncertainties of IWV, BT, and δTs,a.When RE of IWV is within 50%, it produces a very  (a,b) are for H = 0 and 3.0 km, respectively.Figures (d,e) represent the conditions of δT s,a = −10 and 10 K, respectively.

Sensitivity to Water Vapor Content Error
Figure 8a,b shows the DLR error to IWV relative error (RE) at elevation of 0 and 3 km, respectively.We found that satellite retrieved IWVs were sometimes 10 times smaller or larger than the actual data at very small IWVs, and IWV RE was relatively small at large IWVs.Therefore, we set different IWV RE for different conditions; they were set within ± 60% for high IWV conditions (Case 3) and within −90% to 400% for low IWVs (Cases 1 and 2).DLR error was small for IWV RE within ±50% under low-IWV conditions, i.e., −7.8 to 5.4 W/m 2 and −13.7 to 9.5 W/m 2 for Cases 1 and 2 when elevation = 0 km and VZA = 0 • .However, when the IWV RE range was −90 to 400%, DLR had a larger error, −24.5 to 34.8 W/m 2 and −43.5 to 61.7 W/m 2 for Cases 1 and 2. DLR was most sensitive to IWV uncertainty under high-IWV conditions (Case 3), with a DLR error of −74.9 to 58.6 W/m 2 for IWV RE within 60%.Moreover, DLR was slightly more sensitive to IWV at larger VZAs, but the DLR error difference between VZA = 0 • and 50 • could be considered negligible (about 0 to 2 W/m 2 ) unless the IWV or the IWV error was very large.Compared with low elevation, IWV uncertainty produced much larger DLR errors when IWV RE was large for low-IWV conditions at high elevation (Figure 8b).DLR error was −41.9 to 42.3 W/m 2 and −74.3 to 75.1 W/m 2 for Cases 1 and 2, with a DLR RE of −23.8 to 24.1% for both conditions.Nevertheless, DLR errors for high IWV were very similar to those at low altitude.

Sensitivity to Brightness Temperature Error
The BT error was from an instrument calibration error and spatial mismatch between the ground station and the satellite pixel.Compared to the former error, the latter error was very large in rugged areas and under cloud contamination conditions.As shown in Figure 8c, the DLR error showed a linear relationship with the BT error when it was −5 to 5 K.For altitude = 0 km and VZA = 0 • , the DLR error of Case 1 was −13.9 to 14.8 W/m 2 and RE was −7.5 to 7.9% (not shown).The DLR error was more sensitive to BT at high-BT conditions at −21.5 to 22.5 W/m 2 for Case 2 and −29.5 to 31.0 W/m 2 for Case 3.

Sensitivity to Error in Temperature Difference between Ground and Air
Figure 8d,e show the DLR sensitivity to the error of δT s,a when the ground was 10 K cooler or warmer than the surface air.The DLR error linearly increased with the absolute error δT s,a and reached its positive (or negative) maximum when the δT s,a error was −10 K (or 10 K) for δT s,a = −10 K (or 10 K) because δT s,a was set within ±20 K in the DLR calculation.DLR was more sensitive to δT s,a for high-BT and low-IWV conditions (Case 2) than other cases.Furthermore, the DLR error difference between 0 • and 50 • was only 1 to 2 W/m 2 for Cases 1 and 2, but it reached 17.3 W/m 2 for Case 3. The figure indicated that if δT s,a was ignored in DLR estimation, which meant the δT s,a error was 10 K (or −10 K) for δT s,a = −10 K (or 10 K), DLR would be underestimated by 27.3, 48.3, and 25.9 W/m 2 (or overestimated by 23.1, 40.9, and 25.2 W/m 2 ) for Cases 1, 2, and 3 for VZA = 0 • .Therefore, we suggest that δT s,a should be considered for the bare soil, sand, and dry grass surfaces, especially for deserts where higher T s and lower T a and high BT and low IWV frequently exist in the daytime.Also, the δT s,a is non-negligible for an ice surface where T s is much lower than T a .
The above sensitivity analysis indicates that DLR accuracy in actual application is greatly affected by uncertainties of IWV, BT, and δT s,a .When RE of IWV is within 50%, it produces a very large error into DLR retrieval for large IWVs (i.e., 5 cm), and a relatively smaller DLR error for small IWVs (i.e., 0.5 cm).However, IWV RE is usually much larger than 50% under small IWV conditions and therefore leads to a large error of DLR.Moreover, an IWV error leads to larger error of DLR at high elevations than that at low elevations.The BT error caused by the spatial mismatch of pixel and field in rugged areas and by the cloud contamination causes a significant error to DLR retrievals; 10 K error in δT s,a causes an error of about 20 to 40 W/m 2 to DLR.Therefore, δT s,a should be considered under the conditions that T s is very different from T a (δT s,a = −10 K or 10 K), which is related to bare soil, sand, dry grass, and ice surfaces.For the input parameters of IWV and BT, DLR is more sensitive to their errors under high BT and high IWV conditions than under other conditions, whereas for δT s,a , DLR is most sensitive to errors under high BT and low IWV conditions.

MODIS data
The collection 6.0 MODIS products were used, including MOD021KM, MOD03, MOD35_L2 [43], MOD05_L2 [44,45], and MOD07_L2 [46].The input parameters of the new parameterization and the four algorithms for comparison are shown in Table 5.We used near-infrared water vapor from MOD05_L2 in the daytime and the infrared IWV in the nighttime when the near-infrared water vapor data were unavailable.MOD07_L2 furnishes T s , surface pressure, and temperature and moisture profiles at 20 pressure levels.Surface T a was acquired by interpolating MODIS temperature at the lowest pressure level to the surface altitude, assuming a temperature lapse rate of 6.5 K/km [18].Then, δT s,a and the mean temperatures of the lowest two layers for the Gupta2010 algorithm were calculated.The lowest two atmospheric layers were defined as the surface to 780 hPa and 780 to 700 hPa when surface pressure was ~1000 hPa, whereas both layers were ~100 hPa thick over high-altitude regions.The lowest layer was defined as at least 25 hPa above the surface to guarantee that the temperature lapse rate was applied to a layer at least 25 hPa in thickness [47].Li et al. [48] found that the MODIS's profile lacked lower-atmospheric data for high-altitude regions.To evaluate the performance of different atmospheric products in DLR calculation, National Centers for Environmental Prediction (NCEP) data with higher vertical resolution were used and compared with MODIS.The NCEP FNL (Final) Operational Global Analysis data were prepared operationally every six hours at 1 • × 1 • resolution.The NCEP data provided several parameters required by the study, including surface pressure, total water vapor content, 2-m air temperature, surface temperature, and atmospheric profiles.The atmospheric profiles included temperature and relative humidity at 26 levels from 1000 to 10 hPa (1000, 975, 950, 925, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100, 70, 50, 30, 20, and 10 hPa).
The atmospheric parameters for DLR retrieval, including T a , IWV, T s , T e , and atmospheric profiles, were obtained from NCEP, and the other inputs were unchanged.The profiles were obtained by interpolation from the 1 • × 1 • original grid into a 1 km × 1 km mesh.We used the method proposed by [47] to extract the satellite pixel profile and other parameters from the NCEP data, including temporal, spatial, and vertical interpolation.First, two NCEP files representing two adjacent forecast times were loaded and interpolated to the satellite time using linear interpolation with time difference.Then, the profiles at satellite pixel were extracted by bilinear interpolation from the four grid corners surrounding a specific location.Finally, the T a and humidity at the surface level were bilinearly interpolated from the two closest layers to the surface.T s was obtained after temporal and spatial interpolation, whereas T a , T e , and IWV were calculated from the final atmospheric profile.

Ground Measurements
DLR measurements were collected from 14 sites around the globe belonging to either the Chinese Watershed Allied Telemetry Experimental Research (WATER; [49]), the Atmospheric Radiation Measurement (ARM) Program, the international Baseline Surface Radiation Network (BSRN), or the Ameriflux.The information about site location, period, instrument, and measurement uncertainty are presented in Table 6.The sites are in various climate types including polar regions with dry and cold climate, tropical marine climate with high temperature and IWV, and arid regions with low IWV and large δT s,a .The arid sites are marked with asterisks, including the desert surface and the sites that were considered as arid by [22].All sites are a homogeneous surface type, and most of them are in flat and rural regions.The homogeneity of the sites was analyzed using ASTER global digital elevation model (DEM) data with 30-m resolution.The standard deviation of elevation within a 3 km × 3 km is 4-40 m for most sites; it is about 347 m for Dome C (DOM), which is relatively small compared with the site altitude and has little influence on DLR calculation.Therefore, all sites can be regarded as homogeneous in the validation.
The radiation data provided by WATER, ARM, BSRN, and Ameriflux sites were in 10-min, 1-min, 1-min and 30-min average, respectively.For the sites from the former three networks, the data closest to the satellite overpass time were chosen, while for the Ameriflux sites, in-situ DLR was interpolated from the two data points that were closest to the satellite pass time.Besides surface radiations, the three sites of ARM and South Pole (SPO) of BSRN provided discontinuous measurements of IWV.One year of data were used for most sites, whereas field measurements were discontinuous during 2010 at DOM, and only summer and winter data were used for NSA-C1 and SPO because the Terra-MODIS passes overhead as many as eight times per day at these sites.
DLR was measured by Eppley Precision Infrared Radiometer (PIRs), CG3 pyrgeometers of Kipp & Zonen CNR1 Net Radiometers, or Kipp & Zonen CG 4 pyrgeometers at these sites.For the ARM network, the Eppley PIRs were shaded and ventilated to prevent excessive heating of the dome by the sun and were equipped with the Solar Infrared Radiation Station (SIRS) system at SGP-C1 and with the Sky Radiation (SKYRAD) system at Barrow North Slope of Alaska (NSA-C1) and Manus Island Tropical Western Pacific (TWP-C1).Their DLR measurement uncertainties were provided by Stoffel [50,51].PIR accuracy for the BSRN DLR reached 10 W/m 2 through improvement of its calibration [52].DLR uncertainties used the factory accuracies for other sites.The spectral ranges of the CG3, PIR, and CG4 were 5-50, 4.0-50, and 4.5-42 µm, respectively.Although the spectral response functions of the pyrgeometers did not cover the whole infrared spectrum, their readings were calibrated to the total range of longwave emission [53].

Validation Result of MODIS-Derived DLRs
The pixels that the MODIS cloud mask considered to be clear at 95% and 99% probability were selected, and those lacking information about atmospheric profile or water vapor were removed as cloudy.Because the MODIS cloud mask may have substantial uncertainty in polar regions, we determined clear-sky conditions from field measurements at DOM, SPO, and NSA-C1 using the dual-threshold method proposed by [54].The clear-sky threshold and standard deviation threshold of DLR for each month were selected according to the time series of DLR measurements and the clear-sky DLRs simulated by MODTRAN and in-situ atmospheric profiles.The measurements for DLR that were greater than a monthly clear-sky DLR threshold or DLR standard deviation over a 20-min period greater than the monthly standard deviation threshold were removed as cloudy conditions.Furthermore, the pixels with unreasonably low TOA radiance at sites Flagstaff Managed Forest (FMF) and US-GLE Wyoming (GLEES) and those with very low IWV at TWP-C1 were excluded as cloud contaminated.
The DLRs of MODIS clear-sky pixels were calculated and compared with in-situ measurements.As shown in Table 7, the validation results were arranged into three region types-the sites belonging to the arid region, the sites belonging to the high-altitude region (elevation > 2000 m), and other sites belonging to ordinary regions.The statistical results for each site and their averages for each region type are presented.The scatter plots of estimated and ground-measured DLRs for the three region types are shown in Figures 9-11   Table 7. Error statistics for MODIS-derived clear-sky DLR using different algorithms.RMSE, bias, and the square of the correlation coefficients (R 2 ) of all sites and the average values of each type of region are given.RMSE and bias are in units of W/m 2 .The best results are highlighted in bold, and the worst results are marked with asterisks.

NewParam
Tang    DLRs accuracies ranged from 7.8 to 35.7 W/m 2 at the high-altitude sites.As shown in Figure 11a, DLRs had large positive errors at the large DLRs, which mainly occurred at GLEES, and had negative errors at lower DLRs that were in the winter of SPO.Furthermore, some overestimations at Arou (AR) and overestimations and underestimations at FMF were found.The errors at GLEES, AR, and FMF were mainly caused by the uncertainty of IWV.We found MODIS IWV had some values

Performance of New Algorithm
Because MODIS-derived δT s,a may have large uncertainties, we used δT s,a only at the sites in arid and polar regions and the SGP-C1 site where Ts were largely different from Ta and the uncertainties.Specifically, we used dry grass-only coefficients for the arid sites, and all surfaces-average coefficients for other sites.Table 7 indicates that the new algorithm had good results at most sites with accuracies within 25 W/m 2 .The mean RMSEs of DLRs were 20.5, 23.3, and 22.1 W/m2 for ordinary, arid, and high-altitude regions, respectively.
As shown in Figure 9a, the new parameterization performed well under most circumstances in ordinary regions but produced some negative errors at the lower DLRs (≤200 W/m 2 ), which were at NSA-C1, and had some positive and negative errors at the higher DLRs (~400 W/m 2 ), which were at TWP-C1.We further analyzed the error sources by comparing satellite and field measured parameters at NSA-C1 and TWP-C1.Because field measured IWVs were discontinued, only parts of the samples were analyzed.The DLR underestimations of NSA-C1 mainly occurred in the winter season when MODIS IWV had negative errors up to −50% compared to field measurements.The other error sources included uncertainty in δT s,a and the fact that the algorithm could not accurately simulate DLRs for some polar profiles, as mentioned in Section 3.1.The DLR underestimation and overestimation at TWP-C1 were mainly caused by the uncertainty of the MODIS IWV error (from −1.9 to 2.0 cm), whereas uncertainty caused by δT s,a was relatively small.DLR RMSEs were 19.5-26.8W/m 2 for the arid sites.Using δT s,a obviously improved the results at these sites (not shown), especially during daytime.Slight overestimations appeared at the larger DLRs (Figure 10a), which mainly happened during daytime.Also, slight underestimations during nighttime were found (not shown).We analyzed the δT s,a values at sites Huazhaizi (HZZ) and Desert Rock (DRA), and found the in-situ δT s,a had large positive values (10 to 20 K at HZZ and 5 to 20 K at DRA) during daytime and slightly negative values (−5 to 0 K) during nighttime.However, the ranges of MODIS δT s,a errors during daytime and nighttime were −5 to 5 K and 5 to 10 K at HZZ and were −10 to 5 K (mainly negative errors) and 0 to 5 K at DRA.Therefore, the retrieved DLRs had larger negative errors during nighttime at HZZ and larger positive errors during daytime at DRA.We recalculated the DLRs using in-situ δT s,a and found the uncertainty caused by δT s,a was within 20 W/m 2 at both sites.
DLRs accuracies ranged from 7.8 to 35.7 W/m 2 at the high-altitude sites.As shown in Figure 11a, DLRs had large positive errors at the large DLRs, which mainly occurred at GLEES, and had negative errors at lower DLRs that were in the winter of SPO.Furthermore, some overestimations at Arou (AR) and overestimations and underestimations at FMF were found.The errors at GLEES, AR, and FMF were mainly caused by the uncertainty of IWV.We found MODIS IWV had some values much greater than those on adjacent days during nighttime at AR and GLEES.When the pixels with IWV ≥2.1 cm were removed, the DLR accuracies improved from 6.3 ± 35.7 to −3.2 ± 25.7 W/m 2 at GLEES.As for the SPO site, the underestimations were caused by IWV uncertainty, the algorithm limitation in polar winter, and strong temperature inversions.The mean relative error of MODIS-derived IWV was −71.9%, and the thick inversion layer with temperature lapse rate exceeding −20 K/km was frequently found at SPO during winter.

Comparison with Other Algorithms
The new algorithm was compared with the four algorithms in Section 2.1.Considering that the coefficients of Wang-Liang algorithm were derived over the North American continent and cannot represent the atmosphere in tropical and Antarctic regions, its coefficients were calibrated at sites DOM, SPO, and TWP-C1 to compare these algorithms equally.
In the ordinary regions, the RMSE of MODIS DLRs from the new algorithm was 1.0 and 1.5 W/m 2 smaller than those from Wang-Liang and Gupta2010 algorithms and was 13.8 and 5.4 W/m 2 smaller than those from Tang-Li and Zhou-Cess algorithms.Though Wang-Liang and Gupta2010 had similar performance to the new algorithm (Figure 9), Gupta2010 had large positive errors at SGP, and Wang-Liang slightly underestimated the DLRs at Yingke (YK) and overestimated the DLRs at NSA-C1 (Table 7).For the Zhou-Cess algorithm, the DLRs at the lower end were overestimated, and those at the higher end were underestimated for the DLRs within 200-400 W/m 2 (at YK and SGP-C1), while the DLRs ≥400 W/m 2 (at TWP-C1) were underestimated.Tang-Li significantly overestimated the DLRs at all ranges.
For the sites in arid regions, the RMSE from the new algorithm was slightly larger than those from Zhou-Cess and Gupta2010 (~0.5-1.0W/m 2 ), and much smaller than those from the Tang-Li and Wang-Liang algorithms (~8.5-20.3W/m 2 ).Though it performed slightly better, Zhou-Cess showed underestimations overall (Figure 10b), whereas Gupta2010 had obvious negative biases at HZZ and SBO (Figure 10c).The Tang-Li algorithm significantly overestimated DLRs, whereas the Wang-Liang algorithm obviously underestimated DLRs (Figure 10d,e).These findings may be because δT s,a was not considered by either algorithm.The large positive biases of the Tang-Li algorithm mainly occurred during daytime when large δT s,a values were found at these sites.The performance of the Wang-Liang algorithm was consistent with the conclusion of Section 3.1 that DLRs had negative errors for δT s,a < 10 K.
In the high-altitude regions, the new algorithm performed slightly better than the Gupta2010, Wang-Liang, and Tang-Li algorithms (with mean RMSE decreased 1.0 to 3.7 W/m 2 ) and much better than the Zhou-Cess algorithm (with mean RMSE decreased 11.4 W/m 2 ).For the DLR within 200-400 W/m 2 that occurred mainly at sites AR, FMF, and GLEES, the Zhou-Cess and Gupta2010 algorithms had larger positive and negative errors than the new algorithm.Tang-Li slightly overestimated the DLRs overall, while Wang-Liang performed well but had some negative errors (Figure 11).Besides the uncertainty of MODIS IWV, the results of the Gupta2010 and the Zhou-Cess algorithms were also influenced by the uncertainty of MODIS-derived T a .For example, the two algorithms greatly underestimated DLRs at Arou, where we found that T a from MODIS was ~10 K lower than those from the ground measurement.For the lower DLRs in polar regions, Tang-Li overestimated them while the Wang-Liang and the Gupta2010 algorithms underestimated the winter data among them (Figure 11).However, the Zhou-Cess algorithm greatly overestimated the DLRs in polar regions.This may be because the in-situ DLRs and satellite cloud fractions of the two sites were used during the algorithm development [25], and some data of cloudy pixels may not have been discriminated and were fitted as clear sky.
In summary, compared with the Tang-Li and the Wang-Liang algorithms, the new algorithm is based on the Stefan-Boltzmann law and therefore has more reasonable physical meaning.It also performs better over arid regions or other conditions with large δT s,a because the difference between air and surface temperature is accounted for in the algorithm.Compared with the Zhou-Cess and the Gupta2010 algorithms, the new algorithm can avoid the underestimations or overestimations in DLR caused by inaccurate air temperature over high-altitude regions.

Evaluation of Different Atmospheric Products for DLR Retrievals
We found the same problem as Li et al. [48] that MODIS's profile lacked lower-atmospheric data at high-altitude sites AR and YK.NCEP data with higher vertical resolution were used and compared with MODIS for the high-altitude sites.In this section, the atmospheric parameters for DLR retrieval, including T a , IWV, T s , and atmospheric profiles, were obtained from NCEP, and the other inputs were unchanged.

Difference between MODIS-and NCEP-derived Atmospheric Parameters
Figure 12 shows a comparison of MODIS-and NCEP-derived atmospheric profiles at sites AR, GLEES, and SPO.For the AR site, because the surface pressure from MODIS was much smaller than the observed value, the atmospheric information of a thickness ~50 hPa above the surface was missing.Since the surface T a and infrared IWV was obtained from the MODIS atmospheric profile, the information deficiencies in the lower troposphere may have led to the underestimations of T a and nighttime IWV.Conversely, the surface pressure from MODIS was ~60 hPa larger than the observed values at the GLEES site, which may have caused IWV overestimation.The surface-level T a that interpolated from NCEP data was closer to observed values than those from MODIS at the two sites.As for SPO, the surface pressure was similar to the observed value, but the T a interpolated from MODIS lowest atmospheric level by assuming a temperature lapse rate of 6.5 K/km had a larger error than the NCEP-derived T a .MODIS lowest atmospheric level by assuming a temperature lapse rate of 6.5 K/km had a larger error than the NCEP-derived Ta.  Figure 13 shows the histograms of IWV and bias in Ta from the two products at these sites.MODIS-derived Ta had a median bias of about −10 and −4 K for daytime and nighttime at AR, while MODIS lowest atmospheric level by assuming a temperature lapse rate of 6.5 K/km had a larger error than the NCEP-derived Ta.  Figure 13 shows the histograms of IWV and bias in Ta from the two products at these sites.MODIS-derived Ta had a median bias of about −10 and −4 K for daytime and nighttime at AR, while Figure 13 shows the histograms of IWV and bias in T a from the two products at these sites.MODIS-derived T a had a median bias of about −10 and −4 K for daytime and nighttime at AR, while the NCEP-derived T a was better during the daytime.The distribution of bias in T a from NCEP was more centralized toward zero compared to that of MODIS at the GLEES site.The MODIS IWV had some large IWV values at AR and GLEES, which were unreasonable, whereas these very large values were smoothed by NCEP data.For the SPO site, we found the number of underestimated T a from MODIS was deduced, and the IWVs negative biases during winter that were mentioned in Section 4.2 were decreased after using NCEP data compared to those from MODIS.

Difference between MODIS-and NCEP-derived DLRs
DLRs at the high-altitude sites were calculated by the NCEP parameters using the new, the Zhou-Cess, and the Gupta2010 algorithms.Table 8 indicates that using NCEP parameters distinctly improved the results, especially for the Gupta2010 and the Zhou-Cess algorithms.The DLR underestimations at SPO during winter and at AR, as well as the overestimations at GLEES, were substantially corrected.DLR RMSEs decreased by 5.5 and 19.1 W/m 2 at AR and GLEES for Zhou-Cess and decreased by 2.4, 19.0 and 6.7 W/m 2 at AR, SPO, and GLEES for Gupta2010.Because only IWV was improved, the new algorithm was less affected by different sources of atmospheric data than the other two algorithms.The DLR RMSEs decreased at SPO and GLEES (3.1 and 7.8 W/m 2 ) but slightly increased at AR and FMF. Figure 14 displays the scatter points between retrieved and observed DLRs and the distributions of DLR biases from NCEP and MODIS, respectively.Comparing Figure 14 with Figure 11, the dispersed points were diminished for the three algorithms, and the underestimated points at the low DLRs were corrected for Gupta2010.Figure 14d-i indicates that errors from NCEP had narrower distributions than those from MODIS.Furthermore, we compared the two products in other regions and found that using NCEP-derived parameters appreciably improved results of these algorithms under extreme climatic conditions, such as during the winter of NSA-C1 and at site TWP-C1.NCEP data are not recommended for other regions because the coarse resolution causes great uncertainties in atmospheric parameters.Comparing Figure 14 with Figure 11, the dispersed points were diminished for the three algorithms, and the underestimated points at the low DLRs were corrected for Gupta2010.Figure 14 (d)-(i) indicates that errors from NCEP had narrower distributions than those from MODIS.Furthermore, we compared the two products in other regions and found that using NCEP-derived parameters appreciably improved results of these algorithms under extreme climatic conditions, such as during the winter of NSA-C1 and at site TWP-C1.NCEP data are not recommended for other regions because the coarse resolution causes great uncertainties in atmospheric parameters.

Conclusions
In this study, the Yu2013 algorithm that parameterized clear-sky DLR using BT and IWV was improved and extended to MODIS data.The improved parameterization avoids considerable DLR underestimations under high IWV conditions by changing the functional form of water vapor from square-root to logarithmic.An evaluation based on a testing dataset generated from global atmospheric profiles and MODTRAN simulations indicates that the new parameterization was not restricted by sample number and could greatly improve results at large IWV compared with the original parameterization.The total RMSE of DLR decreased by ~1.9 to 3.1 W/m 2 , and DLR accuracy under high IWV (>3.0 cm) improved by 0.5%-16.9%.
Accuracy evaluation using the simulated testing dataset shows that the algorithm gives accuracies between 5.4 and 20.8 W/m 2 when the input parameters are perfectly accurate.A sensitivity analysis indicates that DLR accuracy is greatly affected by uncertainties of IWV, BT, and δT s,a .DLR error is more sensitive to the error in IWV at large IWVs and at small IWVs at high altitudes, to the error of BT at high-BT conditions, and to the error of δT s,a under high-temperature and low-IWV conditions than it is to other conditions.The new algorithm was compared with four state-of-the-art algorithms.The two algorithms parameterized with TOA radiance produced large positive or negative errors, especially for large DLRs when T s was much higher or lower than T a .The new algorithm accounted for δT s,a and made improvements at high IWV, thereby achieving better results than the two algorithms.The new algorithm also gives slightly better results to the two atmospheric parameter-based algorithms for VZA ≤50 • , whereas the latter two algorithms would be greatly affected by the uncertainties of IWV and air and surface temperatures in actual applications.
The new algorithm was applied to MODIS Terra data and extensively validated using one year's ground data from 14 stations around the globe.The algorithm produced favorable results at most sites with accuracies within 25 W/m 2 .However, the algorithm performed poorly under extreme climatic conditions, such as winter in polar regions with a dry and cold climate and in the tropical western Pacific with a warm and moist climate.The poor performance for these conditions was mainly caused by the uncertainty of MODIS-derived IWV and δT s,a and the algorithm limitations under extreme climate.The proposed algorithm had much better results than the TOA radiance-based algorithms in dry-arid regions or other conditions in which T s was substantially different from T a because δT s,a was considered.The algorithm performed 1.0 to 11.4 W/m 2 better than atmospheric parameter-based algorithms in high-altitude regions because DLRs from the latter two algorithms were greatly affected by the inaccurate air temperature.
To evaluate the effect of different atmospheric datasets on DLR accuracy, MODIS and NCEP atmospheric products were compared for DLR estimation.The NCEP-derived atmospheric parameters were more accurate than MODIS data in high-altitude regions because the lowest layers from MODIS profiles sometimes were above or below the actual surface there, which caused the underestimations or overestimations of IWV and T a .In addition, using NCEP data could avoid the very large positive/negative IWV errors and T a errors from MODIS.Therefore, using NCEP data could clearly improve DLR estimation in such regions, especially for the atmospheric parameter-based algorithms (RMSE decreased by 2.2 to 19.1 W/m 2 ).
The study demonstrates that the new parameterization works well in most situations and overcomes the shortcomings of current TOA radiance-based and atmospheric parameter-based algorithms.In future work, the new algorithm will be applied to geostationary satellite data, such as Meteosat-9/10 and GOES13, to obtain diurnal changes in DLR, whereby acquiring accurate IWV and δT s,a is most important during the DLR estimations.Moreover, the algorithm needs further modification to improve performance under extreme climatic conditions.

Figure 1 .
Figure 1.Brightness temperature of Moderate Resolution Imaging Spectroradiometer (MODIS) channel 31 versus atmospheric effective temperature for various δTs,a values.Atmospheric effective temperature was defined by the DLR and atmospheric effective emissivity that was simulated using MODTRAN and 875 TIGR2002 profiles.

Figure 2 .
Figure 2. Flowchart of DLR modeling process and DLR derivation.

Figure 1 .
Figure 1.Brightness temperature of Moderate Resolution Imaging Spectroradiometer (MODIS) channel 31 versus atmospheric effective temperature for various δT s,a values.Atmospheric effective temperature was defined by the DLR and atmospheric effective emissivity that was simulated using MODTRAN and 875 TIGR2002 profiles.

Figure 1 .
Figure 1.Brightness temperature of Moderate Resolution Imaging Spectroradiometer (MODIS) channel 31 versus atmospheric effective temperature for various δTs,a values.Atmospheric effective temperature was defined by the DLR and atmospheric effective emissivity that was simulated using MODTRAN and 875 TIGR2002 profiles.

Figure 2 .
Figure 2. Flowchart of DLR modeling process and DLR derivation.

Figure 2 .
Figure 2. Flowchart of DLR modeling process and DLR derivation.

Figure 3 .
Figure 3. Spectral emissivity curves of seven surface types from ASTER used in MODTRAN simulation.

Figure 3 .
Figure 3. Spectral emissivity curves of seven surface types from ASTER used in MODTRAN simulation.

7 Figure 4 .
Figure 4. Atmospheric emissivity (ε31) defined by Equation (3) using brightness temperature (BT) of MODIS 31st channel, plotted as functions of IWV at the condition of conifer surface, H = 1.435km, view zenith angle (VZA) = 0 o , and δTs,a =0K.The black and gray dots represent TIGR2002 and ECWMF profiles, respectively.Only the TIGR2002 profiles were used in regression in Figure (a), and all the TIGR2002 and the ECWMF profiles were used in Figure (b).Poly2 and poly3 represent quadratic and cubic function, respectively.

Figure 4 .
Figure 4. Atmospheric emissivity (ε 31 ) defined by Equation (3) using brightness temperature (BT) of MODIS 31st channel, plotted as functions of IWV at the condition of conifer surface, H = 1.435 km, view zenith angle (VZA) = 0 • , and δT s,a = 0 K.The black and gray dots represent TIGR2002 and ECWMF profiles, respectively.Only the TIGR2002 profiles were used in regression in Figure (a), and all the TIGR2002 and the ECWMF profiles were used in Figure (b).Poly2 and poly3 represent quadratic and cubic function, respectively.

Figure 5 .
Figure 5. DLRs calculated from testing dataset using the Yu2013 (a and c) and new parameterizations (b and d), respectively.The algorithms coefficients were deriving from TIGR2002 profiles.Figures (a)-(b) show the comparison between the estimated and the MODTRAN-simulated DLRs.Figures (c)-(d) show the residuals of the predicted DLR, which vary with water vapor content for different VZAs, assuming δTs,a = 0 K.

Figure 5 .
Figure 5. DLRs calculated from testing dataset using the Yu2013 (a,c) and new parameterizations (b,d), respectively.The algorithms coefficients were deriving from TIGR2002 profiles.Figures (a,b) show the comparison between the estimated and the MODTRAN-simulated DLRs.Figures(c,d) show the residuals of the predicted DLR, which vary with water vapor content for different VZAs, assuming δT s,a = 0 K.

Figure 6 .
Figure 6.Figures (a-f) illustrate the RMSE and Mean-Bias-Error (MBE) between the actual and estimated DLRs versus VZA and δT s,a for conifer surface and for NewParam (a,d), Tang-Li algorithm (b,e), and Wang-Liang algorithm (c,f), respectively.

Figure 7 .
Figure 7.Comparison of the clear-sky flux calculated from DLR algorithms with those from MODTRAN simulation based on the testing dataset.Figures (a) to (d) represent Tang-Liang, Wang-Li, Gupta2010, and Zhou-Cess algorithms, respectively, and δTs,a is 0 K and 10 K in Figures (a) and (b).

Figure 7 .
Figure 7.Comparison of the clear-sky flux calculated from DLR algorithms with those from MODTRAN simulation based on the testing dataset.Figures (a) to (d) represent Tang-Liang, Wang-Li, Gupta2010, and Zhou-Cess algorithms, respectively, and δT s,a is 0 K and 10 K in Figures (a,b).

Figure 8 .
Figure 8. DLR error versus the error of input parameters such as IWV (a,b), brightness temperature (c), and δT s,a (d,e) under different climate conditions.Cases 1, 2, and 3 represent the conditions of BT = 260 K and IWV = 0.5 cm, BT = 300 K and IWV = 0.5 cm, and BT = 300 K and IWV = 5.0 cm, respectively.Figures(a,b) are for H = 0 and 3.0 km, respectively.Figures(d,e) represent the conditions of δT s,a = −10 and 10 K, respectively. .

Figure 9 .
Figure 9.Comparison between MODIS-estimated DLRs and in-situ measurements at the sites in ordinary regions.Figures (a) to (e) represent the results from the new parameterization (NewParam), Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.The color bar represents data density that ranges from 0 to 100 percent of the data.

Figure 9 .
Figure 9.Comparison between MODIS-estimated DLRs and in-situ measurements at the sites in ordinary regions.Figures (a) to (e) represent the results from the new parameterization (NewParam), Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.The color bar represents data density that ranges from 0 to 100 percent of the data.

Figure 9 .
Figure 9.Comparison between MODIS-estimated DLRs and in-situ measurements at the sites in ordinary regions.Figures (a) to (e) represent the results from the new parameterization (NewParam), Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.The color bar represents data density that ranges from 0 to 100 percent of the data.

Figure 10 .Figure 10 .
Figure 10.Same as Figure 9 but for the sites in the arid regions.Figures (a) to (e) represent the results from NewParam, Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.

Figure 11 .
Figure 11.Same as Figure 9 but for the sites in the high-altitude regions.Figures (a) to (e) represent the results from NewParam, Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.

Figure 11 .
Figure 11.Same as Figure 9 but for the sites in the high-altitude regions.Figures (a) to (e) represent the results from NewParam, Tang-Li, Wang-Liang, Zhou-Cess, and Gupta2010 algorithms.

Figure 12 .
Figure 12.Vertical distribution of air temperatures from the National Centers for Environmental Prediction (NCEP) and MOD07 profiles at sites Arou (AR), US-GLE Wyoming (GLEES), and South Pole (SPO).The solid line, dashed line, and asterisk represent MODIS, NCEP, and field-measured data, respectively.

Figure 13 .
Figure 13.Comparison of MODIS-derived atmospheric parameters with those from NCEP at the sites AR (left column), GLEES (middle column), and SPO (right column).The solid line and dashed line represent parameters from MODIS and NCEP, respectively.

Figure 12 .
Figure 12.Vertical distribution of air temperatures from the National Centers for Environmental Prediction (NCEP) and MOD07 profiles at sites Arou (AR), US-GLE Wyoming (GLEES), and South Pole (SPO).The solid line, dashed line, and asterisk represent MODIS, NCEP, and field-measured data, respectively.

Figure 12 .
Figure 12.Vertical distribution of air temperatures from the National Centers for Environmental Prediction (NCEP) and MOD07 profiles at sites Arou (AR), US-GLE Wyoming (GLEES), and South Pole (SPO).The solid line, dashed line, and asterisk represent MODIS, NCEP, and field-measured data, respectively.

Figure 13 .
Figure 13.Comparison of MODIS-derived atmospheric parameters with those from NCEP at the sites AR (left column), GLEES (middle column), and SPO (right column).The solid line and dashed line represent parameters from MODIS and NCEP, respectively.

Figure 13 .
Figure 13.Comparison of MODIS-derived atmospheric parameters with those from NCEP at the sites AR (left column), GLEES (middle column), and SPO (right column).The solid line and dashed line represent parameters from MODIS and NCEP, respectively.

Figure 14 .
Figure 14.The results of DLRs at the sites in high-altitude regions.Scatter plot between DLRs estimated from NCEP atmospheric parameters and observed DLR are shown in (a-c).Error histograms between observed and retrieved DLRs from MODIS (d-f) and from NCEP (g-i).Bias was computed as retrieved minus observed values.The left, middle, and right panels represent the results from the new parameterization, the Zhou-Cess, and the Gupta2010 algorithms, respectively.

Table 2 .
Error distribution statistics for IWV greater than 3.0 cm and δTs,a = 0K.

Table 2 .
Error distribution statistics for IWV greater than 3.0 cm and δT s,a = 0 K.

Table 3 .
The mean value of actual and calculated DLRs (DLRact and DLRcal), root mean square error (RMSE), bias, maximum error and minimum error between actual DLRs, and those calculated from training datasets for different surface types when δTs,a = 0K and VZA=0 o (unit: W/m 2 ).

Table 5 .
Input parameters and data sources used in this study.

Table 6 .
Description of site conditions.Sites used as examples of dry-arid regions are marked with asterisks.

Table 8 .
Same as Table7but for the results in the high-altitude regions.The clear-sky DLRs were estimated using NCEP atmospheric data.The differences between RMSE from NCEP and that from MODIS are displayed in the brackets.