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

Abstract

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/m² 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/m² for the new parameterization and 19.1 W/m² for other algorithms at high-altitude sites, as compared to MODIS atmospheric products.

1. 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]:

$$DLR=2\pi {\displaystyle {\int }_{v_1}^{v_2}{\displaystyle {\int }_0^1{I}_v(z=0,-\mu )\mu d\mu d\nu }},$$

(1)

where $I_v(z=0,-\mu )$ is the spectral radiance at the surface at a given local zenith angle $\theta$ and $\mu =\cos (\theta )$ and $\nu$ is the wavenumber. Under cloud-free conditions, multiple scatter can be omitted, thus spectral radiance can be expressed as:

$$I_{\lambda }(z=0,-\mu )={\displaystyle {\int }_0^{z_t}{B}_v(z^{\prime })\frac{\partial T_v(0,z^{\prime };-\mu )}{\partial z^{\prime }}d}{z}^{\prime },$$

(2)

where $z_t$ is the top of atmosphere, $B_v(z^{\prime })$ is the Planck function evaluated with the temperature at level $z^{\prime }$, and $T_v(0,z^{\prime };-\mu )$ is the monochromatic transmittance from surface to height $z^{\prime }$.

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² 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.

2. Development of Clear-Sky Algorithm

2.1. 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. Clear-sky downward longwave radiation (DLR) algorithms.

Authors	Abbreviation	Algorithm 1
Zhou et al. [25]	Zhou-Cess	$DLR=a_0+a_1\cdot SULW+a_2\cdot \ln (1+IWV)+a_3\cdot {[\ln (1+IWV)]}^2$
Gupta et al. [22]	Gupta2010	$\begin{array}{l}DLR=(a_0+a_{1}V+a_2{V}^2+a_3{V}^3)\cdot T_e^{3.7},V=\ln (10*IWV)\\ T_e=0.60T_s+0.35T_1+0.05T_2\end{array}$
Tang and Li [30]	Tang-Li	$\begin{array}{ll}DLR=& a_0+a_1{M}_{29}+a_2{M}_{34}+a_3{M}_{33}+a_4{M}_{36}+a_5{M}_{28}+\\ & a_5{M}_{28}+a_6{M}_{31},{{\displaystyle }}^{}{M}_i=\pi L_i\end{array}$,
Wang and Liang [31]	Wang-Liang	$\begin{array}{ll}DLR=& L_{Tair}(a_0+a_1{L}_{27}+a_2{L}_{29}+a_3{L}_{33}+a_4{L}_{34}+b_1{L}_{32}/L_{31}\\ & +b_2{L}_{33}/L_{32}+b_3{L}_{28}/L_{31}+c_{1}H)\end{array}$
Yu et al. [36]	Yu2013	$DLR=\sigma T_k^4(a_{0}IWV+a_1\sqrt{IWV}+a_2)$

¹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₁ and T₂), and H is surface elevation.

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 $\sigma {\epsilon }_a{T}_e^4$, 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:

$${\epsilon }_k=\frac{DLR}{\sigma T_k^4},$$

(3)

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]:

$${\epsilon }_k=a_{0}IWV+a_1\sqrt{IWV}+a_2,$$

(4)

where IWV has units of cm, and a_i (i = 0, 2) represent coefficients. Therefore, DLR is expressed as:

$$DLR=\sigma T_k^4(a_{0}IWV+a_1\sqrt{IWV}+a_2)$$

(5)

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].

2.2. 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₃₁) with a central wavelength of 10.8 μm, thus ε_k is expressed as ε₃₁. 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₃₁ was derived by regression analysis—actually, the relationship between ε₃₁ and IWV was derived. Finally, the algorithm was validated using simulated data and actual data.

2.2.1. 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 δ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⁻¹ at a resolution of 1 cm⁻¹. 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, ε₃₁ 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.

2.2.2. 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 ε₃₁, we intended to develop a robust formulation of atmospheric emissivity ε₃₁. Water vapor is the most important atmospheric gas contributing to clear-sky DLR, and the variation of other gases (such as CO₂ and O₃) have small effects on DLR [32,34]. Therefore, we used IWV to predict atmospheric emissivity ε₃₁. As shown in Figure 4, the atmospheric emissivity ε₃₁ 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 ε₃₁ 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:

$$DLR=\sigma T_{31}^4\left(a_0+a_{1}V+a_2{V}^2+a_3{V}^3\right)$$

(6a)

$$V=\sqrt{\ln (1+IWV})$$

(6b)

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 $\sqrt{\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 ε₃₁ fitted from three different formulations. Using $\sqrt{IWV}$ or ln(1 + IWV) will underestimate or overestimate the atmospheric emissivity ε₃₁ at very large IWV, especially when the number of training profiles is small (Figure 4a). On the contrary, using $\sqrt{\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 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² for δT_s,a = 0 K. The DLR statistical error distribution for IWV greater than 3.0 cm is shown in

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

Sensor Zenith View Angle	Range of Error (W/m2)	Yu2013	NewParam
VZA = 0°	[−20, 15]	99.5%	100%
	[−40, −20]	0.5%	0%
VZA = 30°	[−20, 15]	98.9%	99.5%
	[−40, −20]	1.1%	0.5%
VZA = 60°	[−20, 15]	73.8%	90.7%
	[−40, −20]	25.1%	8.2%
	[−200, −40]	1.1%	1.1%

. For small VZA (<60°), more than 99% of the DLRs had accuracy within 20 W/m², 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.

3. Algorithm Evaluation using Testing Dataset

3.1. 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. 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²).

Surface Types	Conifer	Soil	Wet Land	Concrete	Dry Grass	Sea Water	Fine Snow	All Types
DLRact	321.9	321.9	321.9	321.9	321.9	321.9	321.9	321.9
DLRcal	326.7	324.9	327.6	325.5	323.8	326.8	327.1	326.1
RMSE	5.6	5.4	6.7	5.1	8.4	5.7	6.0	6.2
Bias	1.8	−1.9	3.8	−0.6	−6.0	2.2	2.8	0.3
Max error	16.8	11.7	19.6	13.5	10.4	17.3	18.2	19.6
Min error	−24.0	−25.3	−23.3	−24.9	−26.3	−23.9	−23.7	−26.3

indicates that most surfaces had similar retrieval accuracy (RMSE = 5.1–8.4 W/m²), except of dry grass, which has a large negative error (bias= −6.0 W/m²). 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², respectively. Daytime DLRs had a larger positive error (~3.4 W/m²) than nighttime, and DLRs under temperature inversion conditions had a bias about −3.1 W/m² 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,d gives RMSE and bias variation with VZA for the conifer surface. RMSE increased with VZA from 5.6 to 20.2 W/m² for δT_s,a = 0 K, and the bias was 1.8 to −6.8 W/m². RMSE were 5.4 to 20.0 W/m² and 6.7 to 20.8 W/m² for δT_s,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 Figure 6 and Figure 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² for the Tang-Li algorithm at δT_s,a = 0 K and 18.7 to 31.4 W/m² 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². 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² with VZA increase at δT_s,a = 10 K.

The Gupta2010 and Zhou-Cess algorithms performed well, with RMSEs of 10.3 and 11.2 W/m² (Figure 7). The new parameterization had a smaller RMSE (about 5 W/m²) than the two algorithms for VZA ≤ 50° (comparing Figure 6 and Figure 7) but had a slightly larger RMSE (about 1 to 9 W/m²) 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.

3.2. 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. The definition of each case and the error range of input parameters.

Input Parameter	Value of All Parameters	Error of Input Parameter
IWV	Case 1 (BT = 260 K, IWV = 0.5 cm), δTs,a = 0 K	[−90%, 400%]
	Case 2 (BT = 300 K, IWV = 0.5 cm), δTs,a = 0 K	[−90%, 400%]
	Case 3 (BT = 300 K, IWV = 5.0 cm), δTs,a = 0 K	[−60%, 60%]
BT	Case 1, Case 2, Case 3, δTs,a = 0 K	[−20 K, 20 K]
δTs,a	Case 1, Case 2, Case 3, δTs,a = −10 K	[−5 K, 5 K]
	Case 1, Case 2, Case 3, δTs,a = 10 K	[−20 K, 20 K]

, 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.

3.2.1. 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² and −13.7 to 9.5 W/m² 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² and −43.5 to 61.7 W/m² 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² 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²) 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² and −74.3 to 75.1 W/m² 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.

3.2.2. 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² 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² for Case 2 and −29.5 to 31.0 W/m² for Case 3.

3.2.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² for Cases 1 and 2, but it reached 17.3 W/m² 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² (or overestimated by 23.1, 40.9, and 25.2 W/m²) 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² 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.

4. Validation Results of Satellite Retrieved DLRs

4.1. Satellite and Surface Data

4.1.1. 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. Input parameters and data sources used in this study.

Parameters	Data Source	Algorithms Requiring the Data
Geolocation information and elevation	MOD03, at 1 km	All algorithms
Cloud mask	MOD35, at 1 km	All algorithms
Viewing zenith angle	MOD03, at 1 km	NewParam, Tang-Li, Wang-Liang
TIR channel radiance	MOD021KM, at 1 km	NewParam, Tang-Li, Wang-Liang
Water vapor content	MOD05, at 1 km for daytime and 5 km for nighttime	NewParam, Zhou-Cess, Gupta2010
Atmospheric profiles	MOD07, at 5 km	Gupta2010
Air and surface temperatures	MOD07, at 5 km	NewParam, Zhou-Cess, Gupta2010

. 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].

4.1.2. NCEP Data and Atmospheric Parameter Derivation

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.

4.1.3. 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. Description of site conditions. Sites used as examples of dry-arid regions are marked with asterisks.

Site Label	Geographic Name	Latitude/Longitude (°)	Elevation(m)	Period	Surface Type	Instrument	DLR Accuracy	Network
YK	Yingke oasis station, China	38.850/100.417	1519	2010	Cropland	Kipp & Zonen CNR1	±10%	WATER
AR	Arou freeze/thaw station, China	38.050/100.450	3033	2010	Grass	Eppley PIR	±5%	WATER
HZZ *	Huazhaizi desert station, China	38.767/100.317	1726	2010	Desert	Kipp & Zonen CNR1	±10%	WATER
SGP-C1	South Great Plains, Central Facility, America	36.600/−97.500	318	2010	Grass	Eppley PIR	±2.5% or ±4 W/m2	ARM
NSA-C1	Barrow, North Slope of Alaska, America	71.300/−156.600	7.6	2010	Tundra	Eppley PIR	±6 W/m2	ARM
TWP-C1	Manus Island, Tropical Western Pacific	−2.100/147.400	4	2010	Grass	Eppley PIR	±6 W/m2	ARM
ASP *	Alice Springs, Australia	−23.798/133.888	547	2010	Grass	Eppley PIR	±10 W/m2	BSRN
DOM	Dome C, Antarctica	−75.100/123.383	3233	2010, [1,2]	Glacier	Kipp & Zonen CG4	±3%	BSRN
DRA *	Desert Rock, America	36.626/−116.018	1007	2010	Desert	Eppley PIR	±10 W/m2	BSRN
SBO *	Sede Boqer, Israel	30.905/34.782	500	2010	Desert	Eppley PIR	±10 W/m2	BSRN
SPO	South Pole, Antarctica	−89.983/−24.799	2800	2010, [1,7,8]	Glacier	Eppley PIR	±10 W/m2	BSRN
TAM *	Tamanrasset, Algeria	22.780/5.510	1385	2010	Desert	Eppley PIR	±10 W/m2	BSRN
FMF	Flagstaff Managed Forest, Arizona	35.1426/−111.7273	2160	2010	Forest	Kipp & Zonen CNR1	±10%	Ameriflux
GLEES	US-GLE Wyoming, USA	41.3644/−106.2394	3190	2010	Forest	Eppley PIR	±5%	Ameriflux

. 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² 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].

4.2. 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. Error statistics for MODIS-derived clear-sky DLR using different algorithms. RMSE, bias, and the square of the correlation coefficients (R²) of all sites and the average values of each type of region are given. RMSE and bias are in units of W/m². The best results are highlighted in bold, and the worst results are marked with asterisks.

Site Label	NewParam			Tang-Li			Wang-Liang			Zhou-Cess			Gupta2010
	RMSE	Bias	R2	RMSE	Bias	R2	RMSE	Bias	R2	RMSE	Bias	R2	RMSE	Bias	R2
Sites in ordinary regions (DEM < 2000 m)
YK	23.5	1.2	0.85	23.3	−1.9	0.83	24.1	−12.2	0.86	35.4 *	−21.6	0.74	21.3	−9.3	0.88
SGP-C1	19.4	8.0	0.91	50.2 *	39.8	0.75	21.6	2.4	0.86	24.1	7.8	0.85	25.6	16.9	0.91
NSA-C1	17.7	−1.8	0.95	24.6 *	13.3	0.93	20.3	10.1	0.94	20.6	10.4	0.94	17.7	−6.3	0.95
TWP-C1	22.3 *	−9.2	0.05	19.2	−8.1	0.01	19.9	6.8	0.00	20.5	−16.1	0.18	19.3	−11.3	0.15
Mean	20.5	1.1	0.95	34.3 *	14.4	0.87	21.5	1.8	0.94	25.9	−2.4	0.91	21.7	0.1	0.94
Sites in dry-arid regions
HZZ	24.3	−6.5	0.90	26.7	11.8	0.84	33.2	−9.0	0.82	33.2	−25.1	0.84	26.3	−20.5	0.91
ASP	19.5	0.6	0.82	49.9 *	40.7	0.66	24.0	12.6	0.77	18.0	−6.9	0.85	18.4	−9.6	0.86
DRA	26.8	7.4	0.79	48.6 *	35.7	0.67	30.2	−21.5	0.79	23.7	−2.9	0.73	21.7	−1.9	0.81
SBO	23.9	0.8	0.72	48.6 *	36.4	0.59	27.7	−3.3	0.63	26.1	−13.4	0.66	26.1	−13.8	0.71
TAM	21.9	6.7	0.77	34.4	16.9	0.46	39.2 *	−17.5	0.27	18.1	−3.6	0.81	18.4	3.4	0.84
Mean	23.3	3.3	0.80	43.6 *	29.3	0.61	31.9	−9.0	0.60	22.8	−8.0	0.77	22.3	−1.6	0.81
Sites in high-altitude regions (DEM ≥ 2000 m)
AR	23.7	9.6	0.82	23.7	8.6	0.71	21.4	−4.3	0.74	27.1	−10.4	0.67	27.2 *	−18.8	0.83
DOM	7.8	1.4	0.56	23.3	18.0	0.31	13.4	2.9	0.30	34.8 *	34.0	0.57	12.4	9.8	0.58
SPO	15.5	−10.6	0.83	17.3	0.5	0.77	25.3	−17.3	0.94	29.6 *	27.0	0.78	24.3	−12.6	0.71
FMF	22.5	−8.2	0.80	26.3	−0.1	0.66	19.4	5.3	0.83	23.8	3.0	0.74	17.9	5.3	0.87
GLEES	35.7	6.3	0.62	33.3	12.6	0.65	29.4	10.9	0.72	45.4	14.2	0.39	34.6	8.7	0.61
Mean	22.1	−0.6	0.94	25.8	9.5	0.92	23.1	6.1	0.94	33.5 *	15.3	0.89	23.1	0.1	0.93

, 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 Figure 9, Figure 10 and Figure 11.

4.2.1. 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². 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²), which were at NSA-C1, and had some positive and negative errors at the higher DLRs (~400 W/m²), 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.8 W/m² 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² at both sites.

DLRs accuracies ranged from 7.8 to 35.7 W/m² 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² 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.

4.2.2. 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² smaller than those from Wang-Liang and Gupta2010 algorithms and was 13.8 and 5.4 W/m² 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² (at YK and SGP-C1), while the DLRs ≥400 W/m² (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.0 W/m²), and much smaller than those from the Tang-Li and Wang-Liang algorithms (~8.5–20.3 W/m²). 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²) and much better than the Zhou-Cess algorithm (with mean RMSE decreased 11.4 W/m²). For the DLR within 200–400 W/m² 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.

4.3. 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.

4.3.1. 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.

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.

4.3.2. 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. Same as Table 7 but 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.

Event Label	NewParam			Zhou-Cess			Gupta2010
	RMSE	Bias	R2	RMSE	Bias	R2	RMSE	Bias	R2
AR	25.2(1.5)	7.0	0.78	21.6(−5.5)	−3.2	0.84	24.8(−2.4)	−15.6	0.87
DOM	7.8(0)	3.3	0.63	35.1(+0.3)	34.7	0.73	14.6(+2.2)	13.3	0.74
SPO	12.4(−3.1)	−8.4	0.88	29.4(−0.2)	28.5	0.96	5.3(−19)	0.5	0.96
FMF	26.8(4.3)	−3.4	0.77	20.2(−3.6)	−10.4	0.87	21.8(+3.9)	−14.5	0.88
GLEES	27.9(−7.8)	−11.4	0.74	26.3(−19.1)	2.3	0.73	27.9(−6.7)	−8.4	0.71
Mean	21.3(−0.8)	−2.2	0.94	27.8(−5.7)	10.8	0.96	20.5(−2.6)	−3.1	0.96

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² at AR and GLEES for Zhou-Cess and decreased by 2.4, 19.0 and 6.7 W/m² 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²) 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.

5. 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², 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² 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². 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² 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²).

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.