Impact of Aerosol Vertical Distribution on Aerosol Optical Depth Retrieval from Passive Satellite Sensors

: When retrieving Aerosol Optical Depth (AOD) from passive satellite sensors, the vertical distribution of aerosols usually needs to be assumed, potentially causing uncertainties in the retrievals. In this study, we use the Moderate Resolution Spectroradiometer (MODIS) and Visible Infrared Imaging Radiometer Suite (VIIRS) sensors as examples to investigate the impact of aerosol vertical distribution on AOD retrievals. A series of sensitivity experiments was conducted using radiative transfer models with di ﬀ erent aerosol proﬁles and surface conditions. Assuming a 0.2 AOD, we found that the AOD retrieval error is the most sensitive to the vertical distribution of absorbing aerosols; a − 1 km error in aerosol scale height can lead to a ~30% AOD retrieval error. Moreover, for this aerosol type, ignoring the existence of the boundary layer can further result in a ~10% AOD retrieval error. The di ﬀ erences in the vertical distribution of scattering and absorbing aerosols within the same column may also cause − 15% (scattering aerosols above absorbing aerosols) to 15% (scattering aerosols below absorbing aerosols) errors. Surface reﬂectance also plays an important role in a ﬀ ecting the AOD retrieval error, with higher errors over brighter surfaces in general. The physical mechanism associated with the AOD retrieval errors is also discussed. Finally, by replacing the default exponential proﬁle with the observed aerosol vertical proﬁle by a micro-pulse lidar at the Beijing-PKU site in the VIIRS retrieval algorithm, the retrieved AOD shows a much better agreement with surface observations, with the correlation coe ﬃ cient increased from 0.63 to 0.83 and bias decreased from 0.15 to 0.03. Our study highlights the importance of aerosol vertical proﬁle assumption in satellite AOD retrievals, and indicates that considering more realistic proﬁles can help reduce the uncertainties.


Introduction
Aerosol observations are essential in climate and environmental studies such as quantifying Earth's energy budget and characterizing surface air quality. In recent years, satellite remote sensing has advanced the observations of aerosol properties on both local and global scales. For example, historical sensors like the Advanced Very High Resolution Radiometer [1] on board polar-orbiting meteorological satellites measure the Earth's reflectance in six channels. Modern sensors like the Moderate Resolution Spectroradiometer [2,3], launched on board Aqua in May 2002 and Terra in December 1999, provide aerosol information at up to seven spectral bands, ranging from visible to near-IR wavelengths, and cover the entire globe every 1 to 2 days. Another advanced satellite sensor -the Visible Infrared Imaging Radiometer Suites [4]-was launched onboard the Suomi National Polar-orbiting Partnership (S-NPP) satellite in October 2011. It collects images and measurements of land and atmosphere in 22 spectral bands, covering wavelengths from 412 to 12,050 nm, and provides routine retrievals of land, aerosol and cloud properties. Sophisticated algorithms have been developed to retrieve aerosol property, in particular, aerosol optical depth, from these various sensors. The overall accuracy of Moderate Resolution Spectroradiometer (MODIS)-retrieved Aerosol Optical Depth (AOD) is expected to be ±(0.05% + 15% × AOD) over land [3]. Although global validation proves that the current sensors largely meet these accuracy requirements, regionally, the uncertainties can still be large. For example, over China, only 50.6% of Visible Infrared Imaging Radiometer Suite (VIIRS) AOD fall within the expected accuracy interval, with an overall bias of 0.13 (or~28%) [5]. The MODIS latest Collection 6.1 data also show unsatisfactory performance over Asia. The percentage of AOD retrievals falling within the expected error envelope is~57% compared to 66% globally, and the root mean square error can be as large as 0.176~0.194 over different surfaces [6]. According to previous studies, the AOD errors can be attributed to several main factors, including cloud screening [3], surface reflectance parameterization [7], aerosol model assumptions [7][8][9][10], and aerosol vertical profile assumptions [10,11].
To retrieve aerosol loading from passive satellite platforms, one needs to establish the relationship between satellite observed reflected radiance and AOD. The upward radiation observed by the satellite is composed of the scattered radiance by the atmosphere and reflected radiance by the surface that transmitted through the atmosphere and can be expressed as the following equation: where θ s , θ v , ϕ represent the solar zenith angle, the viewing zenith angle and the relative azimuth angle, respectively; F λ (θ s ), T λ (θ v ) are the transmissions of upward and downward radiance, respectively. S λ is the atmospheric backscattering ratios. ρ sat λ , ρ atm λ and ρ surf λ are the TOA (Top of Atmosphere) reflectance, atmospheric path reflectance and surface reflectance at wavelength λ, respectively. Most of the satellite aerosol retrieval algorithms are built on the look-up table (LUT) strategy (i.e., MODIS, VIIRS, etc.). The LUTs contain a series of ρ sat λ calculated under a variety of geometries, aerosol types and aerosol loadings. Major aerosol properties, including AOD and aerosol model, are retrieved by comparing the satellite observed spectral reflectance with the simulated TOA reflectance stored in LUTs.
In the LUT algorithm, aerosol vertical profiles need to be assumed to calculate the TOA reflectance and all the items on the right-hand side of Equation (1) except ρ surf λ will change with the profile assumption. These are typically described using an exponential distribution or Gaussian distribution, and represented by the scale height h, the altitude below which 1-1/e (63%) of the total AOD is present. Different algorithms may assume different h values. For example, all the spherical aerosol types are assumed to be distributed below 10 km with a layer scale height of 2 km in the MISR (Multiangle Imaging Spectroradiometer) algorithm [12]. The VIIRS and MODIS C6_DT algorithms also applied the exponential distribution with a scale height of 2 km over land [13,14], whereas the Gaussian function is assumed in the MODIS Deep Blue algorithm and POLDER (POLarization and Directionality of the Earth's Reflectances) aerosol retrieval algorithm [15,16]. However, because the real aerosol profile may differ from these assumptions in both the shape and the scale height, the retrieved AOD may be erroneous due to the incorrect aerosol vertical profile used in the algorithm. Wu et al. (2016) suggested that a negative AOD bias of 10% can be caused for elevated smoke or dust layers in the MODIS C6_DT algorithm [10], and, by parameterizing a Gaussian distribution profile with a varying mean height Remote Sens. 2020, 12, 1524 3 of 20 derived from CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation) data, the negative bias can be reduced by 3%-5% over these regions [11].
Nonetheless, there is still a lack of systematic assessment of the impact of aerosol vertical distribution on the accuracy of AOD retrieval. This includes the consideration of different aerosol types, different vertical structures and different underlying surfaces. In particular, as pointed out by previous studies [17,18], AOD retrieval over East Asia, especially China, is still less than satisfactory. It is thus necessary to investigate the aerosol vertical distribution and how this may contribute to the AOD uncertainty over this region. For this purpose, we conduct a series of sensitivity experiments to investigate the impact of aerosol vertical distribution on AOD retrieval accuracy by incorporating different aerosol types and different aerosol profiles. We also examine the effect of replacing lidar retrieved profiles with the exponential profiles. We hope that this study can provide a basis for future algorithm improvements.

Description of Radiative Transfer Models
The impact of aerosol vertical profile on AOD retrieval is first investigated through a series of sensitivity experiments using the 6SV and MODerate resolution atmospheric TRANsmission (MODTRAN) radiative transfer models. The two radiative transfer models will be briefly described below.

6SV Radiative Transfer Model
The vector version of the Second Simulation of the Satellite Signal in the Solar Spectrum (6SV) is a computer code that can accurately simulate the radiative transfer process. It is capable of accounting for the polarization of radiation, elevated targets and modeling the molecular/aerosol/mixed atmosphere. 6SV calculates Rayleigh and aerosol scattering based on the Successive Orders of Scattering Method. It enables the simulation of satellite and aircraft observations and has been widely used in calculating LUTs for satellite aerosol retrieval algorithms, such as MODIS atmospheric correction algorithms and VIIRS retrieval algorithm [19]. Validation has been made against Bréon's Monte Carlo code and Coulson's tabulated values, which indicates agreement to better than 1% [20].

MODTRAN Radiative Transfer Model
MODerate resolution atmospheric TRANsmission (MODTRAN) [21] is a scalar radiative transfer code developed and maintained by the Air Force Research Laboratory (AFRL) and Spectral Science Inc (SSI). The MODTRAN code simulates radiative transfer processes and calculates the atmospheric transmittance and radiance for wavelengths extending from the ultraviolet to the microwave range (0.2~10,000 µm) with a spectral resolution down to 2 cm −1 . Here, we use MODTRAN5.2.2 [22] with higher accuracy to perform experiments with different aerosol types located at different heights.

VIIRS Level 1B Data
In addition to sensitivity experiments, we also perform retrieval experiments with satellite-measured TOA reflectance to test the effect of using observed aerosol vertical profiles. The satellite data used is the Level 1B TOA reflectance data from the VIIRS sensor. VIIRS Level 1B, also known as Sensor Data Records (SDRs), are produced from the Raw Data Records (RDRs), and they contain the calibrated and geolocated TOA brightness temperatures, radiance and reflectance. VIIRS has 22 SDRs, including 16 moderate-resolution bands (M-bands) products with a spatial resolution of 750 m at nadir, five imaging-resolution bands (I-bands) products with a resolution of 350 m at nadir and one Day/Night band (DNB) products with a resolution of 750 m throughout the scan. A previous study showed that the VIIRS reflective solar band radiometric uncertainties in reflectance are within 2% [23]. For simplicity, we perform a dual wavelength retrieval using the 488 and 672 nm channels.

Micro-Pulse Lidar Aerosol Extinction Profiles
The micro-pulse lidar (MPL) is a useful tool for retrieving the vertical distribution of aerosol extinction with a high vertical and temporal resolution. An MPL device was installed and has been operated by Peking University since July 2016 in Beijing (39.99 • N, 116.31 • E), the largest megacity in north China. It has a temporal resolution of 15 seconds and a vertical resolution of 15 m, with a 150 m blind zone. The MPL measurements have been averaged into 30-min intervals to improve the signal-to-noise ratio. Only aerosol extinction retrievals below 5 km under no precipitation conditions are used in this study to ensure the quality of MPL data. In total, 171 days of data are selected. Detailed description on the extinction profile retrieval algorithm can be found in [24]. In this study, we utilize the aerosol extinction coefficient profiles provided in Version 4.20 of the CALIOP Level 2 product [25]. Specifically, the extinction profiles within a 4 • spatial window of the Beijing-PKU site are extracted and averaged. Moreover, basic quality screening strategies are applied. Specifically, if the extinction QC (Quality Control) Flag is greater than 1, extinction uncertainty is greater than 99.9, or, if aerosol layers whose CAD (Cloud-Aerosol Discrimination) scores are greater than −20, data will be screened out.

AERONET Data
The Aerosol Robotic Network (AERONET) is a ground-based, globally distributed network providing long-term, wide-range observations of aerosol optical, microphysical and radiative properties [26]. The AERONET measurements are used here for two purposes: (1) to provide the optical parameters used in the retrieval experiment with observed profiles and (2) to evaluate satellite retrieval results.
AERONET uses a CIMEL Electronique 318 spectral radiometer to measure direct and diffuse sun radiation to retrieve aerosol optical properties. The total uncertainty of AERONET AOD under cloud-free conditions is < ± 0.01 for λ > 440 nm [26]. We use Version 3 Level 2 quality assured data at the Beijing-PKU site, and the AOD is interpolated to 550 nm using a 2nd-order polynomial fit in logarithmic coordinates [27]. In addition to AOD, inversion products from June 2016 to December 2018, including real and imaginary parts of the complex refractive index, are used, in addition to size distribution parameters, to determine the local aerosol models used in the retrieval experiment. The ground-based parameters are calculated as the temporal averages within ±1 h from the satellite overpass time. Apart from the Beijing-PKU site, AERONET observations from other 92 sites are also used in this study to relate the global distribution of AOD and SSA (Single Scattering Albedo) values to our sensitivity results.

Impact of Aerosol Scale Height Assuming Exponential Profile
We first investigate the impact of aerosol scale height in retrieval experiments using a perturbation approach. Specifically, we assume that the true scale height is 2 km, and simulate the TOA reflectance at 488 and 672 µm (corresponding to the M3 and M5 channels of VIIRS, respectively) observed by satellites at nadir view with 200 aerosol loadings ranging from 0.01 to 2.00 using the 6SV radiative Remote Sens. 2020, 12, 1524 5 of 20 transfer model. We then retrieve AOD with these simulated TOA reflectances but perturb the scale height by ±1 km and calculate the retrieval error as: where τ ref is the reference AOD value used to simulate TOA reflectance with scale height set to 2 km, and τ is the retrieved AOD value by matching the TOA reflectance simulated with different scale heights to the reference TOA reflectance. We consider three representative aerosol types, namely desert dust, fine scattering aerosol and fine absorbing aerosols, and four surface reflectances ranging from 0.02 to 0.2. Figure 1 shows the assumed vertical profiles with different scale heights used in the study, from which it is clearly seen that as scale height decreases, more aerosol is concentrated in the lower atmosphere. Figure 2 shows the spectral variability of single scattering albedo (SSA) and asymmetry parameter (g) for the three aerosol types embedded within the 6SV model, namely continental, urban and desert dust aerosols.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 23 We consider three representative aerosol types, namely desert dust, fine scattering aerosol and fine absorbing aerosols, and four surface reflectances ranging from 0.02 to 0.2. Figure 1 shows the assumed vertical profiles with different scale heights used in the study, from which it is clearly seen that as scale height decreases, more aerosol is concentrated in the lower atmosphere. Figure 2 shows the spectral variability of single scattering albedo (SSA) and asymmetry parameter (g) for the three aerosol types embedded within the 6SV model, namely continental, urban and desert dust aerosols.   Figure 3a presents the relationship between AOD error and the aerosol scale height error for three aerosol types assuming AOD = 0.2 and surface albedo = 0.05. Note the true scale height is assumed to be 2 km and the error is defined as the difference between the scale height used in the retrieval and the true scale height. It is clearly seen that the error in the retrieved AOD increases as a function of the error in scale height, and the impact of scale height is the strongest for fine absorbing aerosols, i.e., a scale height error of 1 km can lead to an AOD error exceeding 30%. The effect of negative and positive scale height error is asymmetric, indicating the non-linear response of the AOD We consider three representative aerosol types, namely desert dust, fine scattering aerosol and fine absorbing aerosols, and four surface reflectances ranging from 0.02 to 0.2. Figure 1 shows the assumed vertical profiles with different scale heights used in the study, from which it is clearly seen that as scale height decreases, more aerosol is concentrated in the lower atmosphere. Figure 2 shows the spectral variability of single scattering albedo (SSA) and asymmetry parameter (g) for the three aerosol types embedded within the 6SV model, namely continental, urban and desert dust aerosols.   Figure 3a presents the relationship between AOD error and the aerosol scale height error for three aerosol types assuming AOD = 0.2 and surface albedo = 0.05. Note the true scale height is assumed to be 2 km and the error is defined as the difference between the scale height used in the retrieval and the true scale height. It is clearly seen that the error in the retrieved AOD increases as a function of the error in scale height, and the impact of scale height is the strongest for fine absorbing aerosols, i.e., a scale height error of 1 km can lead to an AOD error exceeding 30%. The effect of negative and positive scale height error is asymmetric, indicating the non-linear response of the AOD  be 2 km and the error is defined as the difference between the scale height used in the retrieval and the true scale height. It is clearly seen that the error in the retrieved AOD increases as a function of the error in scale height, and the impact of scale height is the strongest for fine absorbing aerosols, i.e., a scale height error of 1 km can lead to an AOD error exceeding 30%. The effect of negative and positive scale height error is asymmetric, indicating the non-linear response of the AOD to the aerosol vertical distribution. Moreover, since the AOD retrieval uncertainty is usually required to be within 20% for satellite remote sensing, according to Figure 3a, the scale height uncertainty needs to be confined within ±0.8 km for absorbing aerosols. For scattering aerosols and dust, AOD is not as sensitive to scale height, and the AOD uncertainty is below 3% with ±1 km of scale height variation. The reason for the different behavior of absorbing and scattering aerosols is that, for pure scattering aerosols, their effect is an overall increase in planetary albedo and is only related to column averaged optical properties. For absorbing aerosols, the mechanism is more complicated and involves the balance between aerosol loading and the length of the absorption light path. For a certain aerosol layer at any height, the light will experience the absorption by upper aerosols, scattering and absorption by this aerosol layer, and the re-absorption by the upper aerosol layer. If the scale height is biased low, aerosol loading in the upper layers is less but the absorption light path is longer, and vice versa for high scale height bias. The effect of aerosol loading and light path cancels out at~∆scale height = 0.1 km, where AOD error equals zero (Figure 3a), below or above which the wrong assumption of scale height will both introduce a high bias in the scattered radiation and thus a low bias in AOD. However, the relationship shown in Figure 3a also depends on surface albedo and will be discussed later. The absolute value of AOD also has some impact on the AOD-scale height error relationship, i.e., the relative error gradually decreases with the increase in AOD when AOD < 0.1 but remains stable when AOD further increases (Figure 3b).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 23 to the aerosol vertical distribution. Moreover, since the AOD retrieval uncertainty is usually required to be within 20% for satellite remote sensing, according to Figure 3a, the scale height uncertainty needs to be confined within ±0.8 km for absorbing aerosols. For scattering aerosols and dust, AOD is not as sensitive to scale height, and the AOD uncertainty is below 3% with ± 1km of scale height variation. The reason for the different behavior of absorbing and scattering aerosols is that, for pure scattering aerosols, their effect is an overall increase in planetary albedo and is only related to column averaged optical properties. For absorbing aerosols, the mechanism is more complicated and involves the balance between aerosol loading and the length of the absorption light path. For a certain aerosol layer at any height, the light will experience the absorption by upper aerosols, scattering and absorption by this aerosol layer, and the re-absorption by the upper aerosol layer. If the scale height is biased low, aerosol loading in the upper layers is less but the absorption light path is longer, and vice versa for high scale height bias. The effect of aerosol loading and light path cancels out at ~∆scale height = 0.1 km, where AOD error equals zero (Figure 3a), below or above which the wrong assumption of scale height will both introduce a high bias in the scattered radiation and thus a low bias in AOD. However, the relationship shown in Figure 3a also depends on surface albedo and will be discussed later. The absolute value of AOD also has some impact on the AOD-scale height error relationship, i.e., the relative error gradually decreases with the increase in AOD when AOD < 0.1 but remains stable when AOD further increases (Figure 3b). The relationship in Figure 3a also depends on surface albedo. For absorbing aerosols, the error first increases when surface albedo increases from 0.02 to 0.05. Interestingly, when surface albedo further increases to 0.10, the Δ%AOD-Δscale height relationship becomes reversed, and when it reaches 0.20, the relationship turns nearly flat (Figure 4). A possible explanation of this phenomenon is that the TOA reflectance observed by satellite comprises the reflectance from both the atmosphere and surface. For lower surface albedos, the reflectance from atmosphere is the dominant component of the TOA reflectance observed by satellite. If the scale height is biased high, the absorbing aerosols are higher up in the atmosphere, hence more solar radiation will be absorbed by aerosols and more aerosol loadings are needed to compensate for the TOA reflectance observed by the satellites. However, as the surface albedo increases, the TOA reflectance will be dominated by the reflected radiation from the surface. If we assume that absorbing aerosols are located higher in the atmosphere (positive Δscale height), in the radiative transfer calculation, there will be more radiation for them to absorb and less radiation reflected back to the satellite, which means the retrieval algorithm will produce a larger AOD to compensate for the radiation loss by the extra absorption. The reverse is true for negative Δscale height. On the other hand, when surface albedo is high, the reflected radiation from the surface becomes the dominant term, and a negative Δscale height means more absorption The relationship in Figure 3a also depends on surface albedo. For absorbing aerosols, the error first increases when surface albedo increases from 0.02 to 0.05. Interestingly, when surface albedo further increases to 0.10, the ∆%AOD-∆scale height relationship becomes reversed, and when it reaches 0.20, the relationship turns nearly flat (Figure 4). A possible explanation of this phenomenon is that the TOA reflectance observed by satellite comprises the reflectance from both the atmosphere and surface. For lower surface albedos, the reflectance from atmosphere is the dominant component of the TOA reflectance observed by satellite. If the scale height is biased high, the absorbing aerosols are higher up in the atmosphere, hence more solar radiation will be absorbed by aerosols and more aerosol loadings are needed to compensate for the TOA reflectance observed by the satellites. However, as the surface albedo increases, the TOA reflectance will be dominated by the reflected radiation from the surface. If we assume that absorbing aerosols are located higher in the atmosphere (positive ∆scale height), in the radiative transfer calculation, there will be more radiation for them to absorb and less radiation reflected back to the satellite, which means the retrieval algorithm will produce a larger AOD to compensate for the radiation loss by the extra absorption. The reverse is true for negative ∆scale height. On the other hand, when surface albedo is high, the reflected radiation from the surface becomes the dominant term, and a negative ∆scale height means more absorption of radiation from the surface, thus leading to a positive ∆AOD. However, when surface albedo becomes very high, the signal from aerosol will be dwarfed by the surface and the aerosol vertical distribution will make a negligible impact. The results for scattering aerosols and dust are not sensitive to surface albedo and will thus be skipped here.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 distribution will make a negligible impact. The results for scattering aerosols and dust are not sensitive to surface albedo and will thus be skipped here.

Impact of the Planetary Boundary Layer
In reality, the aerosol vertical profile may differ from the exponential shape. A typical case is the existence of a lower atmosphere planetary boundary layer (PBL). The aerosols tend to distribute relatively uniformly within the PBL, but their concentration can decrease sharply above it. This structure is frequently observed in the Beijing area, and Figure 5 gives several examples from MPL measurements, i.e., the bulk of pollutants is concentrated within the near-surface layer and decreases rapidly from the altitude of 1~2km. Because the PBL is generally not considered in the retrieval algorithm, it is necessary to assess the impact, on AOD retrieval, of neglecting the PBL when it actually exists. We therefore perform the following experiment. Two models are assumed to characterize aerosol vertical distribution (see Figure 6). The first model (blue line) only considers an exponentially decreasing profile starting from the surface, whereas the second model (red line) includes a 1-km-thick boundary layer within which the aerosols are well mixed. We perturb the scale heights from 1~3 km with a 0.1-km step size. Errors are defined by Equation (2) as well, except that the is the reference AOD value used to simulate TOA reflectance with the boundary layer, and is the retrieved AOD value by matching the TOA reflectance simulated without the boundary layer to the reference TOA reflectance.

Impact of the Planetary Boundary Layer
In reality, the aerosol vertical profile may differ from the exponential shape. A typical case is the existence of a lower atmosphere planetary boundary layer (PBL). The aerosols tend to distribute relatively uniformly within the PBL, but their concentration can decrease sharply above it. This structure is frequently observed in the Beijing area, and Figure 5 gives several examples from MPL measurements, i.e., the bulk of pollutants is concentrated within the near-surface layer and decreases rapidly from the altitude of 1~2 km. Because the PBL is generally not considered in the retrieval algorithm, it is necessary to assess the impact, on AOD retrieval, of neglecting the PBL when it actually exists. We therefore perform the following experiment. Two models are assumed to characterize aerosol vertical distribution (see Figure 6). The first model (blue line) only considers an exponentially decreasing profile starting from the surface, whereas the second model (red line) includes a 1-km-thick boundary layer within which the aerosols are well mixed. We perturb the scale heights from 1~3 km with a 0.1-km step size. Errors are defined by Equation (2) as well, except that the τ ref is the reference AOD value used to simulate TOA reflectance with the boundary layer, and τ is the retrieved AOD value by matching the TOA reflectance simulated without the boundary layer to the reference TOA reflectance.   The AOD errors are shown in Figure 7 (note that Figure 7 only shows the results for fine absorbing aerosols because, for scattering aerosols and dust, AOD is not sensitive to vertical distribution assumptions). It is seen that the retrieval error can be as large as~10% for AOD = 0.2 if ignoring the PBL in the retrieval, but will decrease to~5% for AOD > 0.5. Similar to Figure 4, the error is also sensitive to the surface albedo, i.e., when the surface albedo is lower than 0.1, AOD errors are almost negative, but they turn to positive when surface albedo increases. When surface albedo becomes sufficiently large (0.20), the AOD errors become negligible. Even with the same scale height, considering the PBL means that the aerosols will be distributed closer to the surface (see Figure 6) than those without considering the PBL. It is thus easy to interpret the magnitude and sign change of the AOD error using the mechanism discussed in the last section. In addition, similar to Figure 3b, the error is a function of AOD loading: the relative error can reach above 20% when AOD loading is rather small (< 0.1), but it can be confined within 5% for AOD = 1. In summary, our results indicate that the PBL structures can have a nonnegligible impact on the accuracy of AOD retrievals. To constrain AOD uncertainty within 20%, it is better to incorporate a profile with PBL in the retrieval algorithm.
The AOD errors are shown in Figure 7 (note that Figure 7 only shows the results for fine absorbing aerosols because, for scattering aerosols and dust, AOD is not sensitive to vertical distribution assumptions). It is seen that the retrieval error can be as large as ~10% for AOD = 0.2 if ignoring the PBL in the retrieval, but will decrease to ~5% for AOD > 0.5. Similar to Figure 4, the error is also sensitive to the surface albedo, i.e., when the surface albedo is lower than 0.1, AOD errors are almost negative, but they turn to positive when surface albedo increases. When surface albedo becomes sufficiently large (0.20), the AOD errors become negligible. Even with the same scale height, considering the PBL means that the aerosols will be distributed closer to the surface (see Figure 6) than those without considering the PBL. It is thus easy to interpret the magnitude and sign change of the AOD error using the mechanism discussed in the last section. In addition, similar to Figure 3b, the error is a function of AOD loading: the relative error can reach above 20% when AOD loading is rather small (< 0.1), but it can be confined within 5% for AOD = 1. In summary, our results indicate that the PBL structures can have a nonnegligible impact on the accuracy of AOD retrievals. To constrain AOD uncertainty within 20%, it is better to incorporate a profile with PBL in the retrieval algorithm.

Impact of Layered Aerosol Vertical Structure
Another possible scenario of aerosol vertical distribution is the existence of two aerosol layers with different aerosol types, such as elevated smoke layer transported over polluted regions. Figure  8 shows two examples of this phenomenon from CALIPSO observations. The first is from South America, where elevated smoke aerosols at around 3 km are observed over polluted continental aerosols. The second example shows a similar structure with smoke located above polluted dust

Impact of Layered Aerosol Vertical Structure
Another possible scenario of aerosol vertical distribution is the existence of two aerosol layers with different aerosol types, such as elevated smoke layer transported over polluted regions. Figure 8 shows two examples of this phenomenon from CALIPSO observations. The first is from South America, where elevated smoke aerosols at around 3 km are observed over polluted continental aerosols. The second example shows a similar structure with smoke located above polluted dust aerosols over South Africa. We thus investigate the impact of using one exponential profile when the aerosol distribution is actually layered in this section. The MODTRAN model is used for this experiment because the 6SV model does not allow prescribing two aerosol profiles.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 23 aerosols over South Africa. We thus investigate the impact of using one exponential profile when the aerosol distribution is actually layered in this section. The MODTRAN model is used for this experiment because the 6SV model does not allow prescribing two aerosol profiles. In the experiment, we design two ideal layered aerosol structures, as depicted in Figure 9. The first structure has soot lying between 0-1 km and sulfate aerosols (under 50% relative humidity) lying between 3-4 km above the surface (case 1). The second structure is reversed with soot lying above sulfate (case 2). The optical parameters for sulfate and soot aerosols are obtained from the Optical Properties of Aerosols and Clouds (OPAC) database [28], and their SSA and g values at different wavelengths are shown in Figure 10. In the experiment, we design two ideal layered aerosol structures, as depicted in Figure 9. The first structure has soot lying between 0-1 km and sulfate aerosols (under 50% relative humidity) lying between 3-4 km above the surface (case 1). The second structure is reversed with soot lying above sulfate (case 2). The optical parameters for sulfate and soot aerosols are obtained from the Optical Properties of Aerosols and Clouds (OPAC) database [28], and their SSA and g values at different wavelengths are shown in Figure 10.    In both cases, the optical depth fractions for soot and sulfate are 10% and 90%, respectively. We first calculate the TOA reflectance at nadir view for these two cases, and then use these TOA reflectances to retrieve AOD, assuming an exponential profile with a scale height of 2 km. The SSA and g used are equivalent values, assuming the external mixing of the two aerosol components. Figure 11 shows the retrieved AOD error, assuming a mixed exponential profile relative to the two layered cases. First, it is seen that the AOD errors are reversed for the two cases, i.e., errors are positive when soot lies above sulfate and negative when soot lies below sulfate. The explanation of this phenomenon could be that in case 1, sulfate above blocks solar radiation from reaching the black carbon layer, thus reducing its absorptivity and increases TOA reflectance, whereas, for case 2, scattering aerosol underneath will act as a bright surface, which enhances black carbon absorption and thus reduces the TOA reflectance compared to the mixed case. The errors are also quite large and asymmetric, with larger errors for the soot above sulfate case. For case 2, when soot lies above, errors Figure 10. Spectral dependence of single scattering albedo (a) and asymmetry factor (b) for soot and sulfate (50% humidity) obtained from the OPAC database.
In both cases, the optical depth fractions for soot and sulfate are 10% and 90%, respectively. We first calculate the TOA reflectance at nadir view for these two cases, and then use these TOA reflectances to retrieve AOD, assuming an exponential profile with a scale height of 2 km. The SSA and g used are equivalent values, assuming the external mixing of the two aerosol components. Figure 11 shows the retrieved AOD error, assuming a mixed exponential profile relative to the two layered cases. First, it is seen that the AOD errors are reversed for the two cases, i.e., errors are positive when soot lies above sulfate and negative when soot lies below sulfate. The explanation of this phenomenon could be that in case 1, sulfate above blocks solar radiation from reaching the black carbon layer, thus reducing its absorptivity and increases TOA reflectance, whereas, for case 2, scattering aerosol underneath will act as a bright surface, which enhances black carbon absorption and thus reduces the TOA reflectance compared to the mixed case. The errors are also quite large and asymmetric, with larger errors for the soot above sulfate case. For case 2, when soot lies above, errors can reach 28% for AOD = 0.5 and exceed 85% for a heavily polluted situation (AOD = 1.0). For case 1, the error is −18% when AOD = 0.5 and reaches −40% when AOD = 1.0. These results reflect the nonlinear radiative interaction of different aerosol types within the column, and that more diversified profiles will help improve AOD retrievals.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 23 can reach 28% for AOD = 0.5 and exceed 85% for a heavily polluted situation (AOD = 1.0). For case 1, the error is -18% when AOD = 0.5 and reaches −40% when AOD = 1.0. These results reflect the nonlinear radiative interaction of different aerosol types within the column, and that more diversified profiles will help improve AOD retrievals.

AOD Retrieval using Observed Aerosol Vertical Profiles
In the above discussion, we have shown that the AOD retrieval error is highly sensitive to aerosol profile, especially for absorbing aerosols. In addition to sensitivity experiments, we further explore whether the retrieved AOD will improve if we use observed aerosol profiles in the algorithm. We focus on the Beijing-PKU site, because complete aerosol optical and vertical measurements by the sun photometer and MPL lidar are available here.
AOD is retrieved using Level 1B TOA reflectance data from the VIIRS sensor onboard the Suomi-

AOD Retrieval Using Observed Aerosol Vertical Profiles
In the above discussion, we have shown that the AOD retrieval error is highly sensitive to aerosol profile, especially for absorbing aerosols. In addition to sensitivity experiments, we further explore whether the retrieved AOD will improve if we use observed aerosol profiles in the algorithm. We focus on the Beijing-PKU site, because complete aerosol optical and vertical measurements by the sun photometer and MPL lidar are available here.
AOD is retrieved using Level 1B TOA reflectance data from the VIIRS sensor onboard the Suomi-NPP satellite. The geometric conditions (solar, satellite zenith angle and relative azimuth angle) were pre-calculated for the site location and satellite overpass time. Our retrieval algorithm is based on the LUT approach: the LUT is pre-calculated for a set of geometrics, aerosol loading, aerosol model, and surface parameters using the 6SV radiative transfer code. To minimize the impact on the retrieval results by other factors such as aerosol model and surface albedo assumptions, we use aerosol optical properties and surface albedo from the AERONET inversion products. The surface is assumed to be Lambertian and the surface reflectance is obtained from sun-photometer observations. The aerosol models were re-constructed using sun-photometer measurements, spectral SSA and g values of which are shown in Figure 12. For comparison, we generate two sets of LUTs by using the real aerosol extinction profiles retrieved by the MPL during the satellite overpass time and the exponential profiles with a scale height of 2 km, respectively. The spectral reflectance observed by Suomi-NPP is compared with the spectral reflectance from the LUT, and the retrieved AOD is determined by minimizing the residual constructed from the difference in the M3 (488 ) and M5 (672 ) channels. The retrieval results using these two LUTs are subsequently evaluated against AERONET AOD. In total, 111 cases are successfully retrieved for the one-year experiment (July 2017 to June 2018). Figure 13 shows the scatter plots of retrieved AOD using the two sets of vertical profiles against AERONET AOD, and the evaluation statistics are summarized in Table 1. Overall, the results using observed profiles (red dots) exhibit better agreements with AERONET than those using the exponential profile (blue dots). The mean bias is reduced from 0.15 to 0.03, precision is reduced from 0.015 to 0.002, and correlation increased from 0.63 to 0.83. Seasonally, the improvements are the most significant in the winter season, with a greatly reduced mean bias (from 0.38 to −0.05) and RMSE (root mean square deviation) (from 0.09 to 0.01) and a greatly increased correlation (from 0.48 to 0.77) ( Figure 14 and Table 1). This is because absorbing aerosols are the most frequently observed aerosol type in the winter season due to residential heating. The two outliers in the original retrieval result at the upper left corner of Figure 13 result from large differences between the observed and assumed aerosol profiles with highly absorbing aerosols present, i.e., the 440 nm SSA for these two cases are 0.89 and 0.88, respectively. As shown in Figure 15, the winter season has generally lower SSAs compared to the other seasons, with ~65% cases having SSA < 0.90. According to Section 4.1, the impact of vertical distribution is the most significant for absorbing aerosols; therefore, more improvements are seen in the winter season when using the observed profiles. For comparison, we generate two sets of LUTs by using the real aerosol extinction profiles retrieved by the MPL during the satellite overpass time and the exponential profiles with a scale height of 2 km, respectively. The spectral reflectance observed by Suomi-NPP is compared with the spectral reflectance from the LUT, and the retrieved AOD is determined by minimizing the residual constructed from the difference in the M3 (488 µm) and M5 (672 µm) channels. The retrieval results using these two LUTs are subsequently evaluated against AERONET AOD. In total, 111 cases are successfully retrieved for the one-year experiment (July 2017 to June 2018). Figure 13 shows the scatter plots of retrieved AOD using the two sets of vertical profiles against AERONET AOD, and the evaluation statistics are summarized in Table 1. Overall, the results using observed profiles (red dots) exhibit better agreements with AERONET than those using the exponential profile (blue dots). The mean bias is reduced from 0.15 to 0.03, precision is reduced from 0.015 to 0.002, and correlation increased from 0.63 to 0.83. Seasonally, the improvements are the most significant in the winter season, with a greatly reduced mean bias (from 0.38 to −0.05) and RMSE (root mean square deviation) (from 0.09 to 0.01) and a greatly increased correlation (from 0.48 to 0.77) ( Figure 14 and Table 1). This is because absorbing aerosols are the most frequently observed aerosol type in the winter season due to residential heating. The two outliers in the original retrieval result at the upper left corner of Figure 13 result from large differences between the observed and assumed aerosol profiles with highly absorbing aerosols present, i.e., the 440 nm SSA for these two cases are 0.89 and 0.88, respectively. As shown in Figure 15, the winter season has generally lower SSAs compared to the other seasons, with~65% cases having SSA < 0.90. According to Section 4.1, the impact of vertical distribution is the most significant for absorbing aerosols; therefore, more improvements are seen in the winter season when using the observed profiles.     In addition, a similar one-year experiment (July 2017 to June 2018) using regionally averaged extinction profiles by CALIPSO around the Beijing-PKU site is also conducted as this dataset is globally accessible. Similar to Figure 13, the results in Figure 16 using CALIPSO profiles (red dots) more clearly agree with AERONET than those using assumed profiles (blue dots), with the mean bias reduced from 0.20 to 0.14 and correlation coefficient increased from 0.53 to 0.80. The most obvious improvements are also observed in winter, with a greatly reduced mean bias (from 0.43 to 0.01) and RMSE (from 0.14 to 0.01) and a greatly increased correlation (from −0.02 to 0.75) (Figure 17).
It should be noted that, due to the narrow footprint of CALIPSO ground track, the sample size for the CALIPSO experiment is smaller than that in Figures 13 and 14. Moreover, unlike the dualwavelength retrieval algorithm performed in the above experiment, a single-channel (532 μm ) retrieval is performed here since only 532 and 1064 nm extinction profiles are available for CALIPSO. Nonetheless, the improvements in the winter season still imply the feasibility of applying global aerosol vertical profiles in a satellite retrieval algorithm.  In addition, a similar one-year experiment (July 2017 to June 2018) using regionally averaged extinction profiles by CALIPSO around the Beijing-PKU site is also conducted as this dataset is globally accessible. Similar to Figure 13, the results in Figure 16 using CALIPSO profiles (red dots) more clearly agree with AERONET than those using assumed profiles (blue dots), with the mean bias reduced from 0.20 to 0.14 and correlation coefficient increased from 0.53 to 0.80. The most obvious improvements are also observed in winter, with a greatly reduced mean bias (from 0.43 to 0.01) and RMSE (from 0.14 to 0.01) and a greatly increased correlation (from −0.02 to 0.75) (Figure 17).
It should be noted that, due to the narrow footprint of CALIPSO ground track, the sample size for the CALIPSO experiment is smaller than that in Figures 13 and 14. Moreover, unlike the dual-wavelength retrieval algorithm performed in the above experiment, a single-channel (532 µm) retrieval is performed here since only 532 and 1064 nm extinction profiles are available for CALIPSO. Nonetheless, the improvements in the winter season still imply the feasibility of applying global aerosol vertical profiles in a satellite retrieval algorithm.

Discussion
In the previous section, we use radiative transfer models to study the sensitivity of satellite AOD retrievals to the aerosol vertical profile assumptions, and found that the effect is the most obvious for absorbing aerosols and also depends on the total aerosol loading (i.e., AOD). To fully evaluate how aerosol profile assumption may impact on AOD retrieval spatially and temporally, it is necessary to Figure 17. Comparison of the mean bias (upper panel), RMSE (middle panel) and the correlation (bottom panel) between the retrieval using the exponential profile (blue) and the CALIPSO aerosol profile (orange) for the four seasons.

Discussion
In the previous section, we use radiative transfer models to study the sensitivity of satellite AOD retrievals to the aerosol vertical profile assumptions, and found that the effect is the most obvious for absorbing aerosols and also depends on the total aerosol loading (i.e., AOD). To fully evaluate how aerosol profile assumption may impact on AOD retrieval spatially and temporally, it is necessary to discuss the results in a global context, i.e., using the sensitivity results to infer whether aerosol profile assumption is likely a significant source of uncertainty for different regions and different seasons.
We examine the distribution of major aerosol optical parameters for seven major regions, namely North America (NAM), South America (SAM), Europe, North Africa and the Middle East (NAME), South Africa (SAF), Asia and Oceania, whose spatial domains are shown in Figure 18. Figures 19 and 20 display the probability density distributions of AOD and SSA for these seven regions from AERONET data, averaged using measurements from the sites where long-term and continuous observations are available, marked by the red dots on Figure 18. Different features are observed for different regions. In general, NAME and Asia exhibit a right-shifted AOD distribution with a longer tail than the other regions, suggesting that they tend to have a larger mean AOD and more extreme values. Oceania and NAM have the lowest overall aerosol loadings (with the most left-shifted AOD distribution), and the AOD for SAF and SAM lies in between. For regions with larger aerosol loading, the issue about vertical profile assumption has an important effect because a large absolute error can occur. This effect is not negligible even for regions with smaller AOD values because the relative error can be large, as discussed above. The highest SSA are found at the Oceania sites, especially in spring (MAM) and summer (JJA), and the lowest are found in SAF, Asia and SAM. The SSA distribution over SAF is also closely related to the seasonal biomass burning, i.e., they are flat in spring and winter (DJF) but have strong peaks in summer and fall (SON) [29]. The probability density function (PDF) of SSA also exhibits a peak in the fall for the SAM, which may be associated with the peak burning season in August and September in Brazil [29]. The shapes of SSA PDF for Asia show little seasonal variability, with slightly lower values for winter due to heating [30] and higher values for summer due to secondary aerosol production and hygroscopic growth [31,32].
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 23 tail than the other regions, suggesting that they tend to have a larger mean AOD and more extreme values. Oceania and NAM have the lowest overall aerosol loadings (with the most left-shifted AOD distribution), and the AOD for SAF and SAM lies in between. For regions with larger aerosol loading, the issue about vertical profile assumption has an important effect because a large absolute error can occur. This effect is not negligible even for regions with smaller AOD values because the relative error can be large, as discussed above. The highest SSA are found at the Oceania sites, especially in spring (MAM) and summer (JJA), and the lowest are found in SAF, Asia and SAM. The SSA distribution over SAF is also closely related to the seasonal biomass burning, i.e., they are flat in spring and winter (DJF) but have strong peaks in summer and fall (SON) [29]. The probability density function (PDF) of SSA also exhibits a peak in the fall for the SAM, which may be associated with the peak burning season in August and September in Brazil [29]. The shapes of SSA PDF for Asia show little seasonal variability, with slightly lower values for winter due to heating [30] and higher values for summer due to secondary aerosol production and hygroscopic growth [31,32].