Adaption of the MODIS aerosol retrieval algorithm using airborne spectral surface reflectance measurements over urban areas: a case study

Abstract

MODIS (MOderate-resolution Imaging Spectroradiometer) retrievals of aerosol optical depth (AOD) are biased over urban areas, primarily because the reflectance characteristics of urban surfaces are different than that assumed by the retrieval algorithm. Specifically, the operational “dark-target” retrieval is tuned towards vegetated (dark) surfaces and assumes a spectral relationship to estimate the surface reflectance in blue and red wavelengths. From airborne measurements of surface reflectance over the city of Zhongshan, China, were collected that could replace the assumptions within the MODIS retrieval algorithm. The subsequent impact was tested upon two versions of the operational algorithm, Collections 5 and 6 (C5 and C6). AOD retrieval results of the operational and modified algorithms were compared for a specific case study over Zhongshan to show minor differences between them all. However, the Zhongshan-based spectral surface relationship was applied to a much larger urban sample, specifically to the MODIS data taken over Beijing between 2010 and 2014. These results were compared directly to ground-based AERONET (AErosol RObotic NETwork) measurements of AOD. A significant reduction of the differences between the AOD retrieved by the modified algorithms and AERONET was found, whereby the mean difference decreased from $0.27\pm 0.14$ for the operational C5 and $0.19\pm 0.12$ for the operational C6 to $0.10\pm 0.15$ and $-0.02\pm 0.17$ by using the modified C5 and C6 retrievals. Since the modified algorithms assume a higher contribution by the surface to the total measured reflectance from MODIS, consequently the overestimation of AOD by the operational methods is reduced. Furthermore, the sensitivity of the MODIS AOD retrieval with respect to different surface types was investigated. Radiative transfer simulations were performed to model reflectances at top of atmosphere for predefined aerosol properties. The reflectance data were used as input for the retrieval methods. It was shown that the operational MODIS AOD retrieval over land reproduces the AOD reference input of 0.85 for dark surface types (retrieved AOD $=$ 0.87 (C5)). An overestimation of AOD $=$ 0.99 is found for urban surfaces, whereas the modified C5 algorithm shows a good performance with a retrieved value of AOD $=$ 0.86.

Introduction

Aerosol particles influence various aspects of the Earth's ecosystem. By scattering or absorption of solar radiation, aerosol particles influence the Earth's radiation budget, perturbing cloud and precipitation cycles and therefore the global climate . Megacities play a particular role in this regard, because they are a significant source for anthropogenic aerosol particles . Therefore megacities are of interest in global climate and chemical models . Trends of anthropogenic aerosol particles over the 189 largest cities in the world were reported by showing a significant increase over megacities in the Indian subcontinent, the Middle East, and North China.

The column integral of the atmospheric extinction (absorption plus scattering) of solar radiation by aerosol particles is known as the aerosol optical depth (AOD). By definition, AOD is spectral (function of wavelength); however it is often reported in a mid-visible wavelength (e.g., 0.55 $\mathrm{\mu }$m). AOD can be measured from the ground by sun photometers or retrieved from measured radiances, as observed at the top-of-atmosphere (TOA) from spaceborne instruments such as the MODerate-resolution Imaging Spectroradiometer (MODIS) . Since measured upward radiance at the TOA depends on both atmospheric components and the underlying surface, AOD retrieval uncertainty depends on the relative weighting of the two signals. Uncertainty from surface reflectance assumptions dominate for low values of AOD ($<0.15$), whereas aerosol optical assumptions dominate as AOD becomes large $>0.4$ . A dark-target retrieval algorithm, therefore, relies on picking out dark surfaces due to reduced impact on the TOA signal. In general, AOD retrieval over ocean water (which can be very dark) involves lower uncertainties than over land .

However, aerosol retrieval over land is possible if one can measure, model, or otherwise constrain the impact of the surface reflectance on the TOA signal. Vegetated and dark-soiled surfaces tend to have surface reflectance that correlates spectrally, meaning that reflectance in one wavelength can be estimated from reflectance in another. In the original retrieval algorithm, which was used to process up through Collection 4 (C4) of the MODIS aerosol product, it was assumed that surface reflectances in the MODIS blue and red (e.g., 0.47 and 0.65 $\mathrm{\mu }$m) bands were 0.25 and 0.50, respectively, of the value in the shortwave-infrared (2.1 $\mathrm{\mu }$m) band. This spectral relationship was used to constrain the differential equations of the AOD retrieval. As the MODIS product has matured, from C4 to Collection 5 (C5) and to the current operational Collection 6 (C6) versions of the retrieval algorithm , additional information (e.g., vegetation index and scattering angles) has been introduced to better estimate the surface reflectance.

While the resulting dark-target AOD retrieval technique shows good performance over much of the globe, the algorithm developed for dark vegetation continues to show less accurate results over surfaces with high reflectance such as arid regions and urban surfaces . For urban surfaces in particular, and reported a positive bias of retrieved AOD over Mexico City and New York City compared to ground-based sun photometer measurements.

Several approaches have been suggested to consider the problem of unknown surface reflectance in the retrieval of AOD over urban surfaces. found an improvement of the derived AOD when the local surface reflectance was used as derived from satellite data collected during days with low values of AOD. developed an aerosol retrieval algorithm for Hong Kong and the Pearl River Delta region with a spatial resolution of 500 $\mathrm{m}$. They determined the surface reflectance with a minimum reflectance technique which increased the slope of the linear regression coefficient between AODs from MODIS and sun photometer measurements from 0.75 to 0.83. A simplified algorithm to retrieve AODs with a spatial resolution of 500 $\mathrm{m}$ for urban regions such as Beijing was developed by . The results have shown a relative mean bias (RMB; ratio of AOD from MODIS and AOD derived from sun photometer at “Beijing_RADI” AERONET site) of 0.97, whereas the standard MODIS algorithm overestimated the AOD with a RMB of 1.12. used sun photometer measurements and radiative transfer calculations for an improved representation of surface reflectance in the MODIS algorithm. They have shown for the city of Nanjing, China, a seasonal variability of the relationship between the surface reflectance at the visible (0.47 and 0.65 $\mathrm{\mu }$m) and the near-infrared wavelengths (2.11 $\mathrm{\mu }$m).

The present study also focuses on urban areas of China, but instead of using empirical information (e.g., atmospheric correction of TOA radiances) it uses surface reflectance measurements provided by low-flying aircraft. Specifically, airborne spectral measurements taken over Zhongshan (southeast China) have been used to test whether the AOD retrieval is improved when substituting measured values for the empirical surface reflectance relationships. At the same time, since there are differences between the assumptions used for the C5 and C6 versions, the sensitivity of both versions compared.

Section 2 gives a brief description of the surface reflectance assumptions of the MODIS aerosol retrieval algorithm over land (both C5 and C6) and introduces the method to derive the surface reflection properties from airborne measurements over Zhongshan. These properties and the corresponding modification of reflection assumptions in the operational C5 and C6 retrieval algorithms are presented in Sect. 3. In Sect. 4 the modified and operational algorithms are applied to MODIS measurements over Zhongshan and Beijing. Furthermore, AOD values derived from AERONET are compared to MODIS retrievals for Beijing, and a sensitivity study of the reflectance ratios based on radiative transfer simulations is performed. Conclusions are given in Sect. 5.

MODIS AOD retrieval and surface reflectance assumptions

The methodology of the current dark-target retrieval (C5) was detailed by , and the differences between the current C6 and C5 versions of the MODIS aerosol were described by . There are two MODIS instruments: one on Terra satellite (equator overpass in morning) and one on Aqua (overpass in afternoon). With a 2330 $\mathrm{km}$ wide swath, and with pixels of 1 $\mathrm{km}$ (or finer) spatial resolution at nadir, MODIS observes TOA reflectances and radiances in 36 channels (or bands), ranging from visible through thermal-infrared wavelengths. Seven bands (0.47, 0.55, 0.66, 0.86, 1.2, 1.6, and 2.1 $\mathrm{\mu }$m) are used explicitly for dark-target aerosol retrieval over land and ocean and have 500 $\mathrm{m}$ or 250 $\mathrm{m}$ resolution. The retrieval algorithm is the same for both MODIS instruments.

The aerosol retrieval uses a lookup table (LUT) approach, in that a subset of the measured reflectances are compared with pre-computed simulations (using radiative transfer) of representative surface/atmosphere “scenes”. The best-matching solution is assumed to describe the aerosol conditions (including the AOD) of the given scene. Since pixels at 500 $\mathrm{m}$ resolution can be “noisy”, the standard retrieval operates on a 20 $\times$ 20 grouping (10 $\mathrm{km}$ by 10 $\mathrm{km}$) of the original pixels. Pixels containing clouds, snow or ice, and inland water bodies are filtered out, so that the remaining pixels are sorted with respect to their brightness. Pixels with reflectances between 0.01 and 0.25 at 2.1 $\mathrm{\mu }$m wavelength are selected, and the brightest 50 $\mathrm{\%}$ and darkest 20 $\mathrm{\%}$ of these pixels are discarded in order to further reduce possible cloud and surface contamination. If there are at least 12 remaining pixels (out of the original 400) within the 10 $\mathrm{km}$ by 10 $\mathrm{km}$ box, the reflectance in each channel is averaged and assumed to represent the TOA reflectance of this box. Note that, in C6, there is a separate 3 $\mathrm{km}$ resolution retrieve which is not considered here.

Outside of the uncertainties of cloud masking and pixel selection, there are various theoretical sources of error in the AOD retrieval. One is the aerosol “type” used to power the radiative transfer code. The aerosol type is described by an assumed shape, a size distribution, and refractive indices used to calculate scattering and absorption properties. Changes from C5 to C6 included corrections related to the Rayleigh/aerosol LUTs and modifications of the assumed aerosol type models that are prescribed for season and location . The other major uncertainty is related to the assumptions of surface reflectance, and its contribution to the TOA reflectance. While all aspects were considered during the upgrade from C5 to C6, here the sensitivity of the aerosol retrieval on its assumptions related to surface reflectance is studied, focusing on the assumed surface reflectance relationship and whether offline measurements (e.g., from aircraft) help to improve the aerosol retrieval.

Because of atmospheric scattering and absorption, spectral surface reflectance cannot be determined directly from satellite measurements. The original formulation (including the C4 products) of the MODIS aerosol retrieval algorithm was, in fact, based on aircraft measurements of surface reflectance over a few sites. In these studies, determined that for vegetation and dark-soiled surfaces, the surface reflectances in blue (0.49 $\mathrm{\mu }$m), red (0.66 $\mathrm{\mu }$m), and in the near infrared (2.1 $\mathrm{\mu }$m) were related by empirical relationships: ${\mathit{\rho }}_{0.49}={\mathit{\rho }}_{2.1}/\mathrm{4}$ and ${\mathit{\rho }}_{0.66}={\mathit{\rho }}_{2.1}/\mathrm{2}$. While these relationships were based on “prototype” MODIS channels (note small differences in wavelength from actual MODIS wavelengths), these ratios were assumed in operation. Therefore, since the atmosphere (except for under heavy dust conditions) is nearly transparent at 2.1 $\mathrm{\mu }$m, the surface reflectance is approximated by the measured TOA reflectance at 2.1 $\mathrm{\mu }$m (${\mathit{\rho }}_{\text{s,2.1}}$ = ${\mathit{\rho }}_{\text{TOA,2.1}}$), and the values at the other wavelengths can be estimated. The difference between TOA and surface in blue and red channels is related to aerosol.

relaxed the assumption of atmospheric transparency at 2.1 $\mathrm{\mu }$m and reformulated the aerosol retrieval as an inversion of three bands (0.47, 0.65, and 2.1 $\mathrm{\mu }$m). In addition, detailed atmospheric correction was performed over many AERONET sites, so that the surface reflectance relationship became better characterized. In particular, considered more variables for characterizing the wavelength-dependent surface reflectance, including scattering angle ${\mathit{\vartheta }}_{\text{scat}}$ and general “greenness” of the surface. The effect of geographical and seasonal variations on surface vegetation could be diagnosed through the normalized differenced vegetation index ${\text{NDVI}}_{\text{SWIR}}$, combining the measured reflectances ${\mathit{\rho }}^{\mathrm{m}}$ at 1.2 and 2.1 $\mathrm{\mu }$m wavelength: $${\text{NDVI}}_{\text{SWIR}}={\displaystyle \frac{{\mathit{\rho }}_{1.2}^{\mathrm{m}}-{\mathit{\rho }}_{2.1}^{\mathrm{m}}}{{\mathit{\rho }}_{1.2}^{\mathrm{m}}+{\mathit{\rho }}_{2.1}^{\mathrm{m}}}}.$$ The surface reflectance at 2.1 $\mathrm{\mu }$m wavelength is calculated as a function of measured TOA reflectance, the atmospheric backscatter ratio, the transmissivity, and the TOA reflectance given by the precalculated LUTs depending on AOD and aerosol type. According to the surface reflectance at 0.65 and 0.47 $\mathrm{\mu }$m is a function of the surface reflectance at 2.1 $\mathrm{\mu }$m wavelength and is derived by $$\begin{array}{}& & {\mathit{\rho }}_{\mathrm{s},0.65}={\mathit{\rho }}_{\mathrm{s},2.1}\cdot a_{0.65/2.1}+b_{0.65/2.1},& & {\mathit{\rho }}_{\mathrm{s},0.47}={\mathit{\rho }}_{\mathrm{s},0.65}\cdot a_{0.47/0.65}+b_{0.47/0.65}.\end{array}$$ The corresponding slopes a and intercepts b are given by $$\begin{array}{cc}& a_{0.65/2.1}=a_{0.65/2.1}^{\text{NDVI}}+0.002\cdot {\mathit{\vartheta }}_{\text{scat}}-0.27,\\ & b_{0.65/2.1}=-0.00025\cdot {\mathit{\vartheta }}_{\text{scat}}+0.033,\\ & a_{0.47/0.65}=0.49,\\ & & b_{0.47/0.65}=0.005,\end{array}$$ with the scattering angle ${\mathit{\vartheta }}_{\text{scat}}$ and the slope $a_{\text{0.65/2.1}}^{{\text{NDVI}}_{\text{SWIR}}}$ with respect to the vegetation index NDVI${}_{\text{SWIR}}$. Note that due to a coding error, the C5 version's NDVI${}_{\text{SWIR}}$ dependence as given in Eq. (10) from was reversed. This error was corrected in C6 . Values of NDVI${}_{\text{SWIR}}$ dependence for both versions are provided here as Eqs. (5), (6), and (7). $$\begin{array}{cc}& a_{0.65/2.1}^{\text{NDVI}}=\text{0.48 (C5) and }{a}_{0.65/2.1}^{\text{NDVI}}=\text{0.58 (C6)}\\ & & \mathrm{\text{for: NDVI}}<0.25& a_{0.65/2.1}^{\text{NDVI}}=\text{0.58 (C5) and }{a}_{0.65/2.1}^{\text{NDVI}}=\text{0.48 (C6)}\\ & & \mathrm{\text{for: NDVI}}>0.75& a_{0.65/2.1}^{\text{NDVI}}=0.48+0.2\cdot (\text{NDVI}-\text{0.25) (C5) }\\ & a_{0.65/2.1}^{\text{NDVI}}=0.58+0.2\cdot (\text{NDVI}-\text{0.25) (C6)}\\ & & \mathrm{\text{for: }}0.25\le \text{NDVI}\le 0.75\end{array}$$ As noted in the previous section, the aerosol retrieval involves many assumptions about the observations (measured spectral reflectance), aerosol types (size, shape, composition), and other factors. However, for the sensitivity studies presented in the following sections, only the equations that relate to surface reflectance are modified. Therefore, the “stand-alone code” of the MODIS aerosol retrieval algorithms (http://darktarget.gsfc.nasa.gov/reference/code) was employed. These codes provide a single output AOD for a given set of input values, which include measured TOA spectral reflectance and observing geometry. Since the operational MODIS retrieval product (e.g., either C5 or C6) provides this input information, all assumptions about pixel aggregation, cloud masking, and pixel selection are already accounted for. Therefore, this study is focused on comparing the in-field measurements with C5 and C6 assumptions, as well as testing the retrieval sensitivity to the surface assumptions.

Airborne measurements of urban surface reflection properties

Airborne measurements of surface reflection properties were performed over Zhongshan, China. The following subsections introduce the field campaign, the main instruments, and the method to retrieve the surface properties.

Field campaign in Zhongshan, Pearl River Delta, China

The city of Zhongshan, China, is in the vicinity of the megacities Guangzhou, Shenzhen, and Hong Kong. Zhongshan is characterized by a heterogeneous urban surface, representing various artificial materials as well as vegetation and small water bodies. The aerosol loading is often large. In December 2009, spectral reflectance and aerosol properties were measured over Zhongshan, using a combination of aircraft-borne spectroradiometers and surface sun photometers and LIDAR.

Spectral radiation was measured with the Spectral Modular Airborne Radiation measurement sysTem albedometer SMART albedometer;. During the campaign the instrument was equipped with two optical inlets, one for spectral upward irradiance ($F_{\mathit{\lambda }}^{\uparrow }$) and the other for spectral nadir radiance ($I_{\mathit{\lambda }}^{\uparrow }$) measurements (opening angle of radiance inlet $=$ 2.1${}^{\circ }$). Each optical inlet was connected via optical fibers with two multichannel spectrometers. The spectrometers cover a wavelength range between 0.4 and 2.1 $\mathrm{\mu }$m with a full width at half maximum of 1–2 and 9–16 $\mathrm{nm}$. Due to low signal and poor sensitivity, nadir radiances above 2.0 $\mathrm{\mu }$m could not be used in this study. The SMART albedometer was calibrated in-laboratory, before and after the field campaign, with certified radiation sources (1000 $\mathrm{W}$ lamp, integrating sphere, and Lambert reflecting panel) prior and after the field campaign. In addition, in-field calibrations were performed during the campaign. Using Gaussian error propagation, the total uncertainty (one sigma) for the SMART albedometer ranged from 3 to 14 $\mathrm{\%}$ and depended on wavelength.

Concurrently, ground-based data of aerosol properties were collected in the city of Zhongshan with two sun photometers (Cimel CE318-2 and Schulz SP1A) and a LIDAR . From these data the single scattering albedo $\stackrel{\mathrm{̃}}{\mathit{\omega }}$ and asymmetry parameter $g$ of the aerosol particles were retrieved by the inversion method described by . The vertical aerosol extinction profile $b_{\text{ext}}\left(z\right)$ was measured with a combined Raman (387 $\mathrm{nm}$), elastic (355 and 532 $\mathrm{nm}$) backscatter LIDAR. An overview of the applied airborne and ground-based instruments and measurement uncertainties are given in Table  ().

Table of airborne and ground-based measurements in Zhongshan, China, 2009.

Instrument	Measured quantity	Uncertainty
SMART albedometer	$F_{\mathit{\lambda }}^{\uparrow }$, $I_{\mathit{\lambda }}^{\uparrow }$, 0.4–2.1 $\mathrm{\mu }$m	3–14 $\mathrm{\%}$
Sun photometer Schulz SP1A05	AOD	0.01
	0.35–1.0 $\mathrm{\mu }$m, 17 channels
Sun photometer CIMEL CE 318-2	AOD, $\stackrel{\mathrm{̃}}{\mathit{\omega }}$, $g$	AOD: 0.02–0.04
	0.34–1.0 $\mathrm{\mu }$m, 8 channels	$\stackrel{\mathrm{̃}}{\mathit{\omega }}$: 0.03
		$g$: 3–5 $\mathrm{\%}$
Combined Raman backscatter LIDAR	$b_{\text{ext}}\left(z\right)$	0.01–0.02 ${\mathrm{km}}^{-\mathrm{1}}$

The measured data from the sun photometers and LIDAR were used to derive the spectral vertical extinction profile, which is needed for the calculation of surface reflectance and surface albedo from the airborne radiance and irradiance measurements (see Sect. ). Since the LIDAR did not provide information on the vertical aerosol distribution up to 1000 $\mathrm{m}$ height, the vertical extinction profile was extrapolated from 1000 m to the surface in correspondence with the vertically integrated extinction coefficient (AOD) of the sun photometer. The layer between surface and 1000 $\mathrm{m}$ height was separated into two sub-layers. Starting from 1000 $\mathrm{m}$ altitude the extinction was assumed to increase with decreasing height following the measured slope of the LIDAR extinction profile. The second sub-layer with a constant aerosol extinction coefficient was assumed below to represent a well-mixed layer. The sun photometer measurements of AOD, $\stackrel{\mathrm{̃}}{\mathit{\omega }}$, $g$ were extrapolated to wavelengths larger than 1 $\mathrm{\mu }$m by scaling spectra from the literature . $\stackrel{\mathrm{̃}}{\mathit{\omega }}$ and $g$ were assumed to be constant with height. The shape of the vertical profile was assumed to have no spectral dependence, so that extinction in all wavelengths had the same relative profile. Thus, the spectral sun photometer AOD measurements could be scaled to provide extinction profiles at different wavelengths.

Retrieval of surface reflectance and albedo from airborne measurements

The spectral surface reflectance ${\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}$ was calculated according to the atmospheric correction technique, originally developed for the retrieval of the surface albedo ${\mathit{\alpha }}_{\mathrm{s},\mathit{\lambda }}$ : $${\mathit{\alpha }}_{\mathrm{s},\mathit{\lambda }}={\displaystyle \frac{F_{\mathit{\lambda }}^{\uparrow }\left(z_{\mathrm{0}}\right)}{F_{\mathit{\lambda }}^{\downarrow }\left(z_{\mathrm{0}}\right)}},$$ based on airborne measurements of irradiance. Note that used the symbol ${\mathit{\rho }}_{\mathrm{s}}$ for the surface albedo. Atmospheric masking means that a downward-looking radiation sensor measures not only the contribution of reflected radiation from the surface but is also influenced by the multiple scattering within the atmosphere, which needs to be separated. The surface reflectance $${\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}={\displaystyle \frac{\mathit{\pi }\cdot I_{\mathit{\lambda }}^{\uparrow }\left(z_{\mathrm{0}}\right)}{F_{\mathit{\lambda }}^{\downarrow }\left(z_{\mathrm{0}}\right)}}$$ is retrieved by assuming isotropic surface reflectance with $F=\mathit{\pi }\cdot I$. Both surface properties are obtained by an iterative algorithm starting with an initial value of albedo (reflectance) at flight level ${\mathit{\alpha }}_{\mathrm{m}}\left(z_{\mathrm{F}}\right)$ (${\mathit{\rho }}_{\mathrm{m}}\left(z_{\mathrm{F}}\right)$), which is used as a first guess for the surface albedo (reflectance). The spectral upward irradiance and radiance as well as the spectral downward irradiance at flight level $z_{\mathrm{F}}$ and at ground level $z_{\mathrm{0}}$ are calculated with the one-dimensional radiative transfer package libRadtran . The sun photometer and LIDAR measurements provide the input for the aerosol profile. Furthermore, the model uses atmospheric profiles of temperature, pressure, and gas concentrations. The simulated albedo (reflectance) at flight level and ground level is used to calculate the surface albedo (reflectance) for the next iteration step ${\mathit{\alpha }}_{n+\mathrm{1}}$ (${\mathit{\rho }}_{n+\mathrm{1}}$) as follows: $${\mathit{\alpha }}_{n+\mathrm{1}}\left(z_{\mathrm{0}}\right)={\mathit{\alpha }}_{\mathrm{m}}\left(z_{\mathrm{F}}\right)\cdot {\displaystyle \frac{{\mathit{\alpha }}_n\left(z_{\mathrm{0}}\right)}{{\mathit{\alpha }}_n\left(z_{\mathrm{F}}\right)}},$$ $${\mathit{\rho }}_{n+\mathrm{1}}\left(z_{\mathrm{0}}\right)={\mathit{\rho }}_{\mathrm{m}}\left(z_{\mathrm{F}}\right)\cdot {\displaystyle \frac{{\mathit{\rho }}_n\left(z_{\mathrm{0}}\right)}{{\mathit{\rho }}_n\left(z_{\mathrm{F}}\right)}}.$$ The iteration is stopped if the difference between measured $F_{\mathrm{m}}^{\uparrow }$ ($I_{\mathrm{m}}^{\uparrow }$) and modeled upward irradiance $F_{n+\mathrm{1}}^{\uparrow }$ (radiance $I_{n+\mathrm{1}}^{\uparrow }$) is less than 4 $\mathrm{\%}$ (6 $\mathrm{\%}$), which is comparable to the measurement uncertainty of the upward irradiance (radiance): $$\begin{array}{}& & \left|{\displaystyle \frac{F_{n+\mathrm{1}}^{\uparrow }\left(z_{\mathrm{F}}\right)}{F_{\mathrm{m}}^{\uparrow }\left(z_{\mathrm{F}}\right)}},-,\mathrm{1}\right|\le 0.04,& & \left|{\displaystyle \frac{I_{n+\mathrm{1}}^{\uparrow }\left(z_{\mathrm{F}}\right)}{I_{\mathrm{m}}^{\uparrow }\left(z_{\mathrm{F}}\right)}},-,\mathrm{1}\right|\le \mathrm{0.06\hspace{0.17em}.}\end{array}$$ The uncertainties of the retrieved surface albedo and reflectance were determined by Gaussian error propagation of the single error sources (see Table ). For this purpose, a sensitivity study with radiative transfer calculations was performed. Synthetic upward radiances and irradiances at a given flight altitude $z_{\mathrm{F}}$ were calculated with one-dimensional radiative transfer calculations for given aerosol and surface albedo input. Subsequently, the surface albedo and surface reflectance were determined with the method described above, first for the given input values and then with modified input parameters. The variations of the single aerosol parameters were AOD $\pm 30\mathrm{\%}$, $\stackrel{\mathrm{̃}}{\mathit{\omega }}\pm \mathrm{2}\mathrm{\%}$, $g\pm \mathrm{4}\mathrm{\%}$, and $z_{\mathrm{F}}\pm 50\mathrm{m}$. These uncertainties were assumed from the observed variation of the parameters during the field campaign in Zhongshan. The variation of $F_{\mathit{\lambda }}^{\uparrow }$ and $I_{\mathit{\lambda }}^{\uparrow }$ were $\pm (\mathrm{3}$–14) $\mathrm{\%}$ depending on the considered wavelength range. In general, lower AODs result in lower uncertainties of the derived surface properties because, compared to the surface, the relative contribution by the atmosphere to the measured upward radiation is decreased. Therefore, uncertainties of the aerosol properties $g$ and $\stackrel{\mathrm{̃}}{\mathit{\omega }}$ for lower AODs have lower impact on the retrieval results. Table  gives the sensitivities for AOD(0.53 $\mathrm{\mu }$m) $=0.9$ representing an upper limit. In wavelength range I (0.4–0.6 $\mathrm{\mu }$m), the retrieved surface albedo and surface reflectance are sensitive to the assumed optical aerosol. As AOD generally decreases with increasing wavelength, the sensitivity also decreases. At the longest wavelengths (1.8–2.1 $\mathrm{\mu }$m), the uncertainties are dominated by the measurement uncertainties of the SMART instrument, which has lower sensitivity and poor signal-to-noise ratio in this spectral range. A similar study was performed for AOD(0.53 $\mathrm{\mu }$m) $=0.12$ showing uncertainties between (3.2–10.8) $\mathrm{\%}$ (surface albedo) and (6.6–14.4) $\mathrm{\%}$ (surface reflectance).

Relative uncertainties of the different error sources and uncertainty (one sigma) of surface albedo and reflectance after Gaussian error propagation for high aerosol particle optical depth. The roman numerals I–IV indicate the wavelength ranges. I: 0.4–0.6 $\mathrm{\mu }$m, II: 0.6–1.0 $\mathrm{\mu }$m, III: 1.0–1.8 $\mathrm{\mu }$m, IV: 1.8–2.1 $\mathrm{\mu }$m.

Sensitivity (%)
	${\mathit{\alpha }}_{\mathrm{s}}\left(\mathit{\lambda }\right)$, high optical depth					${\mathit{\rho }}_{\mathrm{s}}\left(\mathit{\lambda }\right)$, high optical depth
	(AOD(0.53 $\mathrm{\mu }$m) $=$ 0.9)					(AOD(0.53 $\mathrm{\mu }$m) $=$ 0.9)
	I	II	III	IV		I	II	III	IV
$F_{\mathit{\lambda }}^{\uparrow }$	3.1	3.1	3.5	10.8	$I^{\uparrow }$	6.6	6.6	9.5	14.4
AOD	25	19	9	0.8	AOD	10	5	0.2	1.5
$\stackrel{\mathrm{̃}}{\mathit{\omega }}$	20	7	4.5	0.9	$\stackrel{\mathrm{̃}}{\mathit{\omega }}$	11	4	2.5	0.5
$g$	15	8	4.5	1.5	$g$	8.5	3	2	0.3
$z_{\mathrm{F}}$	1	0.1	0.1	0.1	$z_{\mathrm{F}}$	1	0.1	0.1	0.1
${\mathit{\alpha }}_{\mathrm{s},\mathit{\lambda }}$	35.5	22	11.6	11	${\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}$	18.4	9.7	10	14.5

Results from airborne measurements

Example spectra of surface reflectance and albedo

An 8 $\mathrm{km}$ by 4 $\mathrm{km}$ area of Zhongshan was systematically overflown on 3 December 2009, with 60 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ flight speed in different altitudes as shown in Fig. . The upward radiances and irradiances measured within the marked square are used to calculate the surface albedo and surface reflectance spectra. Measurements at 4100 $\mathrm{m}$ flight altitude correspond to a radiance footprint of 140 $\mathrm{m}$ across and 200–380 $\mathrm{m}$ along the flight direction for an integration time of the SMART albedometer between 1 and 4 sec with an radiance inlet opening angle of 2.1${}^{\circ }$. The solar zenith angle (SZA) ranged between 44 and 55${}^{\circ }$ during the 1.5 $\mathrm{h}$ of the flight. The change of SZA was considered in the derivation of the surface properties.

Flight pattern for the airborne measurements in Zhongshan on 3 December 2009. The area marked with the rectangle shows the position of the measurements which are averaged for the spatially averaged surface albedo and surface reflectance spectra.

[Figure omitted. See PDF]

Mean surface albedo ${\mathit{\alpha }}_{\mathrm{s},\mathit{\lambda }}$ and surface reflectance ${\mathit{\rho }}_{\mathrm{s},\mathit{\lambda }}$ measured over Zhongshan on 3 December 2009. The data were averaged for a time period of 40 $\min$ where the flight altitude was 4100 $\mathrm{m}$ and the SZA ranged between 45.7 and 48.7${}^{\circ }$. The standard deviation is displayed as vertical bars.

[Figure omitted. See PDF]

Figure  presents typical examples of spatially averaged surface albedo and reflectance spectra, derived from the airborne measurements in 4100 $\mathrm{m}$ altitude over Zhongshan. Note that the wavelength ranges with significant water vapor absorption are excluded. By assuming relative homogeneity of the atmosphere, the corresponding standard deviation (displayed as vertical bars) represents a measure of the heterogeneity of the underlying urban surface. Both spectra in Fig.  show spectral signatures typical for vegetation, including the “vegetation step” in the near-infrared range of 0.7–0.9 $\mathrm{\mu }$m. As a general feature, the surface albedo is lower than the surface reflectance for all wavelengths, which is due to the particular geometry of the measurement. The standard deviation of the spatially averaged surface albedo, derived from measurements at 4100 $\mathrm{m}$ flight altitude, is also smaller than that of the surface reflectance. Since albedo is an integral across all directions, its relative uncertainty is lower due to the directional averaging (or blurring). Consequently, the flight altitude and the aerosol conditions will have the most impact on the variability of the retrieved spatially averaged surface albedo from airborne measurements . Resolving the small-scale variability of surface albedo requires low-flying aircraft or ground-based measurements. The surface reflectance derived from radiance measurements at high flight altitudes reproduces the heterogeneity of the surface reflection caused by the narrow viewing angle of the radiance optical inlet.

Revised surface reflectance assumptions for the MODIS AOD retrieval method

According to Eqs. (2) and (3) the surface reflectance at 0.47 $\mathrm{\mu }$m wavelength in the MODIS AOD algorithm is derived from that at 0.65 $\mathrm{\mu }$m based on the surface reflectance at 2.1 $\mathrm{\mu }$m. The airborne measurements of the surface reflectance are not reliable for wavelengths larger than 2.0 $\mathrm{\mu }$m, due to the low sensitivity of the SMART albedometer. The relation between reflectances at 2.0 and 0.65 $\mathrm{\mu }$m cannot be used as approximation to modify Eq. (2) because of the strong differences between 2.0 and 2.1 $\mathrm{\mu }$m, and the finding that reflectances at 2.0 and 0.65 $\mathrm{\mu }$m are poorly correlated ($R=0.62$).

Therefore, our strategy is to modify only the blue/red fit parameter of Eq. (3) to match the airborne observations above Zhongshan. The surface reflectance at 0.47 $\mathrm{\mu }$m as a function of the surface reflectance at 0.65 $\mathrm{\mu }$m for the case study Zhongshan is shown in Fig. a. Linear regressions are performed in two ways, one with a free $y$ intercept (red solid line), the other with $y$ intercept fixed at 0 (slope coefficient only – black solid line). The two regressions are compared with the regression fit assumed for the operational MODIS algorithm (dashed line). All three equations are listed as Eqs. (14), (15), and (16). $$\begin{array}{cc}& & {\mathit{\rho }}_{\mathrm{s},0.47}=0.005+0.49\cdot {\mathit{\rho }}_{\mathrm{s},0.65}\text{MODIS operational}& {\mathit{\rho }}_{\mathrm{s},0.47}=0.02+0.63\cdot {\mathit{\rho }}_{\mathrm{s},0.65}\text{Zhongshan}\\ & & \mathrm{\text{regression with intercept}}& {\mathit{\rho }}_{\mathrm{s},0.47}=0.85\cdot {\mathit{\rho }}_{\mathrm{s},0.65}\text{Zhongshan}\\ & & \mathrm{\text{regression without intercept}}\end{array}$$ Clearly, the measured over-urban surface reflectance relationship for Zhongshan is much different than that assumed by the retrieval algorithm, with the slope coefficient being much larger (0.85 vs. 0.49). The histogram of the slope values is shown in Fig. b. A similar result, with higher values for the relationship between ${\mathit{\rho }}_{\mathrm{s},0.47}$ and ${\mathit{\rho }}_{\mathrm{s},0.65}$ of approximately 0.766 for urban surfaces is already reported by .

(a) Surface reflectance at 0.47 $\mathrm{\mu }$m ${\mathit{\rho }}_{\mathrm{s}}$(0.47 $\mathrm{\mu }$m) as function of surface reflectance at 0.65 $\mathrm{\mu }$m ${\mathit{\rho }}_{\mathrm{s}}$(0.65 $\mathrm{\mu }$m) for Zhongshan. The solid lines represent the linear fits of the measured data with (red) and without (black) an intercept. The dependency of ${\mathit{\rho }}_{\mathrm{s}}$(0.47 $\mathrm{\mu }$m) and ${\mathit{\rho }}_{\mathrm{s}}$(0.65 $\mathrm{\mu }$m) in the operational algorithm is shown as dashed line. (b) Relative frequency distribution for all possible ratios of ${\mathit{\rho }}_{\mathrm{s}}$(0.47 $\mathrm{\mu }$m) and ${\mathit{\rho }}_{\mathrm{s}}$(0.65 $\mathrm{\mu }$m) representing the slope $a_{0.47/0.65}$ for $y$ intercept at 0. The slope of the operational algorithm is shown as dashed line.

[Figure omitted. See PDF]

Implementation of measured surface reflectance ratios into the MODIS algorithm

SMART albedometer measurements of the spectral surface reflectance over Zhongshan are used to modify the surface reflection assumptions in the AOD retrievals with respect to the variability and the spectral signature of urban surfaces. The AOD will be retrieved from MODIS measurements using the modified and the operational algorithms for Zhongshan. Usually, the validation of AOD retrievals is performed by comparisons with AERONET measurements. Since no data were available for Zhongshan, the retrievals will be additionally applied for a larger MODIS data set (2010–2014) that covers the city of Beijing.

In a second approach, the MODIS retrieval methods will be applied for simulated TOA reflectances (representing MODIS measurements) for different surface types. The retrieval results from the operational and the modified algorithms will be compared to the reference AOD value which serves as input for the radiative transfer simulations.

Sensitivity study

The measured spectral reflectance at TOA, from the MODIS products MOD04 (from Terra satellite) and MYD04 (from Aqua satellite), was used as input for the MODIS retrieval stand-alone codes. These reflectances are cloud screened and are considered to be representative of 10 $\mathrm{km}$ by 10 $\mathrm{km}$ pixels (nadir). Note, that these data have already been processed by the operational retrieval (C5 product), and they are being reused to drive the stand-alone versions (C5 and C6). Therefore, the differences in the derived AOD are solely caused by the modification of the stand-alone code, not by differences in the preprocessing of the data between C5 and C6. Here, the main difference between C5 and C6 stand-alone versions is the treatment of the spectral slope $a_{0.65/2.1}$.

The MODIS case was chosen to fit the time and area of the airborne measurements as best as possible in order to minimize differences due to season and viewing geometry. The MODIS data granule from 1 December 2009 was selected as optimum choice with respect to observed location and solar geometry (SZA $=$ 45.7${}^{\circ }$). Although optimal geometry for MODIS retrieval, there was no comparable sun photometer measurement available for this day. However, the sensitivity of the algorithm can be tested for this case. Even though the observed spectral slope, $a_{0.47/0.65}$, measured in Zhongshan, has mean value close to 0.85, Fig. b demonstrates that there is a wide range for measured spectral slope (0.55 to 1.10). This range of slope values is applied to each of the modified algorithm (C5 and C6 based) to determine the effect on the retrieved AOD.

Figure a shows the retrieved AOD (0.55 $\mathrm{\mu }$m) as function of the spectral slope $a_{0.47/0.65}$ measured in Zhongshan for both of the modified retrievals. For this range of spectral slopes, an almost linear relationship between slope and retrieved AOD is observed, with a correlation coefficient of $R=0.99$ for both algorithms. For this particular case, varying the spectral slope between 0.55 and 1.10 results in AOD ranging between 0.87 and 0.75. C6 tends to show slightly more sensitivity to increasing the slope of the surface reflectance, resulting in lower AOD for equivalent slope assumptions. Figure b shows the relative frequency distribution of AOD (0.55 $\mathrm{\mu }$m) as retrieved with the modified C5 and C6 algorithms using different slope parameters $a_{0.47/0.65}$ (Fig. ). The AOD values of 0.84 (C5) and 0.85 (C6) as retrieved by the original MODIS algorithms are displayed as vertical lines. In most runs, the modified algorithms yielded AOD values which were lower than the operational MODIS algorithms. Clearly, modifying the MODIS retrieval code to use measured spectral slopes, leads to revised estimates of the AOD. Since the majority of measured spectral slopes are larger than assumed for the operational retrieval, the result is a lower estimate of AOD. For this case, the sensitivity of the C6 version to assumed blue/red spectral slope is greater than it was for C5.

Testing over Beijing

The previous subsection showed a sensitivity of aerosol retrieval to changing the assumed spectral slope of the surface properties. Now it is tested whether the measured surface ratios actually improve the aerosol retrieval. Since there are such few collocations of AERONET data with MODIS overpass in Zhongshan, it is considered to test the urban surface ratios for use over another city in China. Here the spectral slopes $a_{0.47/0.65}$, derived from the airborne measurements over Zhongshan, were applied to MODIS data over Beijing. For Beijing, there are 120 days of concurrent AERONET data, observed between 2010 and 2014 matching in time ($\pm$30 min) and space (closest MODIS pixel to AERONET station).

(a) Correlation between slope $a_{0.47/0.65}$ used for the AOD retrieval and the resulting AOD at 0.55 $\mathrm{\mu }$m for Zhongshan for the C5 and the C6 algorithm applied for measurements on 1 December 2009. (b) Relative frequency distribution of retrieved AOD for Zhongshan for the modified C5 and C6 code. The corresponding AOD of the unchanged operational retrievals is shown as vertical lines.

[Figure omitted. See PDF]

Figure  presents the scatter plot and regression lines of the relation between AOD derived from MODIS vs. measured by AERONET. Recall that when the retrieval expects a dark surface, a larger observed TOA reflectance is interpreted as enhanced aerosol contribution. Both operational algorithms systematically overestimate the AOD as compared to AERONET, which is consistent with urban surfaces being brighter than expected. However, the C6 version is less biased. Since the operational C6 equations include a larger slope parameter $a_{0.65/2.1}^{\text{NDVI}}$ than C5, the estimated surface reflectance is larger, leading to reduced retrievals of AOD.

Using the measured $a_{0.47/0.65}$-spectral slopes in place of either C5 and C6 assumed slopes will further reduce retrieved AOD. Both modified methods show a reduction of mean differences between AOD(AERONET) and AOD(MODIS) as displayed in Table . However, while the mean bias is reduced in either version, modifying the measured slopes also leads to retrieval of negative AOD values, as shown in Fig. . Modifying the C6 version, which already had reduced bias compared to C5 operational version, tends to lead to even more retrievals of negative AOD. At the same time, the AOD standard deviations have increased (from $\pm 0.12$ to $\pm 0.17$) when using the modified C6 algorithm. One of the reasons might be the different sensitivity of TOA reflectance on changes of AOD over bright and dark surfaces. Over dark surfaces (as assumed in the operational retrieval) changes of AOD result in larger changes of TOA reflectance than over bright urban surfaces (as assumed in the modified retrieval). Consequently, measurement uncertainties have a larger effect on the retrieval uncertainty for an urban surface.

Scatter plot of AOD retrieved from MODIS data using (a) the operational and (b)  modified Collections 5 and 6 method and the AERONET data for Beijing. The dashed line marks the $\mathrm{1}:\mathrm{1}$ relation.

[Figure omitted. See PDF]

Distribution of AOD over Beijing area on 17 May 2012, derived from MODIS measurements by (a) operational C5 method, (b) modified C5 method, (c) operational C6 method, and (d) modified C6 method. Corresponding AOD from AERONET (Beijing, white cross) was 0.55.

[Figure omitted. See PDF]

Mean value and standard deviation of differences between AOD(MODIS) and AOD(AERONET).

Retrieval method	Mean difference
	$\pm$ standard deviation
C5 modified	0.10 $\pm$ 0.15
C5 operational	0.27 $\pm$ 0.14
C6 modified	$-$0.02 $\pm$ 0.17
C6 operational	0.19 $\pm$ 0.12

and have applied not only revised assumptions of the surface reflection based on the MODIS surface reflectance product (MOD09) but also changes in the aerosol model that uses AERONET AOD data to compute the values of $\stackrel{\mathrm{̃}}{\mathit{\omega }}$ and $g$. achieved a correlation with AERONET data with a coefficient of 0.95 and a relative mean bias of 0.97 (AOD(retrieved)$/$AOD(AERONET)). For the modified C5 (C6) algorithm a correlation coefficient of 0.95 (0.93) and a relative mean bias of 1.25 (0.95) was derived (cf. operational C5 (C6): correlation coefficient $=$ 0.95 (0.95), relative mean bias $=$ 1.77 (1.55)). In Fig.  only the mean values of the retrieved AOD distributions derived from the modified methods are used. The spread of the AOD distributions depend significantly on the AOD, higher values of AOD lead to a smaller spread of retrieved AODs (standard deviation $=$ 0.03 for AOD $=$ 1.6 vs. standard deviation $=$ 0.18 for AOD $=$ 0.2) which demonstrates the decreasing sensitivity of the spectral surface reflection assumptions on the retrieval uncertainty with increasing AOD as already shown by and .

The areal distribution of the AOD (colored open squares) over Beijing (17 May 2012) retrieved from MODIS measurements by the operational and modified C5 and C6 methods are presented in Fig. . The AERONET station has measured an AOD of 0.55 for the time of the satellite overpass. The operational algorithms C5 and C6 have shown AODs of about 0.77 (C5) and 0.72 (C6) over Beijing (39.977${}^{\circ }$ N, 116.381${}^{\circ }$ E; white cross). By applying the modified algorithms the AODs over Beijing decreased to values of 0.45 (C5) and 0.23 (C6). For rural areas the modified methods fail to retrieve the AOD. There, negative values are derived indicated by the dark blue squares. However, for urban areas the modified algorithm C5 improves the agreement between AERONET and MODIS AOD.

Sensitivity study of the visible-to-near-infrared reflectance ratios

According to Eq. (2) the MODIS AOD retrieval considers the surface reflectance at 2.1 $\mathrm{\mu }$m in its surface reflectance parameterization. However, 2.1 $\mathrm{\mu }$m is not well determined by airborne measurements due to decreased sensitivity for wavelengths larger than 2.0 $\mathrm{\mu }$m. In order to account for the surface reflectance at 2.1 $\mathrm{\mu }$m and its relation to that at 0.65 $\mathrm{\mu }$m a sensitivity study based on radiative transfer simulations was performed. This study was expanded also to surface types other than urban surfaces to investigate the limitations of the modified C5 and C6 algorithms. Simulated reflectances at TOA for variable surface properties (surface albedo) have been used to generate synthetic MODIS measurements providing the input for the AOD retrieval. The spectral slope and $y$ intercept for the calculation of ${\mathit{\rho }}_{0.65}$ and ${\mathit{\rho }}_{0.47}$ (Eq. ) in the stand-alone code of the aerosol retrieval over land were modified according to the spectra of the different surface types used for the simulation of the TOA reflectances. Figure a shows aircraft-measured surface albedo spectra for four different land surface types, within the wavelength range of 0.4–2.1 $\mathrm{\mu }$m. The four types include our urban measurements taken over Zhongshan, plus previous measurements of Kentucky bluegrass (), oasis, and plowed farmland (). The gaps indicate strong gas absorption bands; they are excluded when deducing the surface albedo from airborne measurements following the procedure described in Sect. 3. A strong increase of surface albedo at 0.7 $\mathrm{\mu }$m due to vegetation is typical for the surface types Kentucky bluegrass (gray dots) and oasis (red dots). Also for the urban surface the impact of vegetation is visible but less prominent in comparison to Kentucky bluegrass and oasis. The surface albedo spectrum of farmland shows the highest albedo values for wavelengths below 0.7 $\mathrm{\mu }$m. The spectral reflectances simulated for TOA as used as synthetic MODIS measurements for the four cases are following the spectral surface albedo with largest differences between farmland (${\mathit{\rho }}_{0.65}$(TOA) $=0.21$) and Kentucky bluegrass (${\mathit{\rho }}_{0.65}$(TOA) $=0.08$). The spectral characteristics of the optical properties of the aerosol particles, as used in the simulation of reflectances at TOA, are presented in Fig. c. For a wavelength of 0.55 $\mathrm{\mu }$m the AOD was set to 0.85 which is comparable to the aerosol conditions during the measurement campaign in Zhongshan. Values of $\stackrel{\mathrm{̃}}{\mathit{\omega }}$(0.55 $\mathrm{\mu }$m) $=0.9$ and $g$(0.55 $\mathrm{\mu }$m) $=0.7$ were applied according to the sun photometer measurements in Zhongshan. The aerosol particles are more absorbing than the moderately absorbing aerosol type ($\stackrel{\mathrm{̃}}{\mathit{\omega }}$(0.55 $\mathrm{\mu }$m) $=0.92$) in the MODIS algorithm but less absorbing than the absorbing aerosol type ($\stackrel{\mathrm{̃}}{\mathit{\omega }}$(0.55 $\mathrm{\mu }$m) $=0.87$). The asymmetry parameter used for the simulations in this study compares with the spheroid dust model ($g$(0.55 $\mathrm{\mu }$m) $=0.7$) and the weak absorbing aerosol type ($g$(0.55 $\mathrm{\mu }$m) $=0.68$).

(a) Surface albedo spectra of different surface types. For urban, oasis, and farmland the measurement uncertainty is given as vertical error bars. (b) Spectral characteristic of the optical aerosol properties used in the model study.

[Figure omitted. See PDF]

First the AOD was retrieved with the operational algorithms C5 and C6 for each of the four surface type cases. Then the spectral slopes $a_{0.65/2.1}$ and $a_{0.47/0.65}$ were modified for each surface type (a unique pair of slopes for each). The $y$ intercepts $b_{0.65/2.1}$ and $b_{0.47/0.65}$ were set to 0, because only one spectrum for each surface type was used. An overview of the slopes and $y$ intercepts used by the operational algorithm and those obtained from the example spectra are presented in Table . The AOD results of the four algorithms for all surface types are listed in Table . Best agreement of the operational algorithm C5 (C6) with the input value of AOD(0.55 $\mathrm{\mu }$m) $=$ 0.85 was found for underlying dark vegetated surfaces with AOD (0.55 $\mathrm{\mu }$m) $=$ 0.87 (0.88) for the case of Kentucky bluegrass and AOD(0.55 $\mathrm{\mu }$m) $=$ 0.85 (0.86) for the oasis surface type. Highest discrepancies were observed for the results from the operational algorithms for the farmland and the urban surface type.

Overview of the slopes and intercepts for the calculation of the surface reflectance ${\mathit{\rho }}_{\mathrm{s},0.46}$ and ${\mathit{\rho }}_{\mathrm{s},0.47}$ following Eqs. (1) and (2) for different surface types used in the operational and modified algorithms C5 and C6. The C6 parameters are given in brackets when they differ from the C5 approach.

	Operational algorithm C5 (C6)				Modified algorithm
Surface type	$a_{0.47/0.65}$	$b_{0.47/0.65}$	$a_{0.65/2.1}$	$b_{0.65/2.1}$	$a_{0.47/0.65}$	$a_{0.65/2.1}$
Kentucky bluegrass	0.49	0.005	0.58 (0.48)	$-$0.00143	1	0.3
Farmland	0.49	0.005	0.48 (0.58)	$-$0.00143	0.43	0.68
Oasis	0.49	0.005	0.51 (0.56)	$-$0.00143	0.49	0.81
Urban	0.49	0.005	0.49 (0.59)	$-$0.00143	0.8	0.8

Retrieved AOD(0.55 $\mathrm{\mu }$m) for different surface types and different slope parameters as listed in Table  for the C5 and C6 algorithms. The reflectance at TOA, input parameter in the retrieval algorithm, was simulated for AOD(0.55 $\mathrm{\mu }$m) $=$ 0.85.

	Operational algorithm		Modified algorithm
Surface type	C5	C6	C5	C6
Kentucky bluegrass	0.87	0.88	0.85	0.84
Farmland	0.71	0.61	0.61	0.64
Oasis	0.85	0.86	0.78	0.82
Urban	0.99	1.01	0.86	0.90

While for the urban surface an overestimation of the AOD was observed with a relative mean bias of 16 and 19 $\mathrm{\%}$, which confirms the outcome of the Beijing study reported in Sect. , the retrieved AOD for the farmland surface type is significantly underestimated (C6 with RMB $=$ 70 $\mathrm{\%}$). Both algorithms C5 and C6 include a threshold for their brightness criteria (2.1 $\mathrm{\mu }$m reflectance $>0.25$) which is exceeded for the farmland case (2.1 $\mathrm{\mu }$m reflectance $=$ 0.26). As stated by the retrieval is error-prone for cases above the brightness criteria. This underestimation can also be found in the modified algorithms with surface modified slopes. This also indicates the increased uncertainty for data which exceed the brightness criteria. For all other surfaces (Kentucky bluegrass, oasis, and urban) the modified algorithms exhibit good performances.

Conclusions

Operational MODIS aerosol retrieval performs poorly over urban surfaces, providing a high-biased estimate of AOD over cities such as Beijing, China. The retrieved AOD significantly depends on surface reflectance assumptions within the retrieval algorithms. Two versions of the operational algorithm (C5 and C6) have different assumptions about the surface reflectance. While both versions tend to show high bias over cities, the C6 version tends to be less so. Here, it is tested whether the aerosol retrieval can be improved by using alternative surface reflectance assumptions for urban surfaces.

First, surface reflectance assumptions were modified in the standard algorithms C5 and C6 using airborne measurements of surface reflectance over Zhongshan, China. These assumptions consider spectral relationships between the surface reflectance at 0.47 $\mathrm{\mu }$m and at 0.65 $\mathrm{\mu }$m. The modified retrieval methods have been applied to MODIS measurements taken during the aircraft campaign in the Zhongshan area. Additionally, the AOD retrieved by the standard and modified algorithms based on MODIS measurements from 2010 to 2014 over Beijing have been compared with AOD measurements of AERONET. In a second approach the sensitivity of the AOD retrieval for several surface types was studied using synthetic MODIS data obtained from radiative transfer simulations.

Airborne measurements of the spectral nadir radiance have been used to retrieve the spectral surface reflectance by aerosol properties obtained from sun photometer and LIDAR of the Zhongshan area. The fit parameters of the linear relationship between the surface reflectance at 0.47 $\mathrm{\mu }$m and at 0.65 in C5 and C6 were varied based on the variability of the measured reflectances over Zhongshan. Due to a lack of sensitivity and a low raw signal above of 2.0 $\mathrm{\mu }$m wavelength of the airborne instrument SMART albedometer (radiance component) only the $a_{0.47/0.65}$ slopes were modified in C5 and C6. The distribution of $a_{0.47/0.65}$ slopes over Zhongshan has revealed a mean value of $0.85\pm 0.09$. This indicates that over urban areas, there is much less spectral dependence of the surface reflectance than is assumed in the operational algorithm where ($a_{0.47/0.65}=0.49$).

For the test case over Zhongshan, the AOD values retrieved by the standard algorithms (AOD $=$ 0.84 (C5) and AOD $=$ 0.85 (C6)) agree with those retrieved by the modified algorithms ($0.83\pm 0.02$ (C5) and $0.82\pm 0.02$ (C6)) within the range of measurement uncertainty. In this test case, presumably the AOD is large enough that the impact of the surface reflectance relationship is not a significant factor in the retrieval.

Since there were no AERONET data available for Zhongshan, the surface reflectance relationship measured over Zhongshan was modified to be used for retrieval of MODIS measurements over Beijing. The modified algorithms were run for MODIS data over Beijing for the period of 2010 to 2014. While the operational algorithms C5 and C6 have shown a significant positive bias compared to AERONET with a mean difference of 0.27 and 0.19, the modified algorithms revealed a reduced mean difference of 0.10 and $-$0.02. While reducing the mean bias over the urban areas, using the urban surface reflection assumptions over non-urban surfaces results in retrieving negative AOD values. Therefore, using assumptions of urban surface spectral reflectance is only applicable for urban areas.

The sensitivity study, based on radiative transfer simulations of the reflectance at TOA (representing MODIS measurements) for given conditions as surface and aerosol properties, confirmed that both operational MODIS aerosol retrievals (C5 and C6) over land reproduce the input AOD(0.55 $\mathrm{\mu }$m) of 0.85 for the surface types Kentucky bluegrass (AOD(0.55 $\mathrm{\mu }$m) $=$ 0.87 (C5) and AOD(0.55 $\mathrm{\mu }$m) $=$ 0.88 (C6)) and oasis (AOD(0.55 $\mathrm{\mu }$m) $=$ 0.85 (C5) and AOD(0.55 $\mathrm{\mu }$m) $=$ 0.86 (C6)). When using modified spectral slopes $a_{0.65/2.1}$ and $a_{0.47/0.65}$ in the C5 algorithm a good performance of the aerosol retrieval was also obtained for the urban surface type with AOD(0.55 $\mathrm{\mu }$m) $=$ 0.86, whereby the modified C6 algorithm tends to overestimate the AOD with AOD(0.55 $\mathrm{\mu }$m) $=$ 0.90. Discrepancies were found for the surface type farmland where the brightness criterium in the MODIS retrieval is exceeded. For all methods (operational and modified) over the farmland surface, the retrieved AOD (0.61–0.71) was significantly lower than the reference value of 0.85.

The operational MODIS aerosol algorithms of collections 5 and 6 retrieve the AOD within the expected uncertainty for densely vegetated surfaces, but shortcomings can be observed for higher reflecting urban areas and farmland. The increased retrieval uncertainty for urban areas can be overcome by using spectral slopes suitable for urban surface within the algorithm. Changing the spectral slope parameter for urban surfaces improved the agreement between AOD(MODIS) and AOD(AERONET). The “stand-alone code” of the MODIS aerosol retrieval algorithm is freely available for the scientific community. For further work it is suggested to change the surface assumption parameters as proposed in this work and to use the modified retrieval algorithm for similar analyses in other urban areas. It is worthwhile to investigate if the surface assumption parameter derived from measurements over Zhongshan also delivers improved aerosol retrievals in other regions than Asia. In general other aircraft or surface-based measurements of the spectral surface reflectance in urban areas can be used to modify the surface assumption parameter in the retrieval code similar to this work. Further studies will validate the MODIS surface reflectance product (MOD09) against aircraft measurements to test whether MOD09 products can be used for the modification of the surface assumptions in the operational aerosol algorithm.

The spectral slopes used in the modified algorithms were derived from airborne nadir radiance measurements. For future investigations airborne imaging spectrometer measurements are suggested to characterize the slope parameters for urban areas also for different viewing geometries than the nadir observation. This will improve the corrections of the operational retrievals due to effect of the bidirectional reflectance distribution function (BRDF), since the anisotropy reflection of the surface is better described by the BRDF than assuming a Lambertian surface .

Acknowledgements

This work was conducted within the scope of the German Research Foundation (DFG) priority program SPP 1233 “Megacities Megachallenge – Informal Dynamics of Global Change” (WE 1900/16-2). This research was funded by the German Research Foundation (DFG, JA2023/2-2). We thank our cooperation partners of the Institute of Remote Sensing and Digital Earth for providing the possibility to join their airborne field experiment in Zhongshan. Furthermore, we thank Hong-Bin Chen and Philippe Goloub for their effort in establishing and maintaining the Beijing AERONET site. Edited by: R. Schofield