Estimation of Net Surface Shortwave Radiation From Remotely Sensed Data Under Dust Aerosol Conditions

Due to the complexity and heterogeneity mechanism of aerosol scattering and absorption, there is no systematic operational net surface shortwave radiation (NSSR) retrieval algorithm for considering the forces of different aerosols. This study is to estimate NSSR under dust aerosol conditions. A parameterized model is proposed to quantify the net shortwave radiative forcing (ASRF) of aerosol. Considering the particularity of dust aerosol on shortwave radiation and its mechanisms, aerosol optical properties and atmospheric water vapor content (WVC) are used to quantify the ASRF, and the previous Tang’s parameterization model for estimating NSSR is refined to estimate NSSR under dust aerosol conditions. On this basis, a more accurate radiative transfer model MODTRAN 5 is primitively utilized to simulate and fit the complex relationship. Moreover, the coefficients for parameterization are recalculated by an improved algorithm under diverse types of aerosols. Finally, the NSSR measurements in the study areas and the NSSR products from the MODIS satellites are utilized to verify the retrieved NSSR. The root-mean-square error (RMSE) is below 11 W/m2 under various land cover types and the research findings can significantly increase the estimation accuracy of NSSR around 40 W/m2 under dust aerosol conditions.


I. INTRODUCTION
Surface shortwave radiation is an important physical quantity in the simulation of global, regional climate, hydrology and land surface processes, and its high-precision remote sensing retrieval is of great practical significance for global and regional climate change, energy balance, ecological environment, atmospheric cycle exploration and other fields [1]. The net surface shortwave radiation (NSSR) is the main factor in the climate circulation and global change, which is measured significantly in Net Surface Radiation (NSR) [2]. With the development of remote sensing satellite technology, highprecision retrieval under the consecutive surface at a large scale has been realized. Extensive algorithms and models have been put forward for estimating the shortwave radiation component of the surface from remote sensing observations, The associate editor coordinating the review of this manuscript and approving it for publication was Wenming Cao . and considerable research results have been obtained [3]- [5]. However, the precision cannot be guaranteed under the presence of clouds or strong aerosols conditions, especially the dust aerosols, the accuracy is affected by the aerosols and clouds through complicated scattering and absorption mechanisms interacted with shortwave radiation [6], [7]. Thus, considering the complicated atmosphere condition to increase the retrieval accuracy of the algorithm is imperatively required.
In the past decades, numerous methods for accurate estimation of NSSR in diverse satellites have been explored and generated including statistical empirical method, physically based method, hybrid method, and the machine learning method [8]- [11]. The statistical empirical approach means that the empirical or statistical relationship between surface shortwave radiation measurements and remote sensing observations is directly established without considering the complicated atmospheric radiation transmission mechanism between the sun, atmosphere and land according to the estimation process of NSSR [10], [11]. As a matter of fact, a majority of algorithms have been developed based on the assumption of clear and cloudless sky conditions [12]- [15]. However, as the components of the surface-atmosphere system, aerosol scatters and absorbs solar radiation energy, causing changes in electromagnetic radiation energy, weakening the total radiation energy reaching the surface, and changing the energy received by the sensor, which is called a direct radiative forcing effect [16]. Meanwhile, aerosol particles change the balance of radiative transfer by changing the cloud droplets in the cloud layer, which is an indirect radiative forcing effect [17], [18]. Aerosols are classified into natural aerosols and man-made aerosols on the basis of sources. Natural source aerosols come from the bursting of ocean bubbles, soil erosion, pollen from biological spores, as well as smoke and dust from volcanic eruptions. Anthropogenic aerosols are the salt aerosols or soot particles generated from the combustion of fossil fuels and the secondary aerosols containing polluting dust. Generally, dust aerosol, as a type of natural-source aerosol, accounts for about half of the aerosol in the troposphere and plays a significant role in the energy balance of earth surface. Consequently, aerosol has a negative direct effect on the global radiation flux, while shortwave radiation reaching the surface has a weakened effect [16]. Consequently, the impact of natural and anthropogenic aerosols on NSSR has become a research focus in recent years. For instance, Wang et al. proposed a WRF-Chem model to simulate tropospheric aerosols in East Asia in March 2005 and found the direct radiation effect of aerosols leading to a 20 W/m 2 reduction in the surface shortwave radiation flux [19]. However, due to the influence of sources, differences in the radiation forcing caused by different types of aerosols exist on regional scale. The more concise approach by considering the aerosol conditions to retrieve NSSR is required urgently.
Many studies have demonstrated that the introduction of aerosol surface radiation estimation method improves the accuracy of retrieval algorithm, while the absence of aerosol radiation impact consideration causes higher errors [20], [21]. For instance, Ceamanos et al. introduced a method of estimating surface radiation (SIRAMix) for aerosol parameters by dynamically estimating the content and type of aerosols, which improved the retrieval accuracy of the downward radiant flux (DSSR) [22]. Generally, the aerosol radiative forcing (ARF) under the complicated atmospheric conditions is higher than the clear skies [23], [24]. Specifically, during super dust storm weather, the surface shortwave net aerosol radiative forcing (ASRF) could result in 144 W/m 2 reduction than non-dusty days [25]. Thus, precisely simulating the effect of aerosol interacted with the radiation balance of shortwave radiation should not be neglected. Furthermore, aerosols are affected by diverse emission sources in various regions, the composition of aerosol is relatively intricate, and aerosol type, aerosol morphology, particle size and vertical characteristics of aerosol particles are arduous to accurately distinguish and obtain [26]. Few studies have taken all the process into account for estimating NSSR, including solar zenith angle variations, absorption of water vapor, absorption of ozone, reflection of clouds, and absorption or scattering of aerosols [27]. Li et al. developed a parameterization model for retrieving NSSR, assuming that aerosols are relatively non-absorptive to atmospheric radiation, although it is unconscionable [3]. To date, a systematic understanding of how aerosol quantitatively contributes to NSSR is still lacking in the relevant research. Scholars have not determined an integrated set of retrieving algorithms of NSSR applicable to all surface types and various aerosol conditions.
The primary objective of this paper is to explore, analyze and quantify the net aerosol shortwave radiative forcing and to develop a highly accurate NSSR retrieval model under the influence of dust aerosol. Therefore, this paper assesses the significance of aerosol correction terms into the parameterization model by Tang et al. [5], in an attempt to investigate the retrieval of NSSR under aerosol conditions with both aerosol absorption and scattering processes considered. In quantifying the model parameters, the concept of Aerosol Optical Depth (AOD) is introduced, extinction coefficient, single scattering albedo and asymmetric factor are utilized to describe the optical properties of aerosols. Furthermore, the MODerate resolution atmospheric TRANsmission model (MODTRAN) 5 is used to simulate numerous atmosphere and geographic conditions [28]. Dust aerosol is considered for investigation in the study. The Moderate Resolution Imaging Spectroradiometer (MODIS) satellite products are chosen to provide accurate, multi-temporal and stable data sources, while ground-based station data are used to verify the accuracy of the new model.
The remainder of this paper is arranged as follows: Section 2 provides the basic data materials and Section 3 presents the proposed method for retrieving NSSR under aerosol conditions. Section 4 describes the validation results combined with the satellite data and ground measurements. The conclusions are drawn in Section 5.

A. SIMULATION DATA
In this work, a large majorities of numerical simulation datasets are constructed instead of processing large samples of the real data and remote sensing data. The atmospheric radiation transfer model, MODTRAN 5, developed by American Air Force Research Laboratory, is used to effectively simulate the reflectance, emission, transmission and scattering under the conditions of various surface landcovers, complicated atmospheric profiles, aerosols compositions, clouds, and observation geometries.
The spectral range of the simulation is set as 300nm-5000nm, the spectral resolution is 1 wavenumber (cm −1 ), and the corresponding wavenumber ranges from 2000 to 33333 cm −1 . In order to acquire shortwave radiation data representing various atmospheric conditions and underlying surface conditions at global scale, as well as MODIS VOLUME 9, 2021 channel radiance data generated at the top of the atmosphere (TOA), 18 profiles from the atmospheric profile library [29] (https://ara.lmd.polytechnique.fr/index.php?page=tigr) of Thermodynamic Initial Guess Retrieval (TIGR) database are selected for MODTRAN simulation, which is supplied by the Atmospheric Radiation Analysis (ARA) group at Laboratoire de Meteorologie Dynamique (LMD). The database contains five major categories of atmospheric conditions (tropical atmosphere, temperate atmospheric, cool temperate and polar atmosphere, in the northern hemisphere polar atmosphere, winter polar atmosphere) representing the total water vapor content (WVC) varied from 0.1 to 8.0 g/cm 2 . Six representative landcover reflectance spectra data as inputs are selected from the MODTRAN default database, including grassland, wetland, sandy loam, broadleaf forest, desert, and urban. Ten different visibilities varied from 3 to 21 km with an interval of 2 km are employed. Furthermore, six view zenith angles (VZAs) (0 • , 33.56 • , 44.42 • , 51.32 • , 56.25 • , and 60 • ), eight solar zenith angles (SZAs) (from 0 • to 70 • with an interval of 10 • ) and three relative azimuth angles (RAAs) (0 • , 60 • , 120 • ) are considered. Specifically, optical characteristics, extinction coefficient, single scattering albedo and asymmetric factor of dust aerosol are obtained from the Optical Properties of Aerosols and Clouds (OPAC) [30], and the radiation dataset simulating dust aerosol conditions in MODTRAN is involved. In this paper, the research scope is limited to the case where the aerosol optical depth (AOD) is less than 1 at a wavelength of 0.55µm, and the AOD is set from 0.1 to 1 with step length of 0.1. In total, 171,072 cases are simulated.

B. SATELLITE DATA
This work involves MODIS data on the Terra and Aqua satellites to verify the feasibility and accuracy of the NSSR estimation model. MODIS is a medium-resolution passive imaging spectroradiometer mounted on the Terra and Aqua platforms of the American polar-orbiting satellites of NASA's Earth Observing System (EOS) [31]. MODIS possesses a total of 36 bands, covering the spectral range from visible light to thermal infrared (0.4µm-14µm) with high confidence noise ratio products. Morning star (Terra) and afternoon star (Aqua) observe the earth once or twice a day, so that the global surface and atmospheric conditions can be consecutively detected. MODIS exists obvious superiority in spatial resolution, spectral resolution and calibration accuracy, which is widely utilized in various scientific research orientations including ocean, land, and atmosphere. Consequently, varieties of products of different levels could satisfy the needs of earth science research.
To derive the NSSR, five MOD series products with a resolution of 1 km are introduced. MOD02 is defined as TOA spectral radiation after radiation calibration, which includes three products, namely, MOD02QKM, MOD02HKM, and MOD021KM [32]. The MOD021KM product is selected in order to match the resolution of all data products when retrieving NSSR. The MOD03 products are the geoloca- tion parameters, including geographic coordinates, ground elevation, solar zenith angle (SZA), viewing zenith angle (VZA), viewing azimuth angle (VAA), and solar azimuth angle (SAA). The MODIS atmospheric precipitation product, MOD05 data, provides atmospheric precipitation data over the earth's surface, which considered as a dependent variable in the parameterization model. Furthermore, the cloud mask data supplied by MOD35 are identified with four different cloud confident coefficients including cloudy, uncertain clear, probably clear, and confident clear circumstance. On account of pixels of images identified with clear confidence have been adopted in our model, other circumstances are regarded as cloudy days. MODIS aerosol product MO04_3K (MODIS Terra/Aqua Aerosol 5 Min L2 Swath 3 km) is the level 2 aerosol product released by NASA, which can be used to obtain the aerosol optical properties from global marine and terrestrial environment of atmosphere, including aerosol type, aerosol equivalent diameter, asymmetric parameters, the Angstrom coefficient, etc. [33]. Additionally, MCD19-A2 [34] Level 2 grid dataset is an AOD product produced by Multi-angle Implementation of Atmospheric Correction (MAIAC) algorithm with 1 km spatial resolution, which is also considered to verify the results.

C. IN SITU MEASUREMENTS
In order to determine the accuracy and feasibility of the NSSR retrieval model considering aerosol changes proposed in this study, an area with a flat and homogeneous underlying surface type affected frequently by aerosol should be chosen as the assessment site. Three study areas with frequent occurrence of dust aerosols are chosen as verification sites with representativeness. The first station dataset contains the observation data of year 2007-2010 from the automatic meteorological station in Yingke irrigation area. The second one is Huazhaizi which is located at a piedmont desert at the middle reaches of Heihe hydrological and atmospheric observation network. The preceding two sites are located in Zhangye, Gansu province whose datasets span the period year from 2008 to 2011 as shown in Table 1. The Arvaikheer site is in Mongolia affected frequently by the influence of the nearby Gobi-desert. The sites in China are adjacent to the Badain Jaran Desert and Tengger Deserts, on the edge of the Gobi-desert, and are frequently disturbed by Gobi-desert sandstorms, with the most frequent and significant impacts occurring in northwest of China in spring [35]. The in-situ measurements in Gansu are obtained from the Watershed Airborne Telemetry Experimental Research [36] (WATER) and Hydrometeorological Observation Network (HiWATER) of the CAS (Chinese Academy of Sciences) Action Plan for West Development Project provided by the website of National Tibetan Plateau Data Center (http://data.tpdc.ac.cn) [37], [38]. Radiation data are observed by pyrgeometers mounted at a height of 4 or 6 meters above ground, respectively, with an observation interval of 10 minutes. The NSSR values are obtained by bilinear interpolation method from the measured values of 10 minutes before and after the satellite overpass local time which approximately represent the NSSR instantaneous values corresponding to MODIS pixels of the observation sites.

III. METHODOLOGY A. DEFINITION OF THE NET AEROSOL SHORTWAVE RADIATION FORCING
In order to investigate the contribution of aerosol on NSSR and to generate a NSSR retrieval model considering the influence of dust aerosol, this work firstly quantifies the net aerosol shortwave radiation forcing. Here ASRF is defined as the NSSR disturbance caused by aerosol depth influence, which is numerically equal to the difference value of NSSR in a certain aerosol type relative to the NSSR in the absence of the aerosol.
in which NSSR aero presents the NSSR under aerosols conditions, NSSR aero_abs is the NSSR corresponding to the absence of a certain aerosol type, the parameters in the NSSR aero_abs retrieval process are different from those in the clear sky model, which are calculated separately. ASRF is the net aerosol shortwave radiation forcing. According to the definition of ASRF, two datasets of NSSR and corresponding MODIS TOA spectral radiation simulation data are generated respectively in this study: one is radiation dataset of distinct aerosols under different AOD conditions, and the other is radiation dataset generated without the corresponding aerosol conditions. Each dataset contains NSSR data, channel TOA radiance data of MODIS sensor; supporting data comprise AOD data, aerosol optical properties data and atmospheric precipitation data.

B. NSSR RETRIEVAL METHOD UNDER THE AEROSOL CIRCUMSTANCE
Many studies have established that for different SZAs, a linear relationship always exists between the net shortwave radiation at the TOA and the NSSR at the surface [39]. Li et al. [3], on the basis of the linear relationship between Cess and Vulis [39], considered SZA as an explicit input parameter, and obtained a parameterized relationship of NSSR that could be calculated only by importing the SZA, atmospheric precipitation and shortwave broadband albedo at the TOA. This method has been widely applied in diverse satellite sensors with good retrieval results. Note that the surface absorption a s is defined as the function of SZA, WVC and the broadband albedo r. A linear relationship between the surface absorption a s and the broadband albedo r is expressed as follows: Studies argue that a s is expressed as a fraction of the flux incident at the TOA, which is defined by the ratio of the NSSR absorbed by the Earth's surface to the solar radiation at the TOA [40]. Where a s is defined by the equation: where µ denotes the cosine of SZA, E 0 is the solar irradiance at the TOA. D is the distance between the sun and the Earth, expressed in an astronomical unit. In equation (2), the intercept α and slope β are respectively expressed as follows: (4) and β = 1 + a 5 + a 6 ln µ + a 7 w z where a 1 -a 7 , and x, y, z are undetermined constant coefficients, w denotes the atmosphere precipitable water at a vertical column (g/cm 2 ). A linear parameterization developed by Tang et al. [5] is introduced in this study, which proposes a relationship between TOA narrowband radiance and the TOA broadband albedo. With the combination of different atmospheric and observation conditions, in order to accommodate the utilization of the narrowband reflectances ρ i , an algorithm conversion relationship from ρ i to the broadband albedo r is represented by (6): where b i (i = 1-5, 7) are the conversion coefficients, ρ 1 -ρ 5 and ρ 7 denote the TOA narrowband reflectances of MODIS bands, c 1i -c 4i are constants for a given SZA from 0 • to 70 • in the interval of 10 • respectively. A new lookup table (LUT) of c 1i -c 4i is reconstructed for calculating. Note that the SZA and RAA do not exist in equation (7), interpolation is introduced. Specifically, to establish relationship between obtained coefficients b 0 -b 7 , a polynomial function as equation (7) is used. ρ i (i = 1-5, 7) could be calculated from the formula given by: where L i is the spectral radiance integrated with the spectral response function of the MODIS bands,Ē band_i is the mean TOA solar irradiance derived from the spectral response function of MODIS satellite. The constant coefficients and the LUT in the equations above would be determined by the least square method. The variables in the model, including the independent TOA VOLUME 9, 2021 albedo r, the surface absorbed coefficient a s , and the narrowband spectral radiance L i , are all simulated in the radiative transfer model MODTRAN 5 under the cases combined by diverse conditions, including land covers, aerosol types, aerosol contents, atmospheric profiles, and different geometric observation conditions.
Nevertheless, the parameterization model for the clear sky conditions was extended by Masuda et al. [41] considering the impacts of the atmospheric properties involving the surface pressure, aerosol type, atmospheric water vapor, the ozone amount and the cloud type. Masuda declared that equation (2) should contain aerosol correction items, ozone correction terms and cloud correction items, in which the aerosol correction items of the coefficient as should be expressed as a function of SZA, AOD and broadband albedo r. However, our study persists that aerosol is not only related to absorption but also related to scattering. It is a complex interactive reflection process leading to positive or negative ASRF. Therefore, the NSSR reaching the surface should be analyzed synthetically according to the different types, characteristics and concentration of aerosol. Li and Masuda introduced an aerosol model named Haze 3, which was defined by Blanchet and List [42] as a moderately absorbing aerosol between the marine aerosol and the continental aerosol since it represented a natural aerosol type of Arctic Summer climate. Thus, this study will focus on a strong scattering aerosol type dust aerosol. Because of the above viewpoints, formulae for quantifying ASRF are obtained, which are involved in the subsequent analysis to determine the constant coefficients. ASRF = NSSR aero − NSSR aero_abs (9) If ASRF is positive, it means that the retrieved NSSR is underestimated in the previous clear sky retrieval model; otherwise, if the ASRF is negative, which means overestimation. According to the study by Lubin and Vogelmann [6], the particle size, optical depth, and relative humidity should be taken into account when aerosol radiation effect is involved in climate or meteorological modeling and it indicates that the aerosol type, optical depth, and atmospheric precipitation are important parameters to quantify impact. The variables which describe the ASRF will be discussed in Section 4.

A. CONSTRUCTION OF NSSR RETRIEVAL MODEL UNDER AEROSOL CONDITIONS
On the basis of the schema for ASRF estimation, the simulated datasets generated by MODTRAN 5 will be divided into two groups: the first group is the datasets of NSSR influenced by dust aerosol, and the other group is the datasets of NSSR generated in the absence of certain type aerosol.
As shown in Fig.1, the relationship between the mean value and standard deviation (SD) of ASRF under the conditions of dust aerosol with the change of different contents of AOD is presented. The solid point is ASRF, and half of the columnar line is the SD of the ASRF.   1 illustrates that ASRF changes from positive forcing to negative forcing with the increase of AOD, and the SD of ASRF increases with the increase of AOD. For NSSR, when AOD is between 0.1-0.3, the ASRF is positive, and between 0.3-1.0, the forcing is negative. The mean value of ASRF varies from -110 W/m 2 to 25 W/m 2 , and the standard deviation varies from 0 W/m 2 to 20 W/m 2 . Thus, it can be concluded that a linear relationship remains between ASRF and AOD. Hence, ASRF can be approximately parameterized using the following formula: where k and b are intercept and slope in Fig.1. The intercept and the slope are affected by atmospheric conditions and composition. Therefore, the changes of k and b with the atmosphere deserve to be investigated with emphasis. Primarily, ASRF is mainly driven by k value which is mainly affected by water vapor, CO 2 and other gases, such as polluted gases. During MODTRAN simulation process, the NSSR variation value is between 1-2 W/m 2 after doubling or halving CO 2 content. However, other polluted gases have small quantity and steady in the atmosphere composition, consequently they have few influence on the retrieval of NSSR than the WVC. The following part is mainly to analyze the influence of water vapor on values change of k and b.
Note that the change of k and b is closely related to the change of water vapor. The ASRF can be expressed by the exponential functional relationship between AOD and WVC. The model for estimating NSSR under clear conditions provided by equation (3) is a method combining parametric model and physical transmission model. The calculation process is complex, but the accuracy could be guaranteed. Thus, the k and b in equation (10) are generated respectively as follows: /d 4 ))) + d 5 ) (11) and b = e 1 × WVC 3 + e 2 × WVC 2 + e 3 × WVC + e 4 (12) where AOD, WVC are aerosol optical depth and atmospheric WVC, respectively, and d 1 , d 2 , d 3 , d 4 , d 5 , e 1 , e 2 , e 3 , and e 4 are regression constant coefficients which are shown in Fig.2. The parameters in the retrieval process are slightly different from those in the clear sky and are recalculated separately in this study. Fig.3 shows the relationship of root-mean-square error (RMSE) between the NSSR and the estimated NSSR for VZAs varied from 0 • to 60 • with an interval of 10 • under different surface types. When VZA = 0 • , the RMSE of the proposed model varies from 8.64W/m 2 to 13.76 W/m 2 . The RMSE tends to fluctuate as the VZA increases. The comparisons of actual NSSR with the NSSR derived by the proposed model is given below in Fig.4 with an observation VZA = 0 • . Fig.4 also shows the RMSE of the estimated value and the actual value of the model considering all types of the surface covers. The RMSE is 10.95 W/m 2 for all types under dust aerosols among 25,920 cases, which is nearly 42.90 W/m 2 improvement compared with the model before.

B. SENSITIVITY ANALYSIS
In order to determine the uncertainty of input parameters and the error of model estimation, sensitivity analysis is  necessary. The total error e( NSSR) of the proposed model contains the uncertainties of the input variables including the MODIS instrument detective noise, the AOD and WVC contributions, and the model accuracy, which is expressed as: where (AOD) and (WVC) are the uncertainties of input parameters AOD and WVC, respectively; δ(L NSSR ), δ(AOD), δ(WVC) are the errors caused by spectral radiance L i , (AOD) and (WVC), and δ(a lg) is the error of the model itself. The uncertainty WVC is 0.5 g/cm 2 for the MODIS  product. The AOD error expression of MODIS can be obtained based on validation test of MODIS aerosol product on a global scale as follows [43]: Fig .5 shows the errors caused by model, spectral radiance and aerosol increase accordingly with the increase of AOD, which indicates the model is more applicable at a higher AOD level. The algorithm error accounts for the majority of the model error, and the overall error e( NSSR) has an increasing trend with the increase of AOD. Theoretically, when AOD = 1, the model error can reach 40 W/m 2 . Hence, changes of aerosol have a great influence on the accuracy of the model. However, the errors caused by multiple input parameters may be counteracted, the actual error of the model may be smaller than the theoretical value illustrated in the figure.

C. VALIDATION WITH THE METEOROLOGICAL STATION DATA
In this work, aerosol optical depth is obtained to participate in model estimation of NSSR. Fig.6 shows the aerosol retrieval values of Yingke, Avaikheer, and Huazhaizi sites. AOD data required for model validation are obtained from MODIS Level 2 aerosol products MOD04_3K/MCD19-A2, with spatial resolutions of 3 km and 1 km, respectively. The aerosol optical depth observed by MODIS is the average value in a region, and the single point value corresponding to the observation site, so it is necessary to match the MODIS observation data with the ground site in a spatial and temporal scale, simultaneously ensure at least five MODIS effective pixels when matching. The MODIS aerosol product data are taken as the average value of satellite observation data within the latitude and longitude range of 0.1 • × 0.1 • centered on the in-situ sites. We used the bilinear interpolation of the measured values of 10 minutes before and after the satellite overpass time to obtain the NSSR data on the ground. Since NSSR changes little in a short time, the value approximately represents the NSSR value matching with the MODIS pixels of the observation site. Cloudless pixels are selected to eliminate the influence of cloud in the NSSR retrieval model. Additionally, to determine whether dust aerosol exists in the pixel, the pixel with aerosol in the center pixel and surrounding pixels at the location of the observation station are selected as the target pixels.
As for the other input parameters in the model, the TOA radiance is taken from the MOD021KM product, which is a Level 1 product that has undergone radiation correction. WVC is derived from the water vapor product MOD05_L2. VZA and other geographic data are obtained from MOD03 product, and cloud mask data are obtained from MOD35 product. The range of AOD in Yingke station varies from 0.17 to 1.31 at the year of 2008-2010 and the WVC varies from 0.29 to 1.81 g/cm 2 . The range of AOD in Arvaikheer station varies from 0.01 to 1.05 at the year of 2003 and the WVC varies from 0.36 to 2.27 g/cm 2 . The AOD of Huazhaizi is from 0.09 to 1.25 and WVC is from 0.37 to 4.0 g/cm 2 , respectively. Fig.6 illustrates the comparisons of the NSSR with estimations using the proposed model (red solid points) and NSSR cl model before improving (blue hollow circles) with the VZA = 0 • in different sites. NSSR cl represents the retrieved NSSR without considering the dust aerosol using conventional parameterized model. In order to evaluate the application ability of the two models in the dust aerosol weather, scatter plots of AOD and WVC with model retrieval errors are given in Fig.7 and Fig.8, respectively.
There exists 21 days of matching data for MOD04 product and 12 days for MCD19 product at Yingke site from year 2008 to 2010 in spring and autumn seasons. The RMSEs between the measured and retrieved NSSR with two AOD products respectively are both 160 W/m 2 as shown in Fig.6 (a) and (b). The biases are 43.5 W/m 2 for the MOD04 products and 84.6 W/m 2 for MCD19 products.  As shown in Fig.7, the boundary values of relative errors utilizing MOD04 in Yingke site exist an improvement of 52.7 W/m 2 compared with the minimum NSSR relative errors of NSSR cl model, and of 82.8 W/m 2 reduction compared with the maximum NSSR relative errors of NSSR cl model. Combined with three-year AOD products, the proposed model can improve the NSSR retrieval accuracy of 43.9 W/m 2 and 42.8 W/m 2 for MOD04 and MCD19 products in dust weather, respectively.
At Arvaikheer site, dust aerosol weather occurs mostly in spring and early summer. The RMSE of NSSR aero from the    Fig.10.
Although the validation sites are adjacent to the desert, dust aerosols are mostly presented in spring and autumn. Thus, aerosol types and contents have seasonal characteristics. The RMSE of the proposed model is obtained by comparing the NSSR calculated by satellite data with the surface measured data under the condition of dust aerosol. However, the RMSEs of the model in the measured data are much larger than ones in the simulated data of MODTRAN. The accuracy may be affected by many uncertain factors. Firstly, water clouds and ice clouds affect the direct and diffuse parts of shortwave radiation to some extent. What's more, aerosol particles are usually regarded as spherical ideally in the simulation, and the deviation of the actual situation may cause errors in the retrieval model. Additionally, dust particles are typically non-spherical, which may affect the RMSE of the proposed model. Shortwave radiation in complex terrain may increase secondary reflections from adjacent slopes, which increases the uncertainty of the model. Besides, the verification process assumes that NSSR value will not change in a short time, which may influence the model accuracy.
Furthermore, the accuracy is affected by the uncertainties of AOD. The profile data of aerosol which describe aerosol optical properties in different layers have not been taken into account in the validation part. Lidar is an effective and reliable method to detect the spatio-temporal distribution of optical physical parameters with high precision compared with the traditional method. CE318 is the sun-photometer which can calculate the transmittance and scattering characteristics in order to measure the aerosol characteristics. However, the research sites do not contain the data above, which causes the uncertainty in the model. Aerosol source, complex composition, particle size distribution and the rapid spatio-temporal variation are different from the ideal conditions leading to the impact of the final retrieval results.
In conclusion, our proposed model performs a good accuracy at the dust aerosol days. When AOD and WVC are qualified our proposed model can improve the retrieval accuracy of NSSR under the influence of dust aerosol and reduce the uncertainty in application of shortwave radiation. When the AOD is less than 0.3, the results of our proposed model are close to those values derived from the model of clear sky, indicating that ASRF contributes little under these conditions. In this paper, the instantaneous NSSR is studied using MODIS data from the overpass satellites in the morning and afternoon, but no time scale extension is involved in. In the future studies, time scale extension methodology would be investigated to infer the change curve of NSSR in the whole day.

V. CONCLUSION
Surface energy balance is a significant part of global change research and is one of the hot issues in remote sensing research field. A model has been generated to estimate instantaneous NSSR under the aerosol days in this study. We produced a model containing an ASRF correction quantified by a function of AOD and WVC, which can weaken the direct and indirect effects of aerosol interacted with radiation. The RMSEs of the proposed model vary from 8.64W/m 2 to 13.76 W/m 2 with the VZA = 0 • under the dust aerosols for different land types, while for VZA = 60 • change from 7.02W/m 2 to 12.06 W/m 2 .
Meteorology stations data measured in Gansu Province, Zhangye, were utilized to validate results of our proposed model. The results indicated that the proposed model could improve the retrieval accuracy of NSSR. The RMSE of Yingke station is 160 W/m 2 which has improved about 40 W/m 2 among dust aerosol days compared with the results retrieved under the clear skies. Similar conclusions are obtained by employing the aerosol products from MOD04 and MCD19-A2 derived from MODIS product data. For Avaikheer and Huazhaizi stations, the proposed model has increased the accuracy for 14.2 W/m 2 and 17 W/m 2 , respectively. The retrieval accuracy has been improved about 20 W/m 2 for the year-round cloudless weather (including aerosol and no aerosol weather) under the atmospheric condition at Yingke station and Huazhaizi station.
It should be noted that the aerosol vertical profiles, heightdependent distribution have not been considered in this study. The vertical distribution profile of aerosol could be measured by lidar and solar photometer. Therefore, the effect of aerosol vertical profile changes on NSSR would be considered in the following subsequent studies. The characterization level of parametric models at different sites is different, resulting in different accuracies. Therefore, some more ground sites should be added to study the universality of the NSSR retrieval model applicable to the influence of dust aerosol.