Estimation of Downwelling Surface Longwave Radiation under Heavy Dust Aerosol Sky

The variation of aerosols, especially dust aerosol, in time and space plays an important role in climate forcing studies. Aerosols can effectively reduce land surface longwave emission and re-emit energy at a colder temperature, which makes it difficult to estimate downwelling surface longwave radiation (DSLR) with satellite data. Using the latest atmospheric radiative transfer code (MODTRAN 5.0), we have simulated the outgoing longwave radiation (OLR) and DSLR under different land surface types and atmospheric profile conditions. The results show that dust aerosol has an obvious “warming” effect to longwave radiation compared with other aerosols; that aerosol longwave radiative forcing (ALRF) increased with the increasing of aerosol optical depth (AOD); and that the atmospheric water vapor content (WVC) is critical to the understanding of ALRF. A method is proposed to improve the accuracy of DSLR estimation from satellite data for the skies under heavy dust aerosols. The AOD and atmospheric WVC under cloud-free conditions with a relatively simple satellite-based radiation model yielding the high accurate DSLR under heavy dust aerosol are used explicitly as model input to reduce the effects of dust aerosol on the estimation of DSLR. Validations of the proposed model with satellites data and field measurements show that it can estimate the DSLR accurately under heavy dust aerosol skies. The root mean square errors (RMSEs) are 20.4 W/m2 and 24.2 W/m2 for Terra and Aqua satellites, respectively, at the Yingke site, and the biases are 2.7 W/m2 and 9.6 W/m2, respectively. For the Arvaikheer site, the RMSEs are 23.2 W/m2 and 19.8 W/m2 for Terra and Aqua, respectively, and the biases are 7.8 W/m2 and 10.5 W/m2, respectively. The proposed method is especially applicable to acquire relatively high accurate DSLR under heavy dust aerosol using MODIS data with available WVC and AOD data.


Introduction
As an important component of global radiation balance, downwelling surface longwave radiation (DSLR) is of major interest for the study of land surface processes and land-atmosphere exchange due to its applicability in hydrology, meteorology and climatology [1,2].It generally has a heat-preservation effect on the earth and is a key variable to clarify the temperature dependency of worldwide ablation rates [3].DSLR is the result of atmospheric scattering, absorption and emission in the entire vertical column, and it can be accurately calculated with complex radiative transfer models for which the concentrations of atmospheric constituents and the vertical distribution of temperature, water vapor, aerosol and clouds as the input parameters [4][5][6][7].It is therefore necessary to evaluate accurately these parameters if one wants to reduce uncertainties in computing their radiative effects [8,9].However, acquisition of these parameters is difficult and they are varying rapidly, at a daily or even hourly scale.Synoptic observations, which have close relationships with DSLR, are used to develop the estimating model.Typically, ground level temperature and humidity are used as parameterization factors for calculation of the DSLR under clear skies.Most of the parameterization schemes are based on empirical relationships between the DSLR and atmospheric temperature and humidity, or both, although some are based on radiative transfer theory [10][11][12][13][14][15].These empirical schemes can perform well under clear sky conditions, especially in climatic conditions similar to those in which they were derived.However, under heavy dust aerosol conditions, their performances degrade severely unless aerosol-dependent corrections are made.
Dust aerosol is considered the main source of troposphere aerosol loading and is a source of great uncertainty in climate change modeling [16].Dust aerosols vary in time in space, and their concentrations alter the microphysical properties of clouds in a process known as the "first indirect" effect [17,18].Aerosol radiative effects in the shortwave spectral range have been widely investigated [19][20][21][22][23][24][25][26].In a recent study, Valenzuela et al. [27] obtained the daily averages of the aerosol shortwave radiative forcing (ASRF: which is the shortwave radiation difference between those with aerosols and the corresponding expected clear-sky conditions) at the surface and at the top of the atmosphere (TOA) during desert dust events from 2005 to 2010 in Granada, Spain.The ASRF mean monthly values ranged from −13 W/m 2 to −34 W/m 2 at the surface and −4 W/m 2 to −13 W/m 2 at the TOA.The negative values of approximately 9 W/m 2 to 21 W/m 2 between the surface and the TOA indicate that desert dust induces a reduction of the net shortwave radiation.This suggests that there is a non-negligible atmospheric absorption of solar radiation that induces a significant warming of the atmosphere.Few studies have investigated the longwave radiative effects associated with desert dust intrusions [28][29][30][31][32].These works note that, on the daily timescale, the longwave effect may compensate for a large fraction (up to 50%) of the substantial decrease of shortwave irradiance recorded at the surface during desert dust episodes.Thus, most studies related to this relevant issue have been conducted over deserts with heavy aerosol load conditions [33][34][35][36][37][38][39].
The frequent occurrence of dust aerosols can exert an influence on the surface longwave radiation fluxes that is comparable to a warming effect caused by the established greenhouse gases and implies that the observed longwave enhancement is climatologically significant.Dust aerosols are effective in reflecting solar shortwave irradiance back to space, thereby "cooling" the global surface [40,41].With equivalent particle diameters on the order of several micrometers [42], they can also effectively reduce land surface longwave emission and re-emit energy at a colder temperature, thereby "warming" the earth [43][44][45].However, numerical simulation of dust aerosol longwave forcing is difficult due to its high spatio-temporal variation and secondary effects, for which it can act as cloud condensation nuclei in the process of emission, transportation and deposition [46].Estimation of the atmospheric downward longwave irradiance under a heavy dust aerosol sky remains a challenging task.
Prior to MODIS, satellite data were limited to reflectance measurements in one or two channels.There was no attempt to retrieve aerosol content on a global scale, though algorithms were developed for use over vegetation or other dark land cover types.The blue channel of MODIS offers the possibility to extend the derivation of optical thickness over land to additional surfaces.Well-calibrated measurements from the MODIS sensor are now available for detecting the optical properties of aerosols, enabling examination of the role of aerosols on longwave radiative forcing [47].As can be seen in Table 1, compared with previous imagers, MODIS has 36 channels with improved spectral, spatial and radiometric resolutions [48] and has been used to identify and retrieve aerosol properties over the global surface [49][50][51].Our study is focused on the retrieval of DSLR, the thermal channels which are set for the detection of longwave radiation were used to develop the new model.The method discussed here employs a relatively simple and computationally efficient formulation that utilizes readily available satellite-derived measurements without the need of a radiative transfer model or profiles of atmospheric states to estimate the DSLR under heavy dust aerosol sky from MODIS data.The sensitivity of the new model was analyzed, and a validation of the new method was conducted using the Watershed Allied Telemetry Experimental Research dataset and Arvaikheer site data [52].The general idea is firstly attempted to quantify the longwave radiative forcing of dust aerosol using its optical parameters and atmospheric parameters, subsequently, by improving the clear-sky DSLR (DSLR cl ) estimation model in [52][53][54][55][56], we generalized a formulation adjust to estimate the DSLR under dust aerosol sky conditions.Section 2 describes the data and test sites used in this study.Section 3 presents the method to derive the DSLR from the TOA radiance and other satellite-derived parameters.The results and sensitivity, as well as other issues related to the proposed model, are analyzed in Section 4. Section 5 describes validations of the proposed method applied to MODIS satellite data.Finally, conclusions are given in Section 6.

TIGR Atmospheric Profiles and Simulated Data
To analyze the longwave radiative forcing of dust aerosol and to develop a model which is available for estimating DSLR under dust aerosol, it is crucial to obtain a dataset that includes a wide range of coincident measurements of DSLR, AOD, surface temperature, and emissivities as well as the channel TOA radiance and outgoing longwave radiation (OLR) simultaneously with the pixel scales.However, it is not easy to synchronously measure such data in field experiments.Therefore, MODTRAN (Moderate Resolution Atmospheric Transmittance and Radiance Code) was used to simulate the DSLR at ground level under different atmospheric profile and land surface conditions, and the corresponding MODIS channel radiances at the TOA were generated [53].
To simulate the TOA radiances under dust aerosol with worldwide atmospheric situations and land surface types, the Thermodynamic Initial Guess Retrieval (TIGR) database [54], which are constructed by the Laboratoire de Meteorologie Dynamique (LMD) and represent a worldwide set of atmospheric conditions from polar to tropical atmosphere, are used.Because the objective of this study focuses on DSLR retrieval in the presence of dust aerosol in cloud-free skies, only the atmospheric profiles under clear skies are considered.Most of the clear sky atmospheric profiles are in dry conditions, to reduce the effect of the mal-distribution and make it more representative, 60 of the total 952 atmospheric profiles are chosen from different latitudinal band to calculate the longwave radiation forcing of aerosol, each of which includes forty levels of pressure, atmosphere temperature, water vapor and ozone.The atmosphere water vapor content (WVC) are from 0.05 g/cm 2 to 6.15 g/cm 2 , and the bottom atmospheric temperatures are from 235 K to 314 K, which can represent the global common atmospheric conditions.Surface emissivities in our MODTRAN simulations were obtained from the ASTER spectral library [55].Since there is no spectral emissivity available beyond the 14 µm wavelength and considering the strong absorption of the atmosphere at spectral wavelengths greater than 14 µm, the surface emissivity used in MODTRAN simulation beyond this wavelength is assumed to be unity.
In MODTRAN, CARD 2 was used to describe the aerosol parameters.IHAZE, ISEASN, IVULCN, and VIS select the altitude and seasonal-dependent aerosol profiles and aerosol extinction coefficients.IHAZE specifies the aerosol model used for the boundary-layer (0 to 2 km) and a default.ISEASN selects the seasonal dependence of the profiles for both the tropospheric (2 to 10 km) and stratospheric (10 to 30 km) aerosols.IVULCN is used to select the profile and extinction type for the stratospheric aerosols and to determine transition profiles through the stratosphere to 100 km.VIS specifies the surface meteorological range overriding the default value associated with the boundary layer chosen by IHAZE.It allows one to define the 550 nm aerosol Rayleigh vertical optical depth when the value is negative.The optical properties (normalized extinction coefficient, single scattering albedo and asymmetry parameter) of dust aerosol are acquired from the Optical Properties of Aerosol and Cloud (OPAC) software package [56], which provides the optical properties of different aerosol and cloud components for wavelengths ranging from visible to TIR and is widely used for general circulation models (GCMs) and climate research.The aerosol extinction coefficient σ ext (λ), for ambient particles with number size distribution n(r), is given by: where Q e (λ, r) is the extinction efficiency of dust aerosol at the wavelength of λ, and r is the radius.We normalized the aerosol extinction coefficient divided by the extinction at the wavelength of 550 nm. Figure 1 shows the normalized extinction coefficient, single scattering albedo and asymmetry parameter of dust aerosol from 6 to 15 µm wavelengths.The lowest normalized extinction coefficient and single scattering albedo are in channel 29, which suggests that there is high atmospheric transmittance and weak scattering in this channel wavelength.In channels 31 and 32, the normalized extinction coefficient and asymmetrical parameters are almost identical.From channels 30 to 36, the normalized extinction coefficient and single scattering albedo decrease but the asymmetry parameters basically remain unchanged.Because dust aerosol can significantly reduce the thermal infrared transmittance of the atmospheric window channel [57], our research scope was focused on aerosol optical depth (AOD) less than 1 in 0.55 µm.For comparison, other common types of aerosols (rural aerosol, urban aerosol, maritime aerosol) were used in our simulation run.To increase the representativeness, reasonable variations of land surface temperature (LST) were set according to the ground level atmospheric temperature (T 0 ), which varied from T 0 − 5 K to T 0 + 15 K in 5 K steps.For the angular dependence of the TOA radiance, six different viewing zenith angles (VZAs) were considered (secant (VZA) = 1.0, 1.2, 1.4, 1.6, 1.8, and 2.0)).After integrating the simulated TOA radiances with the channel response functions of MODIS TIR channels, the MODIS channel radiances were calculated.
representativeness, reasonable variations of land surface temperature (LST) were set according to the ground level atmospheric temperature (T0), which varied from T0 − 5 K to T0 + 15 K in 5 K steps.For the angular dependence of the TOA radiance, six different viewing zenith angles (VZAs) were considered (secant (VZA) = 1.0, 1.2, 1.4, 1.6, 1.8, and 2.0)).After integrating the simulated TOA radiances with the channel response functions of MODIS TIR channels, the MODIS channel radiances were calculated.

Groundbased Radiation Measurements and MODIS Satellite Data
Two test sites were used to determine the quality of the improved model.The Yingke (100 • 24'E, 38 • 51'N) site is in the middle Heihe River Basin in northwest China at an altitude of 1520 m and is surrounded by a large and homogeneous corn field.As a relatively ideal proving ground, it is usually covered by dust aerosol and is selected as a test site in our study.The Arvaikheer (46 • 14'N, 102 • 49'E) site is in Mongolia at an altitude of 1940 m and is covered by grass and usually affected by dust aerosol from the Gobi Desert.The DSLR in the two radiations were measured using pyrgeometers; the sampling interval is five minutes from January 2007 to December 2009 for the Yingke site and from January 2003 to December 2003 for the Arvaikheer site; and the instantaneous value was derived by bilinear interpolation of the data with measured time of fifteen minutes before and after the satellite overpass.Seven MODIS products were are also used: (1) MODIS calibrated radiances (MOD021KM/MYD021KM), which contains geo-located and calibrated earth view observations, and creates tiled MODIS images at 1-km resolution data.The observed thermal TOA radiances are used in our study.(2) MODIS geolocation data (MOD03/MYD03), which contains geodetic latitude and longitude, surface height above geoid, solar zenith and azimuth angles, satellite viewing zenith and azimuth angles, and a land/sea mask for each 1 km sample.The latitude and longitude data are used to perform geometric corrections for other products and the VZAs are used to determine the coefficients of the proposed method.(3) MODIS aerosol datasets (MOD04_L2/MYD04_L2), which include aerosol optical thickness at three visible wavelengths, a measure of the fraction of aerosol optical thickness attributed to the fine mode, and several derived parameters including reflected spectral solar flux at the top of the atmosphere.(4) MODIS water vapor datasets (MOD05_L2/MYD05_L2), which are consist of column water-vapor amounts, and are generated at the 1-km spatial resolution of the MODIS instrument using the near-infrared algorithm during the day, and at 5 × 5 1-km pixel resolution both day and night using the infra-red algorithm when at least nine FOVs are cloud free.(5) MODIS land surface temperature and emissivity (MOD11_L2/MYD11_L2), which provides per-pixel LST and emissivity values at 1 km resolutions.( 6) MODIS cloud detection product (MOD035_L2/MYD035_L2), which provided all masking information, including most recently an internal cloud mask based on spatial variability to identify low clouds and the reflectance in the 1.38-m channel to identify high clouds.
These products are available in hierarchical data format (HDF); MODIS products are provided by NASA's Goddard Space Flight Center (GSFC) and the Atmosphere Archive and Distribution System (LAADS) [58], and the GLASS BBE data are available from the Global Land Cover Facility of the University of Maryland [59].

Methodology
The paper focuses on the analysis of the influence of dust aerosol on the DSLR and the estimation of the DSLR under dust aerosol conditions.The contribution of dust aerosol to the longwave spectral range can be defined as aerosol longwave radiative forcing (ALRF).The DSLR is formulated as a simple, physically based bulk expression that accounts for dust aerosol variability.At a given pixel, the formulation describing the DSLR under dust aerosol conditions may be written as the DSLR under clear sky added to the downwelling radiation forcing generated by the dust aerosol.
where DSLR aero is the DSLR under dust aerosol conditions, DSLR aero− f ree is the corresponding expected clear-sky conditions and ALRF is the downwelling radiation forcing of dust aerosol.
According to Tang and Li [60], the channel radiance received by the satellite thermal infrared sensor at the TOA is approximately the aggregation of the radiance from the earth's surface and each level of the atmosphere.Therefore, an algebraic expression dτ i (p)/dInp is introduced as the channel weighting function.It can depict the magnitude of radiance contribution from different atmospheric vertical layers.In the expression, dτ i (p) and dInp are the channel transmittance difference and the logarithmic atmosphere pressure difference between two adjacent atmosphere levels, respectively.For MODIS, each weighting function of its TIR channel implies that it can detect information at different altitudes in a vertical atmospheric column.The max value means that this atmosphere level has the maximum magnitude of radiance contributing to the sensor channel radiance.Based on the algebraic expression result of each channel, and assuming there is a linear relationship between the DSLR and MODIS channel radiances, for clear sky with no dust aerosols and Lambertian surface conditions, Tang and Li [60] developed a new model to retrieve the DSLR directly from the TOA radiances measured by MODIS TIR channels with the formula: where a 0 to a 6 are the regression coefficients related to VZA and the altitude, and L i (θ)  , it suggests that the retrieved value can be used to help narrow the uncertainties associated with DSLR researches.Furthermore, unlike most of the statistical algorithms that need ground-level atmosphere temperature as an input, it estimates the DSLR directly from the TOA radiance and thus was used in our research.By quantifying the longwave radiative forcing of dust aerosol and improving the clear-sky DSLR estimation models, we generalized a formulation adjust to estimate the DSLR under dust aerosol sky conditions.

Results
According to the definition and the simulation data set, the ALRF under different types and contents of aerosols are shown below.In addition, the impacts of aerosol on each MODIS channels used by the DSLR model are described.
Figures 2 and 3 show the ALRF and OLR influenced by different types and optical depths of aerosols.The symbols indicate the variations and the half-length of vertical bars represent the Standard Deviations (SDs) of the two parameters.Notably, dust aerosol has positive effects on the ALRF but negative effects on the OLR, and both effects are enhanced as the AOD increased; the Standard Deviations of the ALRF and OLR also show similar changes.For maritime aerosol, the ALRF is relatively steady as the AOD increase, the maximum is approximately 3.85 W/m 2 when the AOD is 1.0, with a Standard Deviations of 1.1 W/m 2 .For rural and urban aerosols, the ALRF is 4.85 W/m 2 and 6.22 W/m 2 , respectively, with Standard Deviations of 1.32 W/m 2 and 1.85 W/m 2 , respectively.For dust aerosol, the ALRF increases from 2.09 W/m 2 to 29.8 W/m 2 when the AOD increases from 0.1 to 1.0, with Standard Deviations from 0.66 W/m 2 to 8.3 W/m 2 , respectively.However, the OLR only changes slightly under the four aerosol types and the maximum is -1.26W/m 2 when the AOD = 1.0 for dust aerosol, which means that dust aerosol has little influence on the thermal TOA radiance.In addition, this view has been confirmed in Figure 4, the average rate variation of TOA radiance of each band is different.The TOA radiances of bands 29 and 31 are more sensitive to dust aerosol.The maximum rate variations reach 3.5% and -3.5% for bands 29 and 31, respectively, for AOD = 1.0.Moreover, from these figures, we can see that the ALRF of the four aerosol types can be expressed as linear functions of the AOD with coefficients of determination (R 2 ) of 0.998, 0.998, 0.997 and 0.997.After analyzing the relationships between the ALRF and AOD for other situations, the following linear function can be concluded: where k and b are the slope and offset, respectively, which vary with the condition of the atmospheric profiles and are determined by the atmospheric WVC, atmospheric temperature, CO 2 , O 3 and other constituents.
Remote Sens. 2017, 9, 207 where k and b are the slope and offset, respectively, which vary with the condition of the atmospheric profiles and are determined by the atmospheric WVC, atmospheric temperature, CO2, O3 and other constituents.Figure 5 shows the scattering diagram of the offset b in Equation ( 5) as function of the WVC.The offset had little variation, from −1.45 to 0.3 for dust aerosol, and can be expressed as constants for aerosol types.Therefore, the ALRF is determined by its slope k, which has a close relationship with the compositions of the atmosphere, such as CO2, WVC, trace gases, etc. Wang and Liang [61] found that the daily DSLR increased at an average rate of 2.2 W/m 2 per decade from 1973 to 2008, and that the rising trend has a close relationship with the increasing CO2 concentration.We researched the influence of CO2 on the DSLR and the results show that the DSLR has increased 2.0 W/m 2 as the CO2 concentration doubled, consistent with the study by Wilson and Gea-Banacloche [62].The content of CO2 and other trace gases are steady and relatively small and have a far smaller impact on the DSLR than the WVC.    Figure 5 shows the scattering diagram of the offset b in Equation ( 5) as function of the WVC.The offset had little variation, from −1.45 to 0.3 for dust aerosol, and can be expressed as constants for aerosol types.Therefore, the ALRF is determined by its slope k, which has a close relationship with the compositions of the atmosphere, such as CO 2 , WVC, trace gases, etc. Wang and Liang [61] found that the daily DSLR increased at an average rate of 2.2 W/m 2 per decade from 1973 to 2008, and that the rising trend has a close relationship with the increasing CO 2 concentration.We researched the influence of CO 2 on the DSLR and the results show that the DSLR has increased 2.0 W/m 2 as the CO 2 concentration doubled, consistent with the study by Wilson and Gea-Banacloche [62].The content of CO 2 and other trace gases are steady and relatively small and have a far smaller impact on the DSLR than the WVC.Thus, the WVC is used to parameterize k.For dust aerosol, Figure 6 shows the scattering diagrams of the slope k in Equation ( 5) as a function of the atmospheric WVC.If we take the slope k as a linear function of atmospheric WVC, the RMSE is 3.81 W/m 2 and the determination coefficient (R 2 ) is 0.79.The ALRF can be roughly calculated with WVC and AOD, so Equation ( 2) can be expressed as:  Thus, the WVC is used to parameterize k.For dust aerosol, Figure 6 shows the scattering diagrams of the slope k in Equation ( 5) as a function of the atmospheric WVC.If we take the slope k as a linear function of atmospheric WVC, the RMSE is 3.81 W/m 2 and the determination coefficient (R 2 ) is 0.79.The ALRF can be roughly calculated with WVC and AOD, so Equation ( 2) can be expressed as: where DSLR aero is the DSLR under aerosol sky; DSLR aero− f ree is the DSLR in cloud-free sky when no dust aerosol is present and can be estimated with the DSLR cl model when there are coefficient adjustment; WVC is the atmospheric water vapor content; AOD is the aerosol optical depth; and p, q and b are regression coefficients.Notably, the coefficients in the DSLR clear model may have slight variations because of the slight variations in the OLR between heavy dust aerosol and clear sky.It is well known that the TOA radiance has a close relationship with surface emissivity, which may influence DSLR retrieval.Here, we discussed the accuracy of the improved model on DSLR retrieval for different surface types.Figure 7   It is well known that the TOA radiance has a close relationship with surface emissivity, which may influence DSLR retrieval.Here, we discussed the accuracy of the improved model on DSLR retrieval for different surface types.Figure 7 shows the Root Mean Square Error (RMSE) between the actual DSLR and the model-estimated values for different surface types from VZA = 0 • to VZA = 60 • .When VZA = 0 • , the RMSEs of the proposed model vary from 5.75 W/m 2 to 8.68 W/m 2 for different surface types.For VZA = 60 • , the RMSEs vary from 10.69 W/m 2 to 12.73 W/m 2 .The RMSE of the improved model generally has a close relationship with the surface type and increases as the VZA increases, and the model is VZA dependent.Figure 8 is the comparisons of the actual DSLR with estimations using the proposed model with the training dataset and VZA = 0 • , surface altitude z = 0 km and AOD less than 1.0.The RMSE is 9.16 W/m 2 for all types.For different surface types when AOD = 0.5, as can be seen in Figure 9, the RMSE is from 7.6 W/m 2 to 8.27 W/m 2 .It is well known that the TOA radiance has a close relationship with surface emissivity, which may influence DSLR retrieval.Here, we discussed the accuracy of the improved model on DSLR retrieval for different surface types.Figure 7 shows the Root Mean Square Error (RMSE) between the actual DSLR and the model-estimated values for different surface types from VZA = 0° to VZA = 60°.When VZA = 0°, the RMSEs of the proposed model vary from 5.75 W/m 2 to 8.68 W/m 2 for different surface types.For VZA = 60°, the RMSEs vary from 10.69 W/m 2 to 12.73 W/m 2 .The RMSE of the improved model generally has a close relationship with the surface type and increases as the VZA increases, and the model is VZA dependent.Figure 8 is the comparisons of the actual DSLR with estimations using the proposed model with the training dataset and VZA = 0°, surface altitude z = 0 km and AOD less than 1.0.The RMSE is 9.16 W/m 2 for all types.For different surface types when AOD = 0.5, as can be seen in Figure 9, the RMSE is from 7.6 W/m 2 to 8.27 W/m 2 .

Sensitivity Analysis
To make it convenient to estimate the of the proposed method can be expressed as:

Sensitivity Analysis
To make it convenient to estimate the DSLR aero , we discussed the sensitivity of the proposed model.The performance of the proposed DSLR retrieval model depends on the uncertainty of the input variables, mainly for the MODIS instrument detector noise and calibration error, uncertainties in the WVC and AOD, and in the accuracy of the model itself.Several sources of error must be considered to analyze the sensitivity of the proposed model.The total error e(∆DSLR) of the proposed method can be expressed as: with in which ∆(AOD) and ∆(WVC) are the uncertainty of the WVC and AOD, respectively, and δ(L DSLR ), δ(AOD) and δ(WVC) are the errors induced by the uncertainty of the corresponding parameters.According to Remer et al. [63], on the global scale, the uncertainty of the AOD for land is: 11) δ(AOD) can be calculated using Equation (9).Considering that the MODIS TOA brightness temperature errors are 0.25 K for channels 28, 33, 34 and 35, 0.05 K for channels 29, 31 and 32, and 0.35 K for channel 36, respectively, δ(L DSLR ) can be acquired with Equation (8).Here, we adopted an uncertainty of the WVC for MODIS product of 0.5 g/cm 2 [64], so the corresponding DSLR AOD errors for VZA = 0 • and z = 0 km are as follows.
In Figure 10, it is clear that the δ(L DSLR ), δ(AOD) and δ(WVC) errors caused by the uncertainty of thermal channel radiance, AOD and WVC increased with the increasing of AOD, but the accuracy of the model itself increases, which means that the proposed model is more suitable for heavy dust aerosol.Because the δ(AOD) error has a large proportion and rapid growth, the total error of the proposed algorithm therefore decreased slightly and then increased with the AOD increases, the minimum RMSE corresponds to AOD = 0.4 approximately.It means that the more accurate AOD should be used when estimating the DSLR under dust aerosol with the proposed method.Because the errors from different variables can be neutralized, the actual values are usually smaller than the corresponding values.The total error e(∆DSLR) has a similar tendency for other VZAs and the max value is 16.7 W/m 2 for VZA = 60 • and surface altitude z = 0 km when AOD = 1.0.The improved model has good performance in estimating the DSLR under heavy dust aerosol.
uncertainty of the WVC for MODIS product of 0.5 g/cm 2 [64], so the corresponding AOD DSLR errors for VZA = 0° and z = 0 km are as follows.
In Figure 10, it is clear that the errors caused by the uncertainty of thermal channel radiance, AOD and WVC increased with the increasing of AOD, but the accuracy of the model itself increases, which means that the proposed model is more suitable for heavy dust aerosol.Because the error has a large proportion and rapid growth, the total error of the proposed algorithm therefore decreased slightly and then increased with the AOD increases, the minimum RMSE corresponds to AOD = 0.4 approximately.It means that the more accurate AOD should be used when estimating the DSLR under dust aerosol with the proposed method.Because the errors from different variables can be neutralized, the actual values are usually smaller than the corresponding values.The total error

Validations with In-Situ Measurements
It is important to note that the MODIS aerosol product monitors the AOD over part of the continents.Daily Level 2 data are produced at the spatial resolution of a 10 km×10 km (at nadir)-pixel array.MOD04_L2 and MYD04_L2 contain data collected from Terra and Aqua, respectively.The Dust aerosols are identified with MOD04/MYD04 product and the emissivity method [65].The DSLR with the measured time closest to the MODIS Terra/Aqua overpass time was extracted to represent

Validations with In-Situ Measurements
It is important to note that the MODIS aerosol product monitors the AOD over part of the continents.Daily Level 2 data are produced at the spatial resolution of a 10 km×10 km (at nadir)-pixel array.MOD04_L2 and MYD04_L2 contain data collected from Terra and Aqua, respectively.The Dust aerosols are identified with MOD04/MYD04 product and the emissivity method [65].The DSLR with the measured time closest to the MODIS Terra/Aqua overpass time was extracted to represent the equivalent DSLR for the pixel scale.As it was derived by bilinear interpolation using data with measured time adjacent to the satellite overpass, and DSLR has small changes in short time serials, the error caused by the time inconsistency can be ignored.To eliminate the influence of cloud contamination, only the pixels identified as cloud-free were selected.
The objective of the present work is to estimate the instantaneous DSLR under dust aerosol using the proposed method.Only the pixels identified as clear sky were selected based on the quality control data in the MODIS LST MOD11_L2/MYD11_L2 products, from which the dust aerosol sky data were chosen in this study.The ranges of AOD over the YingKe site are from 0.065 to 0.92 and over the Arvaikheer site are from 0.105 to 0.97, respectively.The ranges of WVC over the two sites are from 0.096 to 4.08 g/cm 2 and from 0.16 to 2.34 g/cm 2 , respectively.Figure 11 shows the scattering diagram of the groundbased measured DSLR and the DSLR extracted from MODIS TOA radiance.The solid circles represent the DSLR estimated using the proposed model and the hollow circles show that estimated using the DSLR cl model.The scattering diagrams of the AOD, as a function of the DSLR at the time of MODIS overpass for the two sites, are also shown in Figure 11.To evaluate the performance of the new method under different AOD and WVC conditions sufficiently, the retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLR cl methods (hollow circles) and those groundbased measurements versus AOD at two test sites for dust aerosol conditions at the moment of MODIS overpass are shown in Figure 12.The retrieval errors versus WVC are shown in Figure 13.
sites, are also shown in Figure 11.To evaluate the performance of the new method under different AOD and WVC conditions sufficiently, the retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLRcl methods (hollow circles) and those groundbased measurements versus AOD at two test sites for dust aerosol conditions at the moment of MODIS overpass are shown in Figure 12.The retrieval errors versus WVC are shown in Figure 13.The performances of DSLRcl model degrade severely compared to previous studies [60,66], and the RMSEs can even reach 40 W/m 2 in Yingke site, which is the reason that ALRF and coefficient adjustments are not considered in the model.Dust aerosol can not only generate ALRF but also affects the TOA radiance observed by each bands, especially for bands 29 and 31, the maximum rate variations could reach ±3.5%.Consequently, it should be taken into account the influence of dust aerosol on developing the DLSR estimation model.As we can see in Figure 12, the proposed method performs better under both small and heavy aerosol conditions for the aerosol-dependent corrections are made.Considering that there are variations of TOA radiance when dust aerosol happens and coefficient adjustment has been made in the proposed method, it also can acquire relatively high accuracy DSLR in some of the larger WVC conditions.However, it is obvious that that the RMSEs of the proposed model at the two sites are higher than those in the simulated data.Considering that the ground data measured at the two sites use a point-based scale, whereas the retrieved value from  In addition to these error factors, the performance of the proposed method is also affected by the quality of the input variable of TOA radiance and uncertainties in the AOD and WVC.The vertical distribution of dust aerosol which has different scattering and absorption in different layers has not been considered in the validation yet.A variable height-dependent distribution can explain that most At the Yingke site, most of the dust aerosols occur in spring and autumn, with 118 days of corresponding data for Terra and 89 days for Aqua.The RMSEs of the proposed model are 20.4W/m 2 and 24.2 W/m 2 for Terra and Aqua, respectively, and the biases are 2.7 W/m 2 and 9.6 W/m 2 , respectively.The retrieval errors are from −40 to 50.8 W/m 2 and from −50.4 to 59.8 W/m 2 for Terra and Aqua, respectively.For the Arvaikheer site, there are 39 days of corresponding data for Terra and 34 days for Aqua, and most dust aerosols occur in spring and early summer or autumn.The RMSEs of the proposed model are 23.2W/m 2 and 19.8 W/m 2 for Terra and Aqua, respectively, and the biases are 7.8 W/m 2 and 10.5 W/m 2 , respectively.The retrieval errors are from −51 to 49.2 W/m 2 and from −39.3 to 38.9 W/m 2 for Terra and Aqua, respectively.
The performances of DSLR cl model degrade severely compared to previous studies [60,66], and the RMSEs can even reach 40 W/m 2 in Yingke site, which is the reason that ALRF and coefficient adjustments are not considered in the model.Dust aerosol can not only generate ALRF but also affects the TOA radiance observed by each bands, especially for bands 29 and 31, the maximum rate variations could reach ±3.5%.Consequently, it should be taken into account the influence of dust aerosol on developing the DLSR estimation model.As we can see in Figure 12, the proposed method performs better under both small and heavy aerosol conditions for the aerosol-dependent corrections are made.Considering that there are variations of TOA radiance when dust aerosol happens and coefficient adjustment has been made in the proposed method, it also can acquire relatively high accuracy DSLR in some of the larger WVC conditions.However, it is obvious that that the RMSEs of the proposed model at the two sites are higher than those in the simulated data.Considering that the ground data measured at the two sites use a point-based scale, whereas the retrieved value from MODIS data using our proposed method corresponds to a spatial resolution of 1 km, the heterogeneity has not been considered in the comparison.If the point-based data are assumed to represent the 1-km spatial scale value for comparison with the value derived from MODIS data, it may lead to large discrepancies if the ground stations are located in a heterogeneous area.Dust particles typically are non-spherical [67], however, if they were dealt as spherical particles, which can also affect the accuracy and the validation of the proposed model.Moreover, we assumed that the DSLR error caused by the time inconsistency can be ignored, but the DSLR may change over the time series, especially when there are winds and thin clouds.
In addition to these error factors, the performance of the proposed method is also affected by the quality of the input variable of TOA radiance and uncertainties in the AOD and WVC.The vertical distribution of dust aerosol which has different scattering and absorption in different layers has not been considered in the validation yet.A variable height-dependent distribution can explain that most of the indetermination in forecasting the aerosol radiative impact derives from differences in the aerosol vertical profile employed in models.Lidar measurements provide powerful tools for detection of atmosphere distribution especially for the description of aerosol characteristics.It is likely unrealistic to have a Lidar at the site soon, so a ceilometer can be a more realistic alternative.It is a pity that we do not have those measurements at the two validation sites yet, as the lack of which can affect the accuracy and the validation of the proposed model.OPAC could be an excellent tool to investigate the sensitivity of radiances on different optical properties, however, there are differences between the dust aerosol at the two test sites and the average conditions of dust aerosol in OPAC.Such as aerosol source and aerosol composition, also the properties of aerosol particles are highly variable in time and space, which can affect the estimated results.
In general, compared with the dust aerosol-free model, the proposed model is particularly applicable to acquire relatively high-precision DSLR under heavy dust aerosol from MODIS data for which qualified WVC and AOD are available.The DSLR estimated using the improved model is higher than that estimated using the dust aerosol-free model, especially in high AOD conditions, which is the reason that the ALRF and correction coefficients are considered in the former.There are little differences between the estimated values with the two methods when the AOD is less than 0.2, in which the dust aerosols have little influence on the estimated DSLR and the ALRF is small in these conditions.

Conclusions
In this paper, a method has been proposed to estimate DSLR from remotely sensed data for the cloud-free skies under heavy dust aerosol conditions.The AOD and atmospheric WVC are used as explicit inputs to reduce the effect of dust aerosol on the estimation of DSLR.Influence of VZA on the method was also discussed using the simulated data.For VZA = 0 • , the RMSEs of the proposed model vary from 5.75 W/m 2 to 8.68 W/m 2 for different surface types.For VZA = 60 • , the RMSEs vary from 10.69 W/m 2 to 12.73 W/m 2 .The RMSE of the proposed model increased as the VZA increased and the model is VZA-dependent.As the effect of surface types on the estimated DSLR was slightly, the proposed method did not account for the influence associated with different land covers.
Some groundbased radiation measurements at the Yingke site, China and the Arvaikheer site, Mongolia, are used to validate the proposed method.The results show that the RMSEs of the retrieved DSLR with the proposed method and the site values are 20.4W/m 2 and 24.2 W/m 2 for Terra and Aqua, respectively, at the Yingke site, and the biases are 2.7 W/m 2 and 9.6 W/m 2 , respectively.For the Arvaikheer site, the RMSEs are 23.2W/m 2 and 19.8 W/m 2 , and the biases are 7.8 W/m 2 and 10.5 W/m 2 , for Terra and Aqua, respectively.The proposed method is particularly applicable to MODIS data to acquire high accurate DSLR under heavy dust aerosol for which qualified WVC and AOD are available.
Because of the limited information of aerosol vertical profiles, height-dependent distribution of aerosol is not considered in our study.It may generate indetermination in forecasting the aerosol radiative and estimating of DSLR.Lidar and ceilometer measurements provide powerful tools for detection of atmosphere distribution especially for the description of aerosol characteristics, which can be a feasible method to overcome the shortage and to radically improve the accuracy of DSLR retrieval from space.It is a topic worthy of further research.

Figure 1 .
Figure 1.The normalized extinction coefficient, single scattering albedo and asymmetry parameter of dust aerosol at 6-15 μm wavelengths.

Figure 1 .
Figure 1.The normalized extinction coefficient, single scattering albedo and asymmetry parameter of dust aerosol at 6-15 µm wavelengths.

Figure 2 .
Figure 2. The aerosol longwave radiative forcing (ALRF) of different aerosol types varies with the aerosol optical depth (AOD).The solid symbols indicate the ALRF and the half-length of the vertical bars represents the Standard Deviations of the ALRF.

Figure 2 .
Figure 2. The aerosol longwave radiative forcing (ALRF) of different aerosol types varies with the aerosol optical depth (AOD).The solid symbols indicate the ALRF and the half-length of the vertical bars represents the Standard Deviations of the ALRF.

Figure 2 .
Figure 2. The aerosol longwave radiative forcing (ALRF) of different aerosol types varies with the aerosol optical depth (AOD).The solid symbols indicate the ALRF and the half-length of the vertical bars represents the Standard Deviations of the ALRF.

Figure 3 .
Figure 3.The variations of outgoing longwave radiation (OLR) at different aerosol conditions as the AOD increase.The symbols indicate the variations of OLR and the half-length of vertical bars represent the Standard Deviations of the OLR.

Figure 3 .
Figure 3.The variations of outgoing longwave radiation (OLR) at different aerosol conditions as the AOD increase.The symbols indicate the variations of OLR and the half-length of vertical bars represent the Standard Deviations of the OLR.Remote Sens. 2017, 9, 207 8 of 19

Figure 4 .
Figure 4. Comparing with the clear sky, the rate variations of top of the atmosphere (TOA) channel radiance at different aerosol conditions as the AOD increase.The symbols indicate the rate variations of channel radiance and the half-length of vertical bars represent the Standard Deviations of the rate variations of TOA channel radiance.

Figure 5 .
Figure 5. Scattering diagram of offset b in Equation (5) as a function of atmosphere water vapor

Figure 4 .
Figure 4. Comparing with the clear sky, the rate variations of top of the atmosphere (TOA) channel radiance at different aerosol conditions as the AOD increase.The symbols indicate the rate variations of channel radiance and the half-length of vertical bars represent the Standard Deviations of the rate variations of TOA channel radiance.

Figure 4 .
Figure 4. Comparing with the clear sky, the rate variations of top of the atmosphere (TOA) channel radiance at different aerosol conditions as the AOD increase.The symbols indicate the rate variations of channel radiance and the half-length of vertical bars represent the Standard Deviations of the rate variations of TOA channel radiance.

Figure 5 .
Figure 5. Scattering diagram of offset b in Equation (5) as a function of atmosphere water vapor content (WVC).
in cloud-free sky when no dust aerosol is present and can be estimated with the DSLRcl model when there are coefficient adjustment; WVC is the atmospheric water vapor content; AOD is the aerosol optical depth; and p, q and b are regression coefficients.Notably, the coefficients in the clear DSLR model may have slight variations because of the slight variations in the OLR between heavy dust aerosol and clear sky.

Figure 5 .
Figure 5. Scattering diagram of offset b in Equation (5) as a function of atmosphere water vapor content (WVC).

Remote Sens. 2017, 9 , 207 9 of 19 Figure 6 .
Figure 6.Scattering diagrams of the slope k in Equation (5) as a function of atmospheric WVC under dust aerosol type.
shows the Root Mean Square Error (RMSE) between the actual DSLR and the model-estimated values for different surface types from VZA = 0° to VZA = 60°.When VZA = 0°, the RMSEs of the proposed model vary from 5.75 W/m 2 to 8.68 W/m 2 for different

Figure 6 .
Figure 6.Scattering diagrams of the slope k in Equation (5) as a function of atmospheric WVC under dust aerosol type.

Figure 6 .
Figure 6.Scattering diagrams of the slope k in Equation (5) as a function of atmospheric WVC under dust aerosol type.

Figure 7 .
Figure 7. Root mean square error (RMSE) between the actual downwelling surface longwave radiation (DSLR) and the DSLR estimated using the proposed models for different surface types from VZA = 0° to VZA = 60° when surface altitude z = 0 km.

Figure 8 .
Figure 8. Comparisons of the actual downwelling surface longwave radiation with estimations using the proposed model with the training dataset and VZA = 0 • , surface altitude z = 0 km and AOD less than 1.0.

Figure 8 .
Figure 8. Comparisons of the actual downwelling surface longwave radiation with estimations using the proposed model with the training dataset and VZA = 0°, surface altitude z = 0 km and AOD less than 1.0.

Figure 9 .
Figure 9. Comparisons of the actual downwelling surface longwave radiation with estimations using the proposed model for different surface types and VZA = 0, surface altitude z = 0 km and AOD = 0.5.
aero DSLR , we discussed the sensitivity of the proposed model.The performance of the proposed DSLR retrieval model depends on the uncertainty of the input variables, mainly for the MODIS instrument detector noise and calibration error, uncertainties in the WVC and AOD, and in the accuracy of the model itself.Several sources of error must be considered to analyze the sensitivity of the proposed model.The total error ) ( DSLR e 

Figure 9 .
Figure 9. Comparisons of the actual downwelling surface longwave radiation with estimations using the proposed model for different surface types and VZA = 0, surface altitude z = 0 km and AOD = 0.5.
tendency for other VZAs and the max value is 16.7 W/m 2 for VZA = 60° and surface altitude z = 0 km when AOD = 1.0.The improved model has good performance in estimating the DSLR under heavy dust aerosol.

Figure 10 .
Figure 10.The aero DSLR error caused by uncertainties of thermal channels, AOD, WVC and the accuracy of the aero free DSLR  model for VZA = 0° and surface altitude z = 0 km.

Figure 10 .
Figure 10.The DSLR aero error caused by uncertainties of thermal channels, AOD, WVC and the accuracy of the DSLR aero− f ree model for VZA = 0 • and surface altitude z = 0 km.

Figure 11 .
Figure 11.Comparisons between the DSLR estimated using the proposed (solid circles) DSLRcl methods (hollow circles) and those groundbased measurements at two test sites for dust aerosol conditions at the moment of MODIS overpass.The scattering diagrams of the AOD (hollow plus signs) as function of the measured DSLR are also shown.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.

Figure 11 . 19 Figure 12 .
Figure 11.Comparisons between the DSLR estimated using the proposed (solid circles) and DSLR cl methods (hollow circles) and those groundbased measurements at two test sites for dust aerosol conditions at the moment of MODIS overpass.The scattering diagrams of the AOD (hollow plus signs) as function of the measured DSLR are also shown.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.Remote Sens. 2017, 9, 207 13 of 19

Figure 12 .
Figure 12.The retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLRcl methods (hollow circles) and those groundbased measurements versus AOD at two test sites for dust aerosol conditions at the moment of MODIS overpass.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.

Figure 12 .
Figure 12.The retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLR cl methods (hollow circles) and those groundbased measurements versus AOD at two test sites for dust aerosol conditions at the moment of MODIS overpass.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.

Figure 13 .
Figure 13.The retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLRcl methods (hollow circles) and those groundbased measurements versus WVC at two test sites for dust aerosol conditions at the moment of MODIS overpass.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.

Figure 13 .
Figure 13.The retrieval errors between DSLR extracted from MODIS data with the proposed (solid circles) and DSLR cl methods (hollow circles) and those groundbased measurements versus WVC at two test sites for dust aerosol conditions at the moment of MODIS overpass.(A) Yingke-terra; (B) Yingke-aqua; (C) Arvaikheer-terra; (D) Arvaikheer-aqua.
is the TOA radiance measured by the MODIS TIR channel.The model has been used to estimate the DSLR over the USA and tested at Bondville, Boulder, Desert Rock, Penn State, and Sioux Falls, which represent different surface types.The RMSEs of the linear model are 23.6 W/m 2 to 39.3 W/m 2 compared with the measured data at the six sites.Considering that the MODIS instrument error is approximately 11 W/m 2