Evaluation of aerosol optical depth and aerosol models from VIIRS retrieval algorithms over North China Plain

The first Visible Infrared Imaging Radiometer Suite (VIIRS) was launched on Suomi National Polar - orbiting Partnership (S - NPP) satellite in late 2011. Similar to the Moderate resolution Imaging Spectroradiometer (MODIS), VIIRS observes top - of - atmosphere spectral reflectance and is potentially suitable for retrieval of the aerosol optical depth (AOD). The VIIRS Environmental Data Record data (VIIRS_EDR) is produced operationally by NOAA, and is based on the MODIS atmospheric correction algorithm. The “MODIS - like” VIIRS data (VIIRS_ML) are being produced experimentally at NASA, from a version of the “dark - target” algorithm that is applied to MODIS. In this study, the AOD and aerosol model types from these two VIIRS retrieval algorithms over the North China Plain (NCP) are evaluated using the ground - based CE318 Sunphotometer three (Beijing), spectral AOD. For three VIIRS_EDR at 550 a positive mean bias (MB) of 0.04 - 0.06 the correlation of 0.83 - 0.86, with the largest MB (0.10 - 0.15) in In VIIRS_ML 550 positive of 0.13 - 0.14 a higher correlation - 0.94) CE318 aerosol model each well aerosol aerosol model used VIIRS_EDR dominant model evaluation overall accuracy rate aerosol model VIIRS_ML over NCP three sites VIIRS_EDR differences in Single Scattering Albedo (SSA) at 670 nm VIIRS_ML CE318 than 0.015, high seasonal differences are found especially over the Xinglong site. The values of SSA from VIIRS_EDR are higher by CE318 over all sites and all assumed aerosol modes, with a positive bias of 0.02 - 0.04 for fine mode, 0.06 - 0.12 for coarse mode and 0.03 - 0.05 for bi - mode at 440nm. of SSA but positive AOD MB of VIIRS_EDR other surface reflectance characterization or cloud contamination) are important sources of error in the VIIRS_EDR algorithm, and their effects on aerosol retrievals may override the effects from non - ideality in these aerosol models.


Introduction
Atmospheric aerosols have important impacts on climate, air quality and human health [1][2][3]. Their properties are highly variable in both space and time. Space-based platforms provide a global view of the aerosol system, unmatched by any other measurement system in terms of the spatial coverage [4]. With the long-history of aerosol products derived from the Moderate resolution Imaging Spectroradiometer (MODIS) and the aging of the instrument, there is programmatic interest in continuing similar aerosol retrieval capabilities. Since the launch of the Visible Infrared Imaging Radiometer Suite (VIIRS) instrument on the Suomi National Polar-orbiting Partnership (S-NPP) satellite in late 2011, there has been great interest in retrieving aerosol properties from VIIRS [5][6][7].
Currently, there are multiple algorithms available for deriving aerosol optical depth (AOD) and other aerosol properties from VIIRS data [6,8]. Here we consider two. The VIIRS Environmental Data Record (VIIRS_EDR) is being produced by the United States National Oceanic and Atmospheric Administration (NOAA) [6]. At the same time, the National Aeronautics and Space Administration (NASA) is considering long-term continuity for developing an aerosol climate data record. For this purpose, Levy et al. are experimenting with a "MODIS-like" dark-target algorithm for use on VIIRS data (VIIRS_ML) [8]. How does each of these algorithms perform for retrieving AOD and other aerosol properties over China, a region of extreme diversity of aerosol sources, compositions, and loadings? Preliminary evaluation of VIIRS_EDR and VIIRS_ML derived AOD, has been performed separately using co-located sunphotometer data from the AErosol Robotic NETwork (AERONET) and other networks [10]. For the period of 1 May 2012 -30 April 2013, Jackson et al. showed VIIRS_EDR underestimated the AOD over land with global mean bias of -0.02 [6]. However, the Suomi NPP VIIRS aerosol data product assessment report showed VIIRS_EDR overestimated the AOD in land of 0.06 to 0.11 for the period of 1 May 2012 -14 Oct 2012 [11]. In addition, larger biases were found in western Asia and India [10]. Since experimental VIIRS_ML has only been available since 2015, the evaluation is limited. Levy et al. briefly validate the VIIRS_ML AOD at 550 nm by comparing it to AERONET observations from March 2013 to February 2014 and showed VIIRS_ML overestimates the AOD over land with global positive bias of 0.005 [8].
The above studies have limited scope in that they only provide estimates of the global expected error over land. They do not focus on evaluating products over regional scales, especially where AERONET data are sparse. Hence, the focus of this study is to evaluate VIIRS_EDR and VIIRS_ML data over the North China Plain (NCP: 114-120°E, 34.5-41°N), one of the most densely populated regions in China that has experienced enormous economic growth in past two or three decades [12][13][14][15]. Indeed, the NCP is one of the most severely polluted areas in the world with frequent heavy haze events in recent years [12,16,17]. Given that retrieved aerosol optical properties are often used as a proxy for assessing climate and air quality in the NCP region, a regional validation of VIIRS aerosol products has important consequences [13,18,19].
The organization of this paper is as follows. We introduce the VIIRS instrument and the two aerosol retrieval algorithms and corresponding data in Section 2. The validation data set and methods for inter-comparison are in Section 3. The AOD evaluation results are presented in Section 4. Section 5 states the aerosol model types and the optical properties comparison between the CE318 sunphotometer and VIIRS, including VIIRS_EDR and VIIRS_ML. The conclusions are presented in Section 6.

VIIRS Satellite data 2.1 What is VIIRS?
VIIRS is a cross-track scanning radiometer with 22 spectral bands covering the visible/infrared spectrum from 0.412 to 12.05 µm. The design and concept of VIIRS operations combine aspects from several legacy instruments, including the NOAA's Advanced Very High Resolution Radiometer (AVHRR), NASA's MODIS, Sea-viewing Wide Field-of-view Sensor (SeaWiFS), and the Department of Defense's Operational Linescan System (OLS) sensors [20]. It has a wider swath (~3000 km) than MODIS, which allows a global sample of the Earth everywhere every day. It flies in a Sun-synchronous near-circular ascending polar orbit 829 km above the Earth with the local equator-crossing time at 13:30.
VIIRS has three types of bands: imagery bands (I-bands), moderate resolution bands (M-bands), and the day-night band [6]. The M-bands (total 16 bands) have 0.742×0.776 km nadir resolution and 1.60×1.58 km at the edge of scan. Other bands are used to create the VIIRS Cloud Mask (VCM), which is used as input to aerosol algorithms, as well as in internal tests to characterize environmental conditions. Most of M-bands are used to derive the aerosol parameters.

Overview of the two retrieval algorithms over land
The VIIRS_EDR and VIIRS_ML algorithms are similar in many ways. Both algorithms start with the satellite measurements of spectral reflectance at the top-of-atmosphere (TOA), and are compared to a look-up table (LUT) to determine the most plausible solutions for aerosol and surface properties. The measured reflectance at the TOA is a summation of scattering events from the surface and the atmosphere. In both algorithms, the aerosol optical properties of the aerosol models are essential for radiative transfer computing to generate the atmospheric LUT that is needed for AOD retrievals.
However, there are also differences. The VIIRS_EDR algorithm has the heritage from the MODIS atmospheric correction algorithm for land surface reflectance, in which the expected surface reflectance ratio at different wavelengths are used as a constraint in retrieving AOD [21,22]. In contrast, VIIRS_ML has the heritage from the MODIS Collection-6 aerosol algorithm, in which surface reflectance ratios at different wavelengths are prescribed and the TOA reflectance is used as the most important constraints [8,9]. Aerosol model type assignment in the two algorithms is also different. In VIIRS_EDR algorithm, the aerosol model type is selected at each pixel for each inversion by using extra blue wavelengths to constrain the aerosol type, while in the VIIRS_ML algorithm, the aerosol model type is assigned to each region and each season prior to retrieval based on the past cluster analysis of AERONET inversions [9,23]. In our assessment of the VIIRS AOD in the NCP region, we will also use aerosol properties from CE318 inversions to evaluate the aerosol model types and associated optical properties in the VIIRS_EDR and VIIRS_ML algorithms, and thereby, analyze one likely source for AOD retrieval uncertainties.
Based on the climatology of AERONET inversion data, the VIIRS_EDR algorithm defines a set of five microphysical aerosol model types. These five models are denoted as dust (for example, observed at Cape Verde), high absorption smoke (African savanna, Zambia), low absorption smoke (Amazonian forest, Brazil), clean urban (Goddard Space Flight Center, Greenbelt, MD), and polluted urban (Mexico City) aerosols. All model types have size distributions defined by bimodal lognormal distributions of spherical particles [6]. As explained by Jackson et al. [6], the retrieval LUT is created by starting with a Mie scattering code, for which aerosol inputs include a real part and an imaginary part of refractive indices and size parameters of aerosol fine and coarse mode (i.e. volume mean radius, standard deviation and volume concentration). During the retrieval, the algorithm selects the aerosol models with the lowest residual which is computed based on deviations between the 412 nm, 445 nm, 488 nm, and 2250 nm surface reflectances predicted from the 672 nm surface reflectance and the computed surface reflectances using the retrieved AOD for that model type.
For the VIIRS_ML, aerosol model types are also derived from AERONET inversion climatology. However, the retrieval algorithm uses that information in a different way. Levy et al. clustered and classified AERONET retrieval products into statistics that represented the most likely aerosol conditions for a particular region and season [23]. These aerosol model types are separated into fine-mode dominated (fine model) and coarse-mode dominated (coarse model), and the fine model is further separated into being strongly absorbing, moderately absorbing and weakly absorbing aerosol models. In the classification, the moderately absorbing aerosol model is set as the default, overwritten only if clear dominance of one of the other two aerosol model types is observed. By clustering, it is shown that the single scattering albedo (SSA) values at 670 nm of three fine models is ~0.85 for strongly absorbing, ~0.90 for moderately absorbing and ~0.95 for weakly absorbing and ~0.95 for the coarse model. The global type classification was updated for Collection-6, by classifying the AERONET data through 2010 [9,23]. Note that the categories of aerosol model type used for VIIRS_ML are not exactly analogous to those used for VIIRS_EDR.

VIIRS_EDR data
The VIIRS_EDR level 2 aerosol products are obtained from NOAA Comprehensive Large Array-data Stewardship System (CLASS) at http://www.nsof.class.noaa.gov. The VIIRS_EDR aerosol parameters are derived primarily from the M-bands of the radiometric channels covering the visible through the shortwave infrared spectral regions (412 nm to 2250 nm). As explained by Jackson et al. [6], the VIIRS_EDR AOD is generated from 8×8 pixel aggregations of the intermediate product (IP), where in turn the IP represents retrieved AOD for each and every native resolution (e.g. 0.75 km) pixel. The pixels with clouds, cloud shadows, snow, ice, subpixel water, bright land surface, fire, sunglint, suspended sediments or shallow water, and large solar zenith angle are screened out using the internal tests and the external VIIRS VCM before proceeding with the aerosol retrieval [6]. The VIIRS_EDR product represents the statistics of the 8×8 aggregation, which is a retrieve then average strategy. Consequently, the resolution of the VIIRS_EDR data is ~6×6 km 2 at nadir (~12.8×12.8 km 2 at the edge of scan). Data screening and aggregation methods can be seen in [6]. VIIRS_EDR AOD has been collected from 2 May 2012 to 31 March 2014 over the NCP. The data between 15 Oct 2012 and 27 Nov 2012 are rejected because an inadvertent error introduced in the operational aerosol code during this period [6]. The AOD at 488 nm, 550 nm and 672 nm and AOD Quality Flags (QF1), as well as Land Model Aerosol Index flag (QF4), are used in this study. The values of QF1 refer to the estimated "quality" of the retrieval product, so that QF1= 0, 1, 2, and 3, represent not value produced, low, medium and high quality, respectively. The values of QF4 refer to which aerosol model type was used in the AOD retrieval, where QF4= 0, 1, 2, 3, and 4, refer to dust, high absorption smoke, low absorption smoke, clean urban, and polluted urban aerosol model, respectively.

VIIRS_ML data
VIIRS_ML data are available from the NASA Atmosphere Science Investigator-led Processing System at the University of Wisconsin (A-SIPS; http://sips.ssec.wisc.edu/). Following the strategy of the Dark-Target retrieval, the VIIRS_ML follows an average, then retrieve once logic [8]. This means that the averaging is upon observations (spectral reflectance) within the box, and following the MODIS protocol, the aerosol retrieval is performed only once. Using 10×10 pixel aggregations, the VIIRS_ML aerosol product is reported at 7.5 km (at nadir) resolution based on M-band pixel resolution. The VIIRS_ML algorithm does the cloud masking by applying the internal spatial variability and reflectance threshold tests (e.g. 3×3 pixel spatial variability and visible/1024/1038 nm tests). The strategies for masking, selecting and aggregating pixels for VIIRS_ML are described in [8]. As for the aerosol model type, Levy et al. note that the aerosol model in NCP region is assumed to be moderately absorbing fine model during Winter (DJF) and Spring (MAM), and weakly absorbing fine model during Summer (JJA) and Autumn (SON) [9]. Since the aerosol type may differ day-to-day, this assumption is meant to be climatologically representative and can lead to errors in instantaneous AOD retrieval. These same assumptions are used for the VIIRS_ML.

Ground-truth data and methods for satellite-sunphotometer comparison 3.1 Sunphotometer data
The ground data used to evaluate the VIIRS aerosol products in this study consists of CE318 sunphotometer (CE318) observations and retrievals. The CE318 instrument performs direct sun extinction measurements at 8 wavelengths ranging from 340 to 1020 nm and sky radiance measurements at 4 wavelengths, i.e., 440, 675, 870, and 1020 nm, respectively. The AOD data were calculated from direct sun observations with an accuracy of 0.01 to 0.02 [25,26]. Refractive index, volume mean radius, volume concentration and single scattering albedo (SSA) retrieved from the CE318 sky measurements characterize the aerosol type. The uncertainties of refractive index are 30-50% for the imaginary part and 0.04 for the real part when AOD at 440 nm (AOD440nm) > 0.4 and solar zenith angle > 50°, and the uncertainties increase for lower AODs [24,27]. SSA uncertainty is estimated to be less than 0.03 for AOD440nm > 0.4 and the uncertainty increases for lower AODs [24,27]. Note that inversion data (size/optics) are sparse compared to direct sun observations of spectral AOD.
CE318 sunphotometer data (including the corresponding inversion products) over three sites in NCP region during the period of 2 May 2012 -31 March 2014 were used. The location and description of the three CE318 sunphotometer sites are provided in Table 1. These three sites can be considered as representative of urban (Beijing), suburban (XiangHe) and regional background (Xinglong) environments, respectively. Affected by Asian monsoons, the NCP region has a moderate continental climate with cold winters and hot summers. Heavy anthropogenic pollution from urbanization, industrial, and agricultural activities mixed with coarse dust particles (most occurring in spring) result in a rather complex nature of aerosol physical and optical properties in the NCP [13]. Notably, the regional background station Xinglong is located at a mountain with the elevation of 970 m which is higher than the other two sites. However, even at this station, urban/industrial and dust aerosol could occur through aerosol regional transportation [28] and secondary aerosol formation. Therefore, as will be shown in our analysis, the complex features of aerosol properties may help explain some of the uncertainties in satellite retrievals of AOD in this region. Data at Beijing and XiangHe are downloaded from AERONET (http://aeronet.gsfc.nasa.gov/) and the data at Xinglong are obtained from China Aerosol Remote Sensing Network (CARSNET) [29]. AERONET level 1.5 inversion data (from sky-light measurements) are used since the level 2 inversion data are very less frequent and unsuitable for data statistics. At the same time, we used the conditions of AOD440nm > 0.4 and solar zenith angle > 50° to constrain the data quality according to [24,27]. The AOD in these two networks are consistent with one another as the correlation coefficients are larger than 0.999 and have a 99.9% significance level [29]. The CARSNET calibration and comparison with AERONET are described in detail in other references [29][30][31].
It is noted that Beijing and XiangHe are part of AERONET stations, and their aerosol data during 2005 were used in the cluster analysis for the VIIRS_ML aerosol model assignment, but the Xinglong site was not used [23]. Moreover, the retrieved aerosol properties in recent years may change from those used in 2005 and 2010 due to the rapid development in the past few years over the NCP region [12,13]. Therefore, using aerosol property data from more ground sites during recent years over this region can be used to help evaluate whether those past analyses from shorter periods (e.g. only one year) and fewer sites are still representative. Notably, none of these three sites were used in the analysis by Dubovik et al. [24], which means they do not characterize typical aerosol properties of smoke, dust, and urban particles that have been adopted in the VIIRS_EDR algorithm. Thus, these three sites are better suited to evaluate the aerosol model type used in both the VIIRS_EDR and VIIRS_ML algorithms over the NCP region.

Method for data matchup
The spatiotemporal collocation between satellite and CE318 measurements follows the method of the Multi-sensor Aerosol Products Sampling System (MAPSS), in which sunphotometer data with ±30 minutes of satellite overpass are compared with satellite data within 25 km radius of the sunphotometer [32]. Minimum requirements for a matchup are at least 2 observations from AERONET and 5 pixels from the satellite. The CE318 AOD at 550 nm and at VIIRS blue (488 nm) and red (672 nm) bands are interpolated from 440 nm, 675 nm, 870nm and 1020 nm by using an established fitting method [33]. The results for comparison of AOD values between VIIRS (VIIRS_EDR and VIIRS_ML) and ground CE318 observations are presented with various statistical parameters, including the number of matchup data (N), the mean bias (MB), root mean squared error (RMSE), correlation coefficient (R), and the percentage of data within the expected error 0.05+0.15AOD (%EE) which is used as the MODIS AOD expected uncertainty over land [34], the slope (Slope) and intercept at y-axis (y-int) of linear regression.

Methods for aerosol model evaluation and aerosol properties analysis
Due to constraints placed on the inversion of CE318 sky-radiance data (AOD440nm > 0.4; solar zenith angle > 50°, etc.), statistics for aerosol optical properties are sparse. To collocate aerosol optical properties from sunphotometer with aerosol models assumed by either VIIRS retrieval algorithm, we require different averaging domains. Here, we use daily-averaged aerosol optical properties retrieved from CE318 sky radiance measurements. Since the aerosol model types used for satellite AOD retrievals may vary spatially, we select only the model type assumed at the pixel that includes the site.
In the extraction of AOD from the VIIRS_EDR, those with quality QF<1 retrievals are rejected (these with QF<1 are not products and are mostly with cloud contamination and sunglint). Seasonal and total frequencies of each aerosol model type occurrence in the three sites are calculated to show the typical aerosol model types used in the VIIRS_EDR land algorithm over NCP sites.
The aerosol model type evaluation of VIIRS (VIIRS_EDR and VIIRS_ML) is based on the SSA comparisons between VIIRS and CE318. The SSA values at four wavelengths (440 nm, 670 nm, 870 nm and 1020 nm) of the five aerosol model types used in VIIRS_EDR can be obtained from Dubovik et al. [24]. The values of SSA at 670 nm (SSA670nm) are 0.98, 0.84, 0.93, 0.97 and 0.88 for dust, high absorption smoke, low absorption smoke, clean urban and polluted urban aerosol model, respectively. Thus, the SSA values at the NCP sites in the VIIRS_EDR retrival can be derived from the aerosol model type assumed in the VIIRS_EDR pixel that includes the site. The VIIRS_ML aerosol model type is assumed globally based on the cluster analysis of SSA670nm derived from all AERONET inversions and it is fixed in each season for the NCP region (i.e. weakly absorbing fine model for Summer and Autumn: SSA670nm ~0.95, moderately absorbing fine model for Spring and Winter: SSA670nm ~0.90) [9,23], so it is unnecessary to extract the aerosol model from VIIRS_ML pixel-by-pixel. Thus, we do the seasonal comparison, that is the seasonal mean SSA670nm values of CE318 inversion are calculated and compared with the seasonal SSA670nm values of the VIIRS_EDR and VIIRS_ML. Only the VIIRS SSAs with the AOD550nm > 0.25 are used to meet the requirement of CE318 AOD440nm > 0.4 [23].
The SSA is also used to classify the aerosol type of the CE318 inversion, which is to evaluate the aerosol model type for each retrieval from the VIIRS algorithms. We firstly collocate the daily matchup data between the VIIRS_EDR and CE318 and between the VIIRS_ML and CE318. To evaluate the aerosol model type of VIIRS_EDR, the CE318 inversions are classified to the five aerosol types as the VIIRS_EDR. The CE318 inversion with Angstrom exponent <0.6 and AOD at 1020 nm >0.3 (according to Dubovik et al. [24]) is classified as dust type. If not dust type, the CE318 inversion with SSA670nm <0.86 is the high absorption smoke, with 0.86< SSA670nm <0.905 is polluted urban, with 0.905< SSA670nm <0.95 is low absorption smoke, and with SSA670nm >0.95 is clean urban aerosol type. As for the evaluation of aerosol model type used in the VIIRS_ML, the CE318 inversions are classified to the four aerosol types as the VIIRS_ML. Use the same way to find out the CE318 data with coarse model (same as the dust). For the rest CE318 data, that with SSA670nm <0.875 is regarded as strong absorbing fine model, 0.875< SSA670nm <0.925 is the moderately absorbing fine model, and SSA670nm >0.925 is the weakly absorbing fine model. The threshold values of the classification are based on the SSA values of the aerosol models used in the VIIRS algorithms. This method is actually using the aerosol size and scattering properties to classify the aerosol type, which has been studied by Giles et al. [35]. After classifying the CE318 aerosol type, the comparisons of aerosol type between the VIIRS_EDR and CE318 and between the VIIRS_ML and CE318 are done to show the accuracy rate of aerosol model type used in the VIIRS_EDR and VIIRS_ML. If the VIIRS aerosol model type is same as that of CE318, the aerosol model type used in the VIIRS is deemed as accurate. The accuracy rate is defined as the ratio of the number of accurate to the number of all daily matchups. The accuracy rate reflects the applicability of aerosol model type used in the VIIRS algorithms.
We also conduct the SSA comparison of different modes (fine, coarse and bi-mode) between the VIIRS_EDR and CE318 retrieval. We compute the SSA for all aerosol modes at 440 nm in the VIIRS_EDR by inputting the aerosol parameters (refractive indices, size parameters and volume concentrations of each mode) into Mie scattering calculation [36]. The reason for using 440 nm is that the aerosol model properties in the VIIRS_EDR algorithm are mostly referred at 440 nm [6]. The input aerosol parameters of the VIIRS_EDR at each site are obtained by extracting the aerosol model type over the site pixel and calculated according to the Table 2 in reference [6]. Although SSA440nm is available from CE318 inversion, for consistency we also use the Mie code to compute the SSA440nm based on the aerosol optical properties inversed from CE318 sky radiances. We have compared the SSAs between the CE318 inversion and the Mie scattering calculation and the result shows that the bias of the two SSA values is very low (less than 0.01). That is because the CE318 inversion also uses Mie scattering calculation to obtain SSA. The resultant SSA440nm of CE318 is compared with that of the VIIRS_EDR SSA440nm. Since the CE318 sky-radiance inversion product is only reliable for AOD440nm>0.4, we also choose the aerosol properties of the VIIRS_EDR when AOD550nm>0.25 [23].   The MB and RMSE in the NCP are both larger than the counterparts in the global assessment statistics [6,10]. Table 3 presents the evaluation results of each NCP site, separately. Clearly, all properties (MB, RMSE, R and %EE) demonstrate the worst performance over Beijing. The MB in XiangHe site is (-0.02)-0.00 with high R of 0.89-0.92, which is more comparable to the global agreement [10]. Beijing is an urban site, while the VIIRS_EDR uses the global surface reflectance ratios as the expected spectral surface reflectance relationship, which may cause the largest error at the Beijing site [6].

Results of AOD inter-comparison 4.1 Evaluation of the VIIRS_EDR AOD at 550 nm
In the previous global evaluation for VIIRS_EDR products, it was recommended to use data with a higher QF [10]. In the NCP, comparisons at both XiangHe (suburban) and Xinglong (rural) support this recommendation (R and %EE increase and MB decreases with increasing QF). However, for the Beijing site, the matchup statistics are poor, and increasing QF does not help. Thus, for sites that are not optimal for aerosol retrieval in the first place (e.g. urban), QF may not be a useful diagnostic.

Evaluation of the VIIRS_ML AOD at 550 nm
Comparing the VIIRS_ML to CE318 ( Table 2), there are 817, 755, and 683 matchups for QF >0, >1, and =3, respectively. Since each algorithm has its own definition of QF, there are different relative contributions of each QF level [6,8]. Overall, the VIIRS_ML AODs show a high correlation with CE318 (R is 0.93-0.94) but overestimates over this region with high MBs (0.13-0.14) and large slope values for the equations of best fit (1.26-1.27). More than half of the VIIRS_ML data are within the expected error envelope (%EE>50%). As compared to VIIRS_EDR, the VIIRS_ML has a higher bias, but has larger correlation with more data within the EE. Following site-by-site comparison for the VIIRS_EDR, we evaluate the VIIRS_ML AOD at each site ( Table 4). Like the VIIRS_EDR, the VIIRS_ML shows the largest bias over the Beijing urban site. Since the urban surface reflectance may be underestimated in the "dark-target" algorithm, this can lead to an overestimation of AOD [37]. The VIIRS_ML AOD over XiangHe performs the best with the highest values of R and %EE but the MB is not the lowest and it is higher than that between the VIIRS_EDR and CE318 (MB of (-0.02)-0.00). However, the R values between the VIIRS_ML and CE318 at all three sites are higher than those found for the VIIRS_EDR and CE318. The values of %EE of the VIIRS_ML over the XiangHe and Xinglong sites are higher but lower over the Beijing site compared to those between the VIIRS_EDR and CE318.
The results of the quality flag analysis of the VIIRS_ML AOD are similar to those for the VIIRS_EDR. Using high quality data leads to the best performance at the XiangHe and Xinglong sites but it is not suitable at the Beijing site.

Evaluation of the VIIRS AOD at red and blue bands
While the AOD is retrieved at 550 nm, neither algorithm uses reflectance at 550 to derive AOD. This is because the Earth's surface tends to be brighter in green wavelengths (e.g. vegetation), and not suitable for aerosol retrieval. The VIIRS_EDR algorithm is based on the calculation of surface reflectance at blue (488 nm) and red (672 nm) and three other bands (412 nm, 445 nm and 2250 nm) [6]. As for the VIIRS_ML algorithm, AOD at 550 nm is inversed by using 488nm, 672nm and 2257nm measured TOA reflectance to find a the optimal solution [8]. Hence, the VIIRS AOD spectral dependence between blue and red bands are also evaluated with the CE318 data. According to quality flag analysis in Section 4.1 and 4.2, we choose the matchup data of QF>1 for both the VIIRS_EDR -CE318 and the VIIRS_ML -CE318. Then we selected the data of VIIRS_EDR, VIIRS_ML and CE318 with the dates in common, that was selecting the matchup data of VIIRS_EDR, CE318 and VIIRS_ML. Figure 1 shows the average AODs of the matchup data between VIIRS_EDR, CE318 and VIIRS_ML at three wavelengths over the three sites. Compared to CE318 measurements, the VIIRS_EDR AOD at 488 nm overestimates over Beijing and Xinglong but performs well for the XiangHe site. As 672 nm, the VIIRS_EDR AOD also overestimates over Beijing, but slightly undervalues AOD at XiangHe. However, the VIIRS_ML AOD overestimates at all the three wavelengths and over all the three sites, especially over Beijing. The biases of each wavelength can reach 0.2. One of the possible reasons for the larger bias over Beijing may be that the VIIRS_ML algorithm is actually the Dark-Target aerosol retrieval algorithm and this algorithm may overestimate the AOD values over bright surfaces such as urban centers by 0.2 [9,37]. The large VIIRS_ML biases may also be related to the errors of aerosol model type (to be discussed in next section).
We also calculate the aerosol Angstrom Exponent (AE) between 488 nm and 672 nm (AE = log (AOD 488nm /AOD 672nm ) log (672nm/488nm) ⁄ ) to describe the AOD spectral dependence. AE is often used as an indicator of aerosol size distribution which is related to aerosol type: AE~0 corresponds to large particles; AE~2 corresponds to small particles. The average AE values in Figure 1 are calculated from each matchup data set. The AE values of the VIIRS_EDR and the VIIRS_ML show large differences when comparing to that from CE318. AE biases of the VIIRS_EDR are 0.06 in Beijing, 0.31 in XiangHe and 0.55 in Xinglong, while the biases for the VIIRS_ML are -0.42 over Beijing, -0.17 in XiangHe and 0.20 over Xinglong. The AE of the VIIRS_EDR is larger than that from CE318 over all the three sites, but the AE of the VIIRS_ML is often lower than CE318 except for over Xinglong site. These indicate that the aerosol size of aerosol model type used in the VIIRS_EDR algorithm is smaller while VIIRS_ML is larger except for Xinglong (the aerosol model type used in VIIRS AOD retrieval will be discussed in next section).

Results of aerosol model and the optical properties inter-comparison
The aerosol model types used in the VIIRS_EDR and the VIIRS_ML over the NCP are evaluated by comparing SSA values with those derived from the CE318 inversion. The comparisons of aerosol optical properties between the VIIRS_EDR and CE318 are also shown in this section to help explain the error from aerosol model used in VIIRS_EDR. Figure 2 shows the frequencies of each aerosol model type used in the VIIRS_EDR (hereafter called M_VIIRS_EDR) AOD retrieval at the Beijing, XiangHe and Xinglong sites during the evaluation period. For M_VIIRS_EDR, the dust and clean urban aerosol models are the two dominant model types used in NCP region. In Beijing, the M_VIIRS_EDR shows that dust and clean urban models account for more than 80% of the aerosols during the evaluation period. The frequency of the polluted urban model is less than 1%. However, Beijing is a mega city with a population of approximately 21 million and 5 million vehicles located in the heavy polluted NCP region. It is undisputed that polluted urban aerosol is the dominant aerosol [38]. Thus, the M_VIIRS_EDR is unsuitable at the Beijing site. Because XiangHe is located 50 km to the east of Beijing, the M_VIIRS_EDR for XiangHe shows similar results to Beijing. However, XiangHe also shows some differences from Beijing: dust and clean urban models decrease and other models increase. As for Xinglong station (regional back ground station), the M_VIIRS_EDR shows more polluted urban and low absorption smoke and less dust models than the Beijing and XiangHe stations. From Beijing to Xinglong, the M_VIIRS_EDR shows that the frequency of polluted urban aerosol models increase, which is inconsistent with the fact that pollution decreases from Beijing to Xinglong according to past study results in the NCP region [13,28,39]. To show the aerosol model type differences between the VIIRS retrievals and CE318 sunphotometer observations, the seasonal values of SSA670nm of the CE318 inversion, VIIRS_EDR and VIIRS_ML in the NCP three sites are shown in Table 5. Comparing the VIIRS_EDR and CE318, the SSA670nm values from the VIIRS_EDR are higher than those from the CE318 during almost all seasons and over all the three sites. This result reflects more frequent weakly absorbing aerosol model type used in VIIRS_EDR retrievals in the NCP sites, as Figure 2 shows more frequency of dust and clean urban aerosol models. The largest difference between the VIIRS_EDR and CE318 at Beijing and XiangHe sites is shown during winter with 0.07 at Beijing and 0.08 at XiangHe. However, the difference at Xinglong during winter is smallest. The largest difference at the Xinglong site is occurred during the spring and summer and the difference value is 0.03, which is less than those at Beijing and XiangHe. These indicate that the M_VIIRS_EDR at Beijing and XiangHe have more errors than at Xinglong. To know how many aerosol types used in VIIRS_EDR are suitable, we calculated the accuracy rate of the M_VIIRS_EDR at the NCP three sites and the results are shown in Figure 3.

Evaluation of the VIIRS_EDR aerosol model type
The accuracy rates of the M_VIIRS_EDR over Beijing, XiangHe, and Xinglong are 0.24, 0.20, and 0.37, respectively (Figure 3a). The accuracy rate of the M_VIIRS_EDR over Xinglong is highest, which is consist with the lowest SSA difference between the VIIRS_EDR and CE318 at Xinglong. As for each model type, although the dust model is more used in the M_VIIRS_EDR (Figure 2), the accuracy rate of the dust model is very small at Beijing and XiangHe and even equal to zero at Xinglong site (Figure 3a and c). The accuracy rate of low absorption smoke is relatively higher than other models of M_VIIRS_EDR. The accuracy rate of the polluted urban aerosol model is practically zero because the frequency of polluted urban model occurred in the M_VIIRS_EDR over Beijing and XiangHe stations is very low (Figure 2) and it is different from the aerosol type of the CE318 of the matchup. All these results indicate that the M_VIIRS_EDR selected more weakly absorbing aerosol models in the NCP sites. It is worth to noting that with more frequencies of clean and dust aerosol models (higher SSA values) used in the VIIRS_EDR retrieval, the VIIRS_EDR AOD should have an underestimation if the surface reflectance characterization in VIIRS_EDR algorithm is perfect. The fact that the VIIRS_EDR AOD in Beijing has a high bias reflects that other factors (e.g. surface reflectance characterization or cloud contamination) are important error sources in the VIIRS_EDR algorithm, and their effects on aerosol retrievals override the effects from non-ideality in aerosol model types. This can be an interesting topic for future studies.

Evaluation of the VIIRS_ML aerosol model type
Using the same way of evaluating the M_VIIRS_EDR, the aerosol model type assumed in the VIIRS_ML (M_VIIRS_ML) in the NCP sites is evaluated and shown in Table 5 and Figure 3. From Table 5, it can be found that the biases of SSA670nm between the VIIRS_ML and CE318 are ≤ 0.04 in almost all sites and all seasons except for Xinglong during the winter. The average bias for all seasons is -0.015, 0.0075 and -0.0075 at Beijing, XiangHe and Xinglong, respectively. These biases are less than the differences between the VIIRS_EDR and CE318 at the corresponding site. That is likely because the M_VIIRS_ML is defined according to the AERONET sunphotometer inversion in local regions while M_VIIRS_EDR is from the predefined aerosol model types at five sites located in other places. Notably, the obvious undervaluation of the VIIRS_ML SSA over Beijing during almost all seasons may cause the overestimation of the VIIRS_ML AOD over Beijing (MB ≥ 0.26 shown in Table 2) versus the VIIRS_EDR.
However, there are some large differences of SSA670nm between the VIIRS_ML and CE318 in some seasons. The moderately absorbing model in spring over Beijing used in the VIIRS_ML may be inappropriate because the CE318 inversion shows weakly absorbing type. The largest difference is found at Xinglong (regional background site). The absolute biases in all seasons over this site are ≥ 0.3 and the highest bias 0.6 is occurred in the winter. The CE318 shows weakly absorbing aerosols in the winter that is in contrast to the moderately absorbing aerosols used by the VIIRS_ML. This may indicate that the M_VIIRS_ML over the NCP may be unsuitable for regional background sites.
In Figure 3, it can be found that the accuracy rates of the M_VIIRS_ML over Beijing, XiangHe, and Xinglong are 0.47, 0.51, and 0.46, respectively (Figure 3b). The average accuracy rate of the M_VIIRS_ML over the NCP region (0.48) is higher than that of M_VIIRS_EDR (0.27). As for each aerosol model type (Figure 3d), the accuracy rate of the weakly absorbing fine model is higher than the moderately absorbing fine model in Beijing and XiangHe. However, in Xinglong, the accuracy rate of the weakly absorbing fine model is lower than the moderately absorbing fine model. The large difference of the accuracy rate of the two model types in Beijing reflects that the moderately absorbing fine model assumed in the spring and winter may require careful consideration because of the dust in the spring and more strongly absorbing aerosols in the winter cannot be neglected [39].
For different sites, the accuracy rates of the M_VIIRS_EDR over Beijing and XiangHe are higher than their corresponding values of the M_VIIRS_EDR. From Beijing to Xinglong, the accuracy rate of the M_VIIRS_EDR decreases, while M_VIIRS_ML varies little with lowest value over Xinglong site. This reflects that the M_VIIRS_EDR is unsuitable over the NCP urban and suburban sites, while the M_VIIRS_ML is suitable over urban and suburban sites. The highest accuracy rate of the M_VIIRS_EDR over the NCP sites is only 0.37. Thus, the aerosol model type selection in the VIIRS_EDR algorithm is inappropriate in the NCP region, which may cause an important error in the AOD inversion.

Inter-comparison of aerosol properties between CE318 and the VIIRS_EDR
Since the M_VIIRS_EDR performs less well in the NCP region based on above analysis and M_VIIRS_ML is defined according to the CE318 inversion in local regions, we only compare the aerosol properties between CE318 and the VIIRS_EDR. Table 6 shows the averages of the SSA440nm of fine, coarse and bi-modal aerosols derived from CE318 (AOD440nm>0.4) and the VIIRS_EDR (AOD550nm>0.25) at the three sites over the time period of evaluation (2 May 2012-31 March 2014. The values of SSA440nm from the CE318 are less than that of the VIIRS_EDR for all modes and all sites, which indicates that VIIRS_EDR overestimates SSA440nm values over all sites (consistent with SSA670nm in Table 5). The difference of SSA440nm in the coarse mode (i.e., 0.09 in Beijing, 0.12 in XiangHe, and 0.06 in Xinglong) is larger than SSA440nm in fine mode (i.e., 0.02 in Beijing, 0.04 in XiangHe and 0.03 in Xinglong). This indicates that aerosol properties in the coarse mode in the VIIRS_EDR need to be revised for the NCP region. The biases of bi-mode SSA440nm are 0.03 (Beijing), 0.05 (XiangHe) and 0.03 (Xinglong). The overestimation of SSA but largest positive AOD MB of the VIIRS_EDR over Beijing site indicate again that surface reflectance positive bias overpowers the negative bias due to SSA of aerosol model type over Beijing. The largest overestimation of SSA and negative MB of AOD (in Table 3) are occurred at XiangHe, which may indicate that errors from the aerosol model type are overpowering at the XiangHe site. Since the SSA is calculated by inputting the aerosol parameters in Table 2 from Jackson et al. [6] to the Mie scattering code and large differences of SSA between the VIIRS_EDR and CE318 are found above, it is necessary to compare the aerosol properties between the VIIRS_EDR and CE318. The constrains of CE318 AOD440nm>0.4 and VIIRS_EDR AOD550nm>0.25 are also used and the daily aerosol optical properties retrieved from CE318 sky radiance measurements are averaged before the followed analysis. Table 7 shows the average aerosol physical properties from CE318 and the VIIRS_EDR. The refractive index is RI with real part is RI(r), while imaginary part is RI(i). Volume mean radius, standard deviation and volume concentration are r, σ and V, respectively. The fine and coarse mode aerosols are shown by f and c in bracket pairs. Figure 4 shows the seasonal comparison of normalized aerosol physical properties from CE318 and the VIIRS_EDR at the three sites. Each parameter is normalized between [0.1 1] to well show the difference of CE318 and the VIIRS_EDR. The length of each radius in the circle equals 1 and each radius direction stands for one aerosol parameter.  Distinctly differences of various parameters between CE318 and the VIIRS_EDR can be found. For the total averages in three sites ( Table 7), CE318 shows high values of RI(r) 1.47-1.49 and RI(i) 0.0083-0.013, which indicates that aerosol in the NCP is more absorptive than that used in the VIIRS_EDR, which corresponds the overestimation of SSA for the VIIRS_EDR ( Table 5).
The r(f) values from CE318 (0.19-0.22) are slightly higher than the values of the VIIRS_EDR (0. 18-0.19), while the r(c) values of CE318 (2.70-2.84) are significantly lower than the values of the VIIRS_EDR (3.15-3.51), which may cause more dust aerosol for the VIIRS_EDR (Figure 2) and also reflects the dust aerosol model in Cape Verde maybe different from the dust in Asia [40,41]. The σ(f) in CE318 is higher than the VIIRS_EDR but lower for σ(c). The V(f) values of CE318 are lower than that of the VIIRS_EDR except for Beijing station. The V(c) of the VIIRS_EDR is comparable with that of CE318 except in XiangHe, which consists with the lowest accuracy rate of the M_VIIRS_EDR in XiangHe (Figure 3a).
In Figure 4, the shapes generated by 8 parameters (1 -RI(r), 2 -RI(i), 3 -r(f), 4 -σ(f), 5 -V(f), 6 -r(c), 7 -σ(c), and 8 -V(c)) for CE318 and the VIIRS_EDR differ from each other and each of them (CE318 and VIIRS_EDR) shows a different seasonal variation. For Beijing, CE318 shows similar shapes in summer and autumn but significantly different in spring and winter; higher RI(r) in spring and winter and higher RI(i) in winter. While the VIIRS_EDR shows two pairs of similar shapes; spring-summer and autumn-winter. As for XiangHe station, CE318 shapes are similar to Beijing's but the VIIRS_EDR shows some difference from Beijing's; autumn is not similar to winter but similar to spring and summer. As for Xinglong, CE318 shapes show a distinct difference in winter compared to XiangHe's due to the lower RI(r), RI(i) and r(c) and higher r(f), while the VIIRS_EDR shapes are similar to XiangHe's in all seasons. The significant differences of the aerosol microphysical properties between the VIIRS_EDR and CE318 over the NCP indicate that the aerosol microphysical properties in the VIIRS_EDR algorithm are not suitable over the NCP region. Furthermore, this also reflects that the five models based on the five AERONET stations in the reference [24] can not be applied globally.

Conclusions
Using CE318 data at three sites in the NCP region to evaluate the VIIRS AOD, aerosol model types and aerosol optical properties used in the VIIRS_EDR and VIIRS_ML algorithms at three sites over the NCP, we conclude: a. The VIIRS_EDR AOD at 550 nm has a positive MB of 0.04-0.06 with R of 0.83-0.86. Among the three sites, the bias at Beijing is largest with MB of 0.10-0.15, RMSE of 0.27-0.30 and low %EE of 36.8%-38.8%. The quality flags analysis shows that using the high quality products of AOD at XiangHe and Xinglong are recommended but not at Beijing site. The VIIRS_ML AOD overestimates more over the NCP region with higher positive MB of 0.13-0.14 but shows higher correlation (0.93-0.94) with ground-based AOD. The results of evaluation of the VIIRS_ML for each site and quality flags analysis are found to be similar to that for the VIIRS_EDR. b. The aerosol model types used in the VIIRS_EDR AOD retrieval in the three sites are mostly dust and clean urban aerosol models (the frequencies of these two models in three sites are all larger than 50%) with less frequency of polluted urban aerosol models used (less than 1%). The accuracy rates of the M_VIIRS_EDR over the Beijing, XiangHe and Xinglong sites (0.24, 0.20 and 0.37) are lower than that of the M_VIIRS_ML (0.47, 0.51 and 0.46) during the evaluation period. c. The values of SSA440nm from CE318 are less than the VIIRS_EDR for all modes and sites, with differences of 0.02-0.04 for fine mode, 0.06-0.12 for coarse mode and 0.03-0.05 for bi-modes. The overestimation of SSA but positive AOD mean bias of the VIIRS_EDR indicate that other factors (e.g. surface reflectance characterization or cloud contamination) are important error sources in the VIIRS_EDR algorithm, and their effects on aerosol retrievals override the effects from non-ideality in aerosol model types. The differences of SSA670nm between VIIRS_ML and CE318 in the NCP are mostly less than 0.015 but high seasonal differences are also found. The undervaluation of SSA used in the VIIRS_ML algorithm over the NCP causes the overestimation of AOD, especially at Beijing site.
We recommend that the aerosol model types and the microphysical properties in the VIIRS_EDR algorithm in NCP region are not representative and need to be refined. The AOD bias in Beijing is largest but we do not find the lowest accuracy rate of the M_VIIRS_EDR. In addition, the higher values of the VIIRS_EDR SSA versus CE318, which should lead to lower AODs from the satellite inversion, are inconsistent with the positive MBs of AOD in the NCP region. All these points indicate that there are other error sources that need to be examined in the AOD retrieval for the VIIRS_EDR algorithm, especially for Beijing site. Future studies should investigate these potential sources of error including the surface reflectance and cloud contamination.