Interference of Heavy Aerosol Loading on the VIIRS Aerosol Optical Depth ( AOD ) Retrieval Algorithm

Aerosol optical depth (AOD) has been widely used in climate research, atmospheric environmental observations, and other applications. However, high AOD retrieval remains challenging over heavily polluted regions, such as the North China Plain (NCP). The Visible Infrared Imaging Radiometer Suite (VIIRS), which was designed as a successor to the Moderate Resolution Imaging Spectroradiometer (MODIS), will undertake the aerosol observations mission in the coming years. Using the VIIRS AOD retrieval algorithm as an example, we analyzed the influence of heavy aerosol loading through the 6SV radiative transfer model (RTM) with a focus on three aspects: cloud masking, ephemeral water body tests, and data quality estimation. First, certain pixels were mistakenly screened out as clouds and ephemeral water bodies because of heavy aerosols, resulting in the loss of AOD retrievals. Second, the greenness of the surface could not be accurately identified by the top of atmosphere (TOA) index, and the quality of the aggregation data may be artificially high. Thus, the AOD retrieval algorithm did not perform satisfactorily, indicated by the low availability of data coverage (at least 37.97% of all data records were missing according to ground-based observations) and overestimation of the data quality (high-quality data increased from 63.42% to 80.97% according to radiative simulations). To resolve these problems, the implementation of a spatial variability cloud mask method and surficial index are suggested in order to improve the algorithm.


Introduction
Atmospheric aerosols are solid and liquid particles that are suspended in the air and are often related to dust, smoke, soot, and sea salt.Climate models indicate that aerosols can significantly impact the radiation budget of the Earth [1], cloud formation [2], and precipitation [3].However, the uncertainty associated with the average climate impacts of aerosols remains large [4,5].Furthermore, aerosols can also impact human health in heavily polluted regions [6][7][8].The aerosol optical depth (AOD) is a basic optical property of aerosol research and has been broadly applied in climate research and atmospheric environmental observations.Satellite remote sensing has the advantage of observing and quantifying aerosol systems at a global scale from space [9].The Visible Infrared Imaging Radiometer Suite (VIIRS) aboard the Suomi National Polar-orbiting Partnership (Suomi-NPP) spacecraft was launched in October of 2011 [10].This instrument was largely built on the success of the Moderate Resolution Imaging Spectroradiometer (MODIS), which has successfully retrieved AOD for more than 15 years [11].The VIIRS was designed to have many similar features as its predecessors, and its aerosol algorithm was also based on the MODIS Dark-Target algorithm [12].VIIRS will perform tasks for climate and air quality applications after MODIS completes its mission.
The accuracy and availability of AOD products under polluted atmospheric environments, such as heavy polluted areas in China, are limited.Increasing fossil fuel consumption and biomass burning in China [13] have caused severe air pollution events and have worsened the atmospheric environment in northern China [14,15].Haze is an atmospheric phenomenon in which aerosol particles obscure the clarity of the sky and decrease the visibility below 10 km.Frequent haze events can be detected by ground-based observations [16,17] such as the AErosol Robotic NETwork (AERONET) [18], Chinese Sun Haze-meter Network (CSHNET) [19,20], and others.In situ observations have shown that the haze frequency and affected area have significantly increased over recent decades [21,22].However, limited in situ observations and uneven distributions could introduce considerable uncertainty.Satellite observations can provide wide spatial coverage and long-term data records.From a temporal perspective, Zhang et al. [23] used the Absorbing Aerosol Index (AAI) to show that the haze over northern and eastern China follows an increasing trend that is similar to the pattern that is observed from MODIS AOD [24,25].Climate Data Records (CDR) and regular air quality observations require an accurate, consistent, and wide-coverage AOD product [26,27].However, certain retrieval vacancies exist over areas of heavy aerosol loading, so the VIIRS AOD products are not acceptable under hazy conditions.Previous research has attempted to improve the ability to retrieve hazy AODs with MODIS data [11,28]; however, the quality of these products remains insufficient under polluted conditions.
AOD products may be influenced by unsuitable hypotheses for cloud masks and pixel selection and poor data quality assurance.The success of aerosol retrieval depends on the ability to screen out unsuitable pixels.The most important step is accurate cloud masking.The standard MODIS cloud mask (MxD35) is considered too cloud conservative and not clear-sky sufficient for aerosol retrieval [29,30].Therefore, after applying the MODIS Collection_4 algorithm, Martins et al. (2002) developed a new independent cloud mask that was mainly based on a spatial variability test.With VIIRS, the cloud mask in aerosol retrieval is based on a VIIRS cloud mask product (VCM), which is similar to MxD35 [12,31].Although the VCM product performs well as evaluated by MODIS and CALIPSO data [32], this product still has flaws when used in AOD retrieval.Furthermore, the ephemeral water body test calculates the top of the atmosphere (TOA) normalized difference vegetation index (NDVI) and excludes pixels below a certain threshold [12].However, heavy aerosol loading likely affects the calculation of this parameter.The data quality for the product's aggregation strategy depends on the pixel number and greenness.The greenness is defined by another TOA NDVI [33] that is minimally affected by the AOD.However, under hazy conditions, the hypothesis must be reconsidered.Therefore, under heavy aerosol loading conditions, AOD products can be affected by cloud masking, the ephemeral water body test, and quality assurance issues.
This article focuses on analyzing how these three factors influence the AOD retrieval algorithm under polluted atmospheric conditions and provides feasible advice.The data from the retrieval algorithm and analytical method are described in Section 2.Then, a radiative transfer simulation and certain examples are used to illustrate the results of a cloud mask, ephemeral water body test, and data aggregation in Section 3. Section 4 analyzes the causes of these impacts and provides a quantitative evaluation.Finally, Section 5 provides concluding remarks.

North China Plain
The study area is mainly located on the North China Plain (NCP), which is the largest alluvial plain in China.The NCP is surrounded by the Yanshan and Taihang mountains to the north and west and vast deserts in north-western China and Mongolia.Because its climate is characterized by both humid winds from the Pacific and dry winds from the interior of the Asian continent, the composition of the particulate matter in the air is complex and includes dust, sea salt, and industrial matter.The NCP is the most polluted area in China and one of the most polluted areas in the world because of weather conditions, terrain influences, and pollutant emissions (Figure 1).and vast deserts in north-western China and Mongolia.Because its climate is characterized by both humid winds from the Pacific and dry winds from the interior of the Asian continent, the composition of the particulate matter in the air is complex and includes dust, sea salt, and industrial matter.The NCP is the most polluted area in China and one of the most polluted areas in the world because of weather conditions, terrain influences, and pollutant emissions (Figure 1).During periods of intense atmospheric pollutant emissions and calm and steady weather, particulate concentrations are extremely high, resulting in hazy days and inhibiting AOD satellite retrievals.Figure 2 shows the VIIRS AOD product on three polluted days over the NCP.Certain areas with heavy aerosol loading, such as the invalid values in the red ellipses in Figure 2, lack satelliteretrieved AOD information.The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Ground-Based Observations
AERONET uses sky scanning spectral radiometers to make ground-based observations of atmospheric aerosol optical properties and precipitable water [18,34].The network provides a longterm and continuous dataset for satellite product validation.The level 2.0 AOD datasets have undergone cloud screening, calibration checks, and quality assurance.The level 2.0 dataset of Beijing_CAMS station from 2013 to 2016 was used in this study.
Because AERONET does not measure the 550 nm band, the AOD at 500 nm is used instead, During periods of intense atmospheric pollutant emissions and calm and steady weather, particulate concentrations are extremely high, resulting in hazy days and inhibiting AOD satellite retrievals.Figure 2 shows the VIIRS AOD product on three polluted days over the NCP.Certain areas with heavy aerosol loading, such as the invalid values in the red ellipses in Figure 2, lack satellite-retrieved AOD information.
Remote Sens. 2017, 9, 397 3 of 15 and vast deserts in north-western China and Mongolia.Because its climate is characterized by both humid winds from the Pacific and dry winds from the interior of the Asian continent, the composition of the particulate matter in the air is complex and includes dust, sea salt, and industrial matter.The NCP is the most polluted area in China and one of the most polluted areas in the world because of weather conditions, terrain influences, and pollutant emissions (Figure 1).During periods of intense atmospheric pollutant emissions and calm and steady weather, particulate concentrations are extremely high, resulting in hazy days and inhibiting AOD satellite retrievals.Figure 2 shows the VIIRS AOD product on three polluted days over the NCP.Certain areas with heavy aerosol loading, such as the invalid values in the red ellipses in Figure 2, lack satelliteretrieved AOD information.The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Ground-Based Observations
AERONET uses sky scanning spectral radiometers to make ground-based observations of atmospheric aerosol optical properties and precipitable water [18,34].The network provides a longterm and continuous dataset for satellite product validation.The level 2.0 AOD datasets have undergone cloud screening, calibration checks, and quality assurance.The level 2.0 dataset of Beijing_CAMS station from 2013 to 2016 was used in this study.
Because AERONET does not measure the 550 nm band, the AOD at 500 nm is used instead, The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Ground-Based Observations
AERONET uses sky scanning spectral radiometers to make ground-based observations of atmospheric aerosol optical properties and precipitable water [18,34].The network provides a long-term and continuous dataset for satellite product validation.The level 2.0 AOD datasets have undergone cloud screening, calibration checks, and quality assurance.The level 2.0 dataset of Beijing_CAMS station from 2013 to 2016 was used in this study.
Because AERONET does not measure the 550 nm band, the AOD at 500 nm is used instead, denoted AOD AERONET .We defined AOD AERONET > 0.6 as polluted.The protocol requires at least six AERONET measurements within a 3-h period centered on the satellite overpass time.In total, 187 days were selected as polluted days.The VIIRS EDR data matchup requires retrievals within a 27.5-km-radius circle that is centered on the AERONET station [35].If at least 20% of the pixels of all the potential retrievals are found to be clouds, the record is defined as being affected by cloud.The ephemeral water body and over-range tests are performed identically to the cloud test.

Satellite Data
The VIIRS instrument aboard the Suomi-NPP spacecraft was launched in October of 2011 and was designed to have similar capabilities as MODIS.Suomi-NPP orbits with a similar equator crossing time as Aqua.VIIRS data from 2013 to 2016 were used in this study, including sensor TOA reflectance (ρ TOA ), cloud mask data, geolocation data, and AOD data.The TOA measurement data were level 1b Sensor Data Records (SDR), including moderate-resolution bands (M-bands) with a spatial resolution of 750 m and imagery bands (I-bands) with a spatial resolution of 375 m.The cloud mask data were pixel-level Intermediate Product (IP) data.The AOD data were 6-km-resolution level 2 Environmental Data Records (EDR) aggregated as 8 × 8 IP retrievals.All of these data were downloaded from the National Oceanic and Atmospheric Administration (NOAA)'s Comprehensive Large Array-data Stewardship System (CLASS) website (http://www.class.ncdc.noaa.gov/saa/products/welcome).
MODIS is a key sensor aboard the Terra and Aqua satellites, which were launched in 2000 and 2002, respectively.MOD09 is the eight-day surface reflectance (ρ S ) product of MODIS/AQUA.The dataset was fused to the monthly average surface reflectance in January 2015 by using a minimum method [36].These mature surficial reflectance data were used as input data in the radiative simulation.The dataset was downloaded from The Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) managed by the National Aeronautics and Space Administration (NASA) (https://ladsweb.nascom.nasa.gov/search/).

Radiative Transfer Simulation
The Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer model (RTM) provides accurate simulations of satellite and plane observations [37].The new vector version (6SV) of this code can work in both scalar and vector modes [38].
In this article, we used the 6SV RTM to simulate the radiative transfer procedure.The angle, aerosol type and target altitude were not important factors when analyzing the influences of high AOD values.Therefore, the satellite zenith angle, solar zenith angle, and relative azimuth angle were set to 30 • , 30 • and 60 • , respectively.The aerosol type was assumed to be continental, and the target altitude was set to 0.
In Section 3.1, ρ TOA was simulated under each AOD (in 550 nm) from 0 to 3 using different ρ s values for VIIRS bands M1 and M3.
In Section 3.2, three typical land covers-soil, vegetation, and water-are selected.The surface reflectance information for several land cover types is listed in Table 1.We used 6SV to calculate ρ TOA in the VIIRS bands I1 and I2 and then obtained the TOA NDVI.In Section 3.3, the MODIS surface reflectance product was used to simulate the TOA reflectance in the NCP area (113 • E-116 • E, 34 • N-39 • N) under three different atmospheric conditions.The air is assumed to be clean when AOD = 0.1, it is lightly polluted when AOD = 1, and it is heavily polluted when AOD = 2.

Cloud Mask Algorithm
The VIIRS cloud mask depends on an external identification result, the VCM-IP product, which is not robust enough for aerosol retrievals.The VCM technique incorporates several cloud detection tests to determine whether a pixel is obstructed by a cloud, and the VIIRS pixels are assigned a label depending on the cloud confidence level, i.e., confidently cloudy, probably cloudy, probably clear, or confidently clear [39,40].These tests include reflectance, brightness temperature (BT), brightness temperature difference (BTD), and spatial tests using M-band and I-band data.
The spatial variability in the reflectance at the TOA is suitable for a cloud mask that is devoted to the retrieval of aerosol data [41].The spatial test uses the absolute standard deviation of every 3 × 3 pixel (3 × 3 STD) threshold to identify clouds.The 3 × 3 STD (σ) is calculated as follows: where ρ i is the TOA reflectance of each pixel and ρ is the average TOA reflectance of all nine pixels.In this study, we calculated the 3 × 3 STD in bands M1 (0.412 µm) and M3 (0.486 µm).

Ephemeral Water Body Test Method
The presence of surface water over land can affect retrieval algorithms; thus, an ephemeral water body detection test was applied to overcome this deficiency.This test is based on the TOA NDVI values calculated using bands I1 (0.638 µm) and I2 (0.862 µm) with the following equation [12] where ρ I1 and ρ I2 are the TOA reflectances of bands I1 and I2, respectively.The TOA NDVI threshold is 0.1.If the TOA NDVI value of a pixel is less than the threshold, the pixel is identified as an ephemeral water body [12].

EDR Product Aggregation Strategy
The VIIRS AOD EDR product was constructed by aggregating 8 × 8 arrays of pixel levels retrieved from AOD IP.The overall quality (high, medium, and low) of the EDR depends on the pixel number, which is based on the IP quality within the EDR cell.
The quality of IP retrievals may be affected by the zenith and azimuth angles, clouds, pixel greenness, and other factors.Among these parameters, greenness was the only factor affected by aerosol loading.Thus, we only considered greenness and assumed that the other factors remained unchanged.The brightness index was used to identify when a pixel was dominated by a bright surface, less vegetated soil, or vegetation-dominated soil [12].This index is also referred to as NDVI SWIR [33] and is calculated by using the equation where ρ M8 and ρ M11 are the TOA reflectance values of bands M8 and M11, respectively.The VIIRS algorithm identifies a bright pixel when A vegetation-dominated pixel corresponds to NDV I SW IR > 0.2, and the other pixels with intermediate values are defined as less vegetated."Vegetation-dominated" is a necessary condition for "high" quality, and less vegetated conditions should be sufficient to assign the "degraded" quality flag [42].
The EDR product was labelled "high/medium/low" quality, depending on the number of AOD IP retrievals of different quality in the 8 × 8 EDR cell.The aggregation logic is displayed as a flowchart in Figure 3 according to the Algorithm Theoretical Basis Document (ATBD) [42].

Cloud Mask
The aerosol retrieval only works under cloud-free conditions.Thus, cloudy pixels must be identified and removed.The VIIRS AOD retrieval depends on the cloud mask when using the information from the VCM input.
We obtained cloud mask information in the form of "Confidently Cloudy" and "Probably Cloudy" from the AOD IP (IVAOT) quality flag (QF) data.To reveal inadequate cloud tests under heavy aerosol loading, we chose two locally hazy days-23 December 2013 and 18 March 2016which are shown in Figure 4. On 23 December 2013 (Figure 4a), the NCP was covered by heavy haze but was almost cloud free.However, the VCM result indicates clouds in the area.A similar indication was found on 18 March 2016 (Figure 4b).These pixels would be excluded, which would result in no AOD retrieval.The VCM performs well in clear areas, although heavy aerosol loading may mislead the tests.

Cloud Mask
The aerosol retrieval only works under cloud-free conditions.Thus, cloudy pixels must be identified and removed.The VIIRS AOD retrieval depends on the cloud mask when using the information from the VCM input.
We obtained cloud mask information in the form of "Confidently Cloudy" and "Probably Cloudy" from the AOD IP (IVAOT) quality flag (QF) data.To reveal inadequate cloud tests under heavy aerosol loading, we chose two locally hazy days-23 December 2013 and 18 March 2016-which are shown in Figure 4. On 23 December 2013 (Figure 4a), the NCP was covered by heavy haze but was almost cloud free.However, the VCM result indicates clouds in the area.A similar indication was found on 18 March 2016 (Figure 4b).These pixels would be excluded, which would result in no AOD retrieval.The VCM performs well in clear areas, although heavy aerosol loading may mislead the tests.
heavy aerosol loading, we chose two locally hazy days-23 December 2013 and 18 March 2016which are shown in Figure 4. On 23 December 2013 (Figure 4a), the NCP was covered by heavy haze but was almost cloud free.However, the VCM result indicates clouds in the area.A similar indication was found on 18 March 2016 (Figure 4b).These pixels would be excluded, which would result in no AOD retrieval.The VCM performs well in clear areas, although heavy aerosol loading may mislead the tests.In the visible channel, aerosols show a highly homogeneous spatial structure that can be easily separated from most clouds.Thus, the spatial variability test is efficient at masking clouds during aerosol retrieval [41].We used a 6SV radiative transfer core to simulate ρTOA in the VIIRS bands M1 In the visible channel, aerosols show a highly homogeneous spatial structure that can be easily separated from most clouds.Thus, the spatial variability test is efficient at masking clouds during aerosol retrieval [41].We used a 6SV radiative transfer core to simulate ρ TOA in the VIIRS bands M1 and M3.As shown in Figure 5, the difference in ρ TOA values is smaller under high AOD conditions than under low AOD conditions.Therefore, the ρ TOA values under heavy aerosol conditions are more homogeneous than those under clear sky conditions.The high degree of homogeneity under heavy aerosol conditions makes these pixels easier to distinguish from clouds.
Remote Sens. 2017, 9, 397 7 of 15 and M3.As shown in Figure 5, the difference in ρTOA values is smaller under high AOD conditions than under low AOD conditions.Therefore, the ρTOA values under heavy aerosol conditions are more homogeneous than those under clear sky conditions.The high degree of homogeneity under heavy aerosol conditions makes these pixels easier to distinguish from clouds.We pre-selected more than 200,000 pixels from the RGB image in eastern Asia and classified them into three groups-clouds, haze, and clear sky-by visual interpretation.Figure 6 shows the statistical STD histograms of every 3 × 3 set of pixels in the VIIRS bands M1 and M3.The upper row is the frequency distribution diagram, and the bottom row is the cumulative frequency distribution diagram.In Figure 6a,b, the frequency peak of haze was located farther from that of clouds than that of clear sky.In Figure 6d, the differences among clouds, haze, and clear sky are more significant.Therefore, clouds and haze are easily separated in standard deviation histograms.Based on the histogram in Figure 6, the thresholds were defined as the separator between clouds and clear sky.We hope to reserve clear sky pixels as much as possible on the basis of the majority of cloud pixels being screened out.The red vertical lines are the suggested thresholds based on the 3 × 3 STD test for VIIRS bands M1 and M3.The threshold in M1 is 0.005, and the threshold in M3 is 0.01.These thresholds generally do not exclude aerosol pixels (less than 2% of these samples) and only allow a small amount little cloud contamination (less than 5%).We pre-selected more than 200,000 pixels from the RGB image in eastern Asia and classified them into three groups-clouds, haze, and clear sky-by visual interpretation.Figure 6 shows the statistical STD histograms of every 3 × 3 set of pixels in the VIIRS bands M1 and M3.The upper row is the frequency distribution diagram, and the bottom row is the cumulative frequency distribution diagram.In Figure 6a,b, the frequency peak of haze was located farther from that of clouds than that of clear sky.In Figure 6d, the differences among clouds, haze, and clear sky are more significant.Therefore, clouds and haze are easily separated in standard deviation histograms.Based on the histogram in Figure 6, the thresholds were defined as the separator between clouds and clear sky.We hope to reserve clear sky pixels as much as possible on the basis of the majority of cloud pixels being screened out.The red vertical lines are the suggested thresholds based on the 3 × 3 STD test for VIIRS bands M1 and M3.The threshold in M1 is 0.005, and the threshold in M3 is 0.01.These thresholds generally do not exclude aerosol pixels (less than 2% of these samples) and only allow a small amount little cloud contamination (less than 5%).
is the frequency distribution diagram, and the bottom row is the cumulative frequency distribution diagram.In Figure 6a,b, the frequency peak of haze was located farther from that of clouds than that of clear sky.In Figure 6d, the differences among clouds, haze, and clear sky are more significant.Therefore, clouds and haze are easily separated in standard deviation histograms.Based on the histogram in Figure 6, the thresholds were defined as the separator between clouds and clear sky.We hope to reserve clear sky pixels as much as possible on the basis of the majority of cloud pixels being screened out.The red vertical lines are the suggested thresholds based on the 3 × 3 STD test for VIIRS bands M1 and M3.The threshold in M1 is 0.005, and the threshold in M3 is 0.01.These thresholds generally do not exclude aerosol pixels (less than 2% of these samples) and only allow a small amount little cloud contamination (less than 5%).

Ephemeral Water Body Test
We extracted the results of the ephemeral water body test from the "IVAOT" QF data from 13 January 2014 and 10 March 2014, as shown in Figure 7.In these days of low precipitation, the NCP area could not contain as large a range of ephemeral water bodies as the algorithm test.These identification errors only occur in heavy aerosol loading areas.
Remote Sens. 2017, 9, 397 8 of 15 We extracted the results of the ephemeral water body test from the "IVAOT" QF data from 13 January 2014 and 10 March 2014, as shown in Figure 7.In these days of low precipitation, the NCP area could not contain as large a range of ephemeral water bodies as the algorithm test.These identification errors only occur in heavy aerosol loading areas.Figure 8 shows the TOA NDVI simulation results for different land cover types.As seen in the vegetation line, the NDVI decreases with increasing AOD, and the NDVI value is always greater than 0. The TOA NDVI of water exhibits an almost constant negative pattern at any AOD range.Additionally, the TOA NDVI value of soil is obviously influenced by AOD: the values are greater than zero at low AOD values and close to or even less than zero at higher AOD values.Figure 8 shows the TOA NDVI simulation results for different land cover types.As seen in the vegetation line, the NDVI decreases with increasing AOD, and the NDVI value is always greater than 0. The TOA NDVI of water exhibits an almost constant negative pattern at any AOD range.Additionally, the TOA NDVI value of soil is obviously influenced by AOD: the values are greater than zero at low AOD values and close to or even less than zero at higher AOD values.Figure 8 shows the TOA NDVI simulation results for different land cover types.As seen in the vegetation line, the NDVI decreases with increasing AOD, and the NDVI value is always greater than 0. The TOA NDVI of water exhibits an almost constant negative pattern at any AOD range.Additionally, the TOA NDVI value of soil is obviously influenced by AOD: the values are greater than zero at low AOD values and close to or even less than zero at higher AOD values.
Figure 8. TOA NDVI simulation results for six types of land cover.The satellite zenith angle was 30°, the solar zenith angle was 30°, and the relative azimuth angle was 120°.In this simulation, the aerosol type was assumed to be continental, and the AOD ranged from 0 to 3.
The threshold of the VIIRS ephemeral water body test is 0.1.If the TOA NDVI value of a pixel is less than the threshold, then the pixel is identified as ephemeral water body.As the simulation in Figure 8 shows, certain soil pixels tended to be identified as ephemeral water bodies at high AOD levels., the solar zenith angle was 30 • , and the relative azimuth angle was 120 • .In this simulation, the aerosol type was assumed to be continental, and the AOD ranged from 0 to 3.
The threshold of the VIIRS ephemeral water body test is 0.1.If the TOA NDVI value of a pixel is less than the threshold, then the pixel is identified as ephemeral water body.As the simulation in Figure 8 shows, certain soil pixels tended to be identified as ephemeral water bodies at high AOD levels.

Available Retrievals
Based on the NOAA AOD distribution image (Figure 2), the retrievals failed in certain areas with heavy aerosol conditions.The cause of this failure was demonstrated by radiative simulation.Here, we use AERONET ground-based observation data from 2013 to 2016 to quantify this deficiency.The high AOD values were not retrieved because of the mistaken identification as clouds or ephemeral water bodies.Moreover, this process might have occurred because the AOD retrieval range (0-2) was exceeded.Table 2 lists the number of retrieval results compared with the AERONET dataset.Clouds, ephemeral water bodies, and the retrieval ranges are denoted as a, b, and c, respectively.Of the 187 AERONET observations, only 67 days have complete retrievals, and 71 days have no retrievals at all.The no-retrieval days comprised 37.97% of the dataset.Additionally, 49 days (26.20%) feature partial retrievals.The factors that were responsible for the lack of retrieval varied and included clouds, ephemeral water bodies, retrieval ranges, or combinations of multiple factors.In the partial retrieval group, the area with the 27.5-km-radius circle centered on the Beijing_CAMS station contains at least 20% pixels with retrieval data and at least 20% pixels with no retrieval data.We cannot confirm whether the clouds were real or whether pixels were mistakenly identified as clouds.However, some misidentification had to be included for clarity.In this example, the interference percentage of clouds and ephemeral water bodies was approximately equal.In summary, at least 37.97% of the high AOD data were not retrieved at the Beijing_CAMS station because of heavy aerosol loading.

Quality Assurance
The reason why NDVI SWIR was selected as a quality indicator was that the radiance in these bands is only slightly influenced by aerosol loading [33].However, this hypothesis is inappropriate when the AOD is high.
Figure 9a shows the MODIS 1.23 µm surface reflectance of the NCP in January 2015.Then, we used the reflectance to calculate NDVI SWIR when AOD = 0.1, 1 and 2, which are denoted NDVI SWIR _0.1, NDVI SWIR _1 and NDVI SWIR _2, respectively.Figure 9b is the difference between NDVI SWIR values under AOD = 0.1 and AOD = 1 (NDVI SWIR _1-NDVI SWIR _0.1). Figure 9c shows the difference between the NDVI SWIR values under AOD = 0.1 and AOD = 2 (NDVI SWIR _2-NDVI SWIR _0.1).As shown in Figure 9b,c, the NDVI SWIR values of aerosol loads of AOD = 1 and 2 were much higher than those of AOD = 0.1.These differences can also be observed in Figure 10, which is a histogram of NDVI SWIR values under different aerosol loading conditions.The frequency peak moves towards higher NDVI SWIR values as the AOD increases.The pixels with higher NDVI SWIR values probably tended to be identified as vegetation-dominated pixels.Although longer wavelengths are less influenced by aerosols, the influence cannot be neglected when the AOD is high.Therefore, the TOA NDVI SWIR is not suitable for land cover detection at high AOD levels.Every pixel was classified as "less vegetated" or "vegetated-dominated" according to the conditions in Section 2.5.3.Figure 9d-f   The EDR data QF depends on the quality and quantity of the IP data.Because this work was only concerned with surficial dominant types, we calculated the QF with 6-km-resolution EDR cells (including 8 × 8 pixels) using the aggregation strategy referred to in Section 2.5.3.Figure 11 shows the EDR QF when AOD = 0.1, 1, and 2. In Figure 11a, the data quality is better because the AOD is low.However, in Figure 11b,c, more pixels are identified as vegetation-dominated areas (see Figure 9e,f) when the AOD is higher, resulting in higher-quality EDR data.Certain EDR cells were identified as high quality in the central and western areas at AOD = 2, as shown in Figure 11c, whereas the cells were identified as medium quality at AOD = 0.1.The EDR data QF depends on the quality and quantity of the IP data.Because this work was only concerned with surficial dominant types, we calculated the QF with 6-km-resolution EDR cells (including 8 × 8 pixels) using the aggregation strategy referred to in Section 2.5.3.Figure 11 shows the EDR QF when AOD = 0.1, 1, and 2. In Figure 11a, the data quality is better because the AOD is low.However, in Figure 11b,c, more pixels are identified as vegetation-dominated areas (see Figure 9e,f) when the AOD is higher, resulting in higher-quality EDR data.Certain EDR cells were identified as high quality in the central and western areas at AOD = 2, as shown in Figure 11c, whereas the cells were identified as medium quality at AOD = 0.1.Table 3 lists the numbers and percentages of high-and medium-quality pixels under different aerosol loads, which correspond to the data in Figure 11.The high-quality EDR data percentage increased from 63.42% at AOD value of 0.1 to 80.91% at AOD value of 2.  Table 3 lists the numbers and percentages of high-and medium-quality pixels under different aerosol loads, which correspond to the data in Figure 11.The high-quality EDR data percentage increased from 63.42% at AOD value of 0.1 to 80.91% at AOD value of 2.

Analysis and Discussion
The AOD retrieval algorithm has three key basic scientific problems: cloud masks, aerosol models, and surface reflectance.Furthermore, the algorithm faces complex issues under high aerosol loading conditions.We analyzed and investigated these influences from the perspective of cloud masks, ephemeral water body tests, and data quality.These main effects influenced the retrieval availability and data quality.

Impact on Retrieval Availability
Accurate pixel selection is an important step before AOD inversions.Inappropriate pixels, such as clouds, bright surfaces, ice, snow, sun glint, among others, must be screened out.The cloud mask represents the most uncertain part of the AOD retrieval.Certain cloud mask methods [43][44][45] may have difficulty differentiating between clouds and high aerosol concentrations [41].The VCM identification results also exhibit this limitation under high aerosol loading conditions, such as those shown in Figure 4.This type of error diminishes the retrieval area because of the masking of certain hazy pixels as clouds.
Identifying ephemeral water bodies is also an important component of pixel selection.However, when aerosol loading is high, the TOA radiative characteristics are obviously affected by the atmosphere and result in inaccurate identifications.As shown in Figure 8, the "soil 2" and "soil 3" lines decreased to values close to 0 when the AOD is high (AOD > 2.5).For less vegetated NCP areas, the NDVI of bare soil may decrease close to or less than zero under the influence of high AOD values.Therefore, these pixels would be screened out as ephemeral water body pixels.In certain sub-tropical and temperate climate regions, these mistaken identifications will occur under heavy aerosol loading when the vegetation coverage decreases in winter.
The influence of these two issues discussed above is focused on the available amount of AOD retrieval data.Cloud masks and ephemeral water body tests are influenced by heavy aerosol loading, resulting in the elimination of certain suitable pixels before retrieval.Therefore, AOD datasets tend to lack information in areas with high AOD values, resulting in underestimations in the spatial and temporal average AOD.This influence would lead to lower concentrations compared to ground-based observations and chemistry transport model (CTM) simulation results, which is very important for radiation and climate research as well as for air quality monitoring.

Impact on Data Quality
Using the Dark-Target retrieval algorithm, the AOD accuracy depends on surficial type and whether the type is consistent with the dark dense vegetation hypothesis.Therefore, the greenness parameter (TOA NDVI SWIR ) is an important indicator for determining the data quality.The NDVI SWIR is only slightly affected by the atmosphere at low AOD values, but the effects at high AOD values cannot be ignored.An investigation into the effects of different AOD values on data quality was performed for the NCP.The high-quality data percentage increased from 63.42% (AOD = 0.1) to 80.97% (AOD = 2).In this case, the data quality was different for each AOD assumption.However, in a realistic situation, the data quality should not change with increasing AOD.Furthermore, the data quality should be lower because of the lower accuracy of RTM simulations under high AOD conditions.Therefore, the quality degree is overestimated when the AOD is high.For scientific climate research that uses high-quality AOD data, overestimated data may introduce uncertainty errors that are related to high AOD values and may influence the uniformity of the data standard.
Aerosol models and radiative simulations are other influential factors.The aerosol model that is used in satellite observations is based on cluster analysis with global long-term AERONET data [46,47].However, the aerosol components are complex and change when pollution episodes occur in extremely polluted countries or areas.Therefore, aerosol models are unlikely to be truly representative of the optical conditions viewed by satellites.Furthermore, compared to the Monte Carlo code, the 6SV RTM performs worse when AOD = 0.8 than when AOD = 0.2 [48].Additionally, the error would be further magnified at higher AOD values because of multiple scattering.

Proposed Solution
Feasible methods are recommended for resolving these problems.On one hand, we should consider using a spatial variability test method that is effective for cloud tests with aerosol retrieval.According to the histogram analysis, the difference between clouds and haze is enhanced using 3 × 3 STDs.High aerosol concentration areas can be easily distinguished from clouds with this test.The suggested thresholds are 0.005 for the M1 band and 0.01 for the M3 band.On the other hand, the problems of ephemeral water bodies and data quality are both caused by atmospheric interference.Therefore, we recommend using the surficial parameters to overcome this issue.We could use a simplified retrieval algorithm to calculate the approximate values of AOD and use them as crude atmospheric corrections.Additionally, pre-calculated surface parameters could be used instead of TOA parameters to provide more accurate estimates of ephemeral water bodies and vegetation coverage and diminish the interference of high aerosol loading.Additionally, researchers should note that the quality of high AOD data is not credible when using satellite AOD products and should use these data carefully.

Conclusions
As a key optical and physical parameter of aerosols, AOD is critical for environment and climate research.However, under heavy aerosol loading conditions, cloud masks and ephemeral water body tests decrease the amount of available retrieval data, and the algorithm overestimates the data quality.Candidate pixels were identified as clouds or ephemeral water bodies because of heavy aerosol loading, indicating that certain crucial research areas that are associated with high AOD would not be retrieved.This retrieval coverage limitation can be attributed to incorrect cloud masks and water body tests.According to the statistical results from the AERONET Beijing_CAMS station, at least 37.97% of the high AOD data were not retrieved.Additionally, more high-quality (80.97% in radiative simulation) AOD data were retrieved under polluted atmospheric conditions because of improper TOA NDVI tests.These factors restrict the use of the AOD product.In environmental air quality assessments, a lack of high values may lead to underestimations in the spatial and temporal mean AOD values.This study used VIIRS AOD data as an example to qualitatively and semi-quantitatively analyze the limitations and deficiencies of AOD data under heavy aerosol loading.This work was necessary to improve algorithms and data applications.

Figure 1 .
Figure 1.Annual average AOD distributions over the research area in 2015 (MODIS Collection 6 Deep Blue AOD at 550 nm).The right-hand figure shows the NCP, which is marked by a black square frame in the left-hand figure.

Figure 2 .
Figure 2. National Oceanic and Atmospheric Administration (NOAA)'s VIIRS AOD products (all data quality) over hazy areas.The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Figure 1 .
Figure 1.Annual average AOD distributions over the research area in 2015 (MODIS Collection 6 Deep Blue AOD at 550 nm).The right-hand figure shows the NCP, which is marked by a black square frame in the left-hand figure.

Figure 1 .
Figure 1.Annual average AOD distributions over the research area in 2015 (MODIS Collection 6 Deep Blue AOD at 550 nm).The right-hand figure shows the NCP, which is marked by a black square frame in the left-hand figure.

Figure 2 .
Figure 2. National Oceanic and Atmospheric Administration (NOAA)'s VIIRS AOD products (all data quality) over hazy areas.The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Figure 2 .
Figure 2. National Oceanic and Atmospheric Administration (NOAA)'s VIIRS AOD products (all data quality) over hazy areas.The AOD products were overlaid on the true color image, and no retrieval areas were set as transparent.Some AOD values are invalid, which are marked with red ellipses, because of heavy haze events.

Figure 4 .
Figure 4. Two-day VIIRS true color images (a,b) and NOAA cloud mask result (c,d) over the NCP on 23 December 2013 and 18 March 2016.The cloud pixels are represented in blue in the cloud mask result.

Figure 4 .
Figure 4. Two-day VIIRS true color images (a,b) and NOAA cloud mask result (c,d) over the NCP on 23 December 2013 and 18 March 2016.The cloud pixels are represented in blue in the cloud mask result.

Figure 5 .
Figure 5. ρTOA simulation for the M1 (a) and M3 (b) bands under different AOD values (ranging from 0 to 3).The different lines represent several surface reflectance values.

Figure 5 .
Figure 5. ρ TOA simulation for the M1 (a) and M3 (b) bands under different AOD values (ranging from 0 to 3).The different lines represent several surface reflectance values.

Figure 6 .
Figure 6.Histograms (a,b) and cumulative histograms (c,d) of the ρTOA STD in the VIIRS bands M1 (left column) and M3 (right column) for three types of pixels, including clouds (blue), haze (grey), and clear sky (green).The red lines are the suggested thresholds of the spatial variability test, which are 0.005 for M1 and 0.01 for M3.

Figure 6 .
Figure 6.Histograms (a,b) and cumulative histograms (c,d) of the ρ TOA STD in the VIIRS bands M1 (left column) and M3 (right column) for three types of pixels, including clouds (blue), haze (grey), and clear sky (green).The red lines are the suggested thresholds of the spatial variability test, which are 0.005 for M1 and 0.01 for M3.

Figure 7 .
Figure 7. VIIRS true color image on 13 January 2014 (a) and 10 March 2014 (b) and the corresponding ephemeral water body test results (c,d) over the NCP.The ephemeral water body pixels are represented in blue.

Figure 7 .
Figure 7. VIIRS true color image on 13 January 2014 (a) and 10 March 2014 (b) and the corresponding ephemeral water body test results (c,d) over the NCP.The ephemeral water body pixels are represented in blue.

Figure 7 .
Figure 7. VIIRS true color image on 13 January 2014 (a) and 10 March 2014 (b) and the corresponding ephemeral water body test results (c,d) over the NCP.The ephemeral water body pixels are represented in blue.

Figure 8 .
Figure 8. TOA NDVI simulation results for six types of land cover.The satellite zenith angle was 30• , the solar zenith angle was 30 • , and the relative azimuth angle was 120 • .In this simulation, the aerosol type was assumed to be continental, and the AOD ranged from 0 to 3.

Figure 9 .
Figure 9. (a) MODIS surface reflectance at 1.23 µm over the NCP.The NDVISWIR values were simulated by using the surface reflectance under aerosol conditions of AOD = 0.1, 1 and 2. (b) Difference in the NDVISWIR simulation values between AOD = 1 and 0.1.(c) Difference in the NDVISWIR simulation values between AOD = 2 and 0.1.The surface type identification results under different aerosol loads of (d) AOD = 0.1, (e) AOD = 1 and (f) AOD = 2 were also identified.

Figure 9 .
Figure 9. (a) MODIS surface reflectance at 1.23 µm over the NCP.The NDVI SWIR values were simulated by using the surface reflectance under aerosol conditions of AOD = 0.1, 1 and 2. (b) Difference in the NDVI SWIR simulation values between AOD = 1 and 0.1.(c) Difference in the NDVI SWIR simulation values between AOD = 2 and 0.1.The surface type identification results under different aerosol loads of (d) AOD = 0.1, (e) AOD = 1 and (f) AOD = 2 were also identified.

Figure 9 .
Figure 9. (a) MODIS surface reflectance at 1.23 µm over the NCP.The NDVISWIR values were simulated by using the surface reflectance under aerosol conditions of AOD = 0.1, 1 and 2. (b) Difference in the NDVISWIR simulation values between AOD = 1 and 0.1.(c) Difference in the NDVISWIR simulation values between AOD = 2 and 0.1.The surface type identification results under different aerosol loads of (d) AOD = 0.1, (e) AOD = 1 and (f) AOD = 2 were also identified.

Figure 10 .
Figure 10.Histograms of NDVISWIR when AOD = 0.1, 1, and 2. The red lines represent the NDVISWIR frequency peak under the three atmospheric conditions.
show the classification results for AOD = 0.1, 1 and 2, respectively.The vegetation-dominated areas under AOD = 2 conditions are much larger than those under AOD = 0.1.

Figure 10 .
Figure 10.Histograms of NDVI SWIR when AOD = 0.1, 1, and 2. The red lines represent the NDVI SWIR frequency peak under the three different atmospheric conditions.

Table 1 .
Surface reflectance information of five types of land cover.

Table 2 .
Number of VIIRS AOD EDR retrievals compared with the AERONET Beijing_CAMS station.

Table 3 .
EDR data quality statistics under three aerosol loading conditions.

Table 3 .
EDR data quality statistics under three aerosol loading conditions.