A surface reflectance scheme for retrieving aerosol optical depth over urban surfaces in MODIS Dark Target retrieval algorithm

The MODerate resolution Imaging Spectroradiometer (MODIS) instruments, aboard the two Earth Observing System (EOS) satellites Terra and Aqua, provide aerosol information with nearly daily global coverage at moderate spatial resolution (10 and 3 km). Almost 15 years of aerosol data records are now available from MODIS that can be used for various climate and air-quality applications. However, the application of MODIS aerosol products for airquality concerns is limited by a reduction in retrieval accuracy over urban surfaces. This is largely because the urban surface reflectance behaves differently than that assumed for natural surfaces. In this study, we address the inaccuracies produced by the MODIS Dark Target (MDT) algorithm aerosol optical depth (AOD) retrievals over urban areas and suggest improvements by modifying the surface reflectance scheme in the algorithm. By integrating MODIS Land Surface Reflectance and Land Cover Type information into the aerosol surface parameterization scheme for urban areas, much of the issues associated with the standard algorithm have been mitigated for our test region, the continental United States (CONUS). The new surface scheme takes into account the change in underlying surface type and is only applied for MODIS pixels with urban percentage (UP) larger than 20 %. Over the urban areas where the new scheme has been applied (UP > 20 %), the number of AOD retrievals falling within expected error (EE %) has increased by 20 %, and the strong positive bias against ground-based sun photometry has been eliminated. However, we note that the new retrieval introduces a small negative bias for AOD values less than 0.1 due to the ultra-sensitivity of the AOD retrieval to the surface parameterization under low atmospheric aerosol loadings. Global application of the new urban surface parameterization appears promising, but further research and analysis are required before global implementation.

retrievals) include the TOA spectral reflectance used for each retrieved AOD value, the data 1 provided in the C6 output file (e.g. MxD04_L2) can be recycled through the S-MDT 2 algorithm to retrieve the same AOD value as provided within the C6 product. Thus, the S-3 MDT can be easily modified to test different assumptions within the retrieval, including 4 surface reflectance assumptions. The insights gleaned from the S-MDT exercises, can then be 5 transferred back to modify the full (operational) retrieval code, and tests (and statistics) can be 6 performed for global data. The land surface is too variable to apply an explicit model to describe its spectral optical 10 properties. Thus MDT uses an empirical parameterization based on only three MODIS bands.
11 Kaufman et al., (1997aKaufman et al., ( , 2002 noted that for most vegetated and dark-soiled land surfaces, 12 observations showed that surface reflectance in a blue wavelength (0.47 µm) and red 13 wavelength (0.65 µm) were about ¼ and ½ of the surface reflectance in a shortwave infrared surfaces. As a land surface transitions from natural vegetation to manmade structures such as 1 buildings and roads, the global SWIR-VIS relationship is violated. This was noted in Levy et 2 al., (2007b, 2009). In addition, Castanho et al., 2007 analyzed the SWIR-VIS relationship 3 over Mexico City and found that the ratios of SWIR-red are much higher (0.73-0.76) than 4 what is assumed in the MDT global algorithm (0.50 -0.55). They also found that the SWIR-5 VIS relationship strongly depends on differences in urbanization fraction. Using the global 6 values for surface ratio when the actual ratio is much higher will underestimate the visible 7 surface reflectance and the resulting aerosol contribution will be overestimated. Munchak et around each AERONET site. While, there have been attempts to improve the MDT aerosol 12 retrieval on the local or even regional scale, no attempt has been made to date to develop a 13 general improved surface reflectance parameterization that can be applied to the global 14 retrieval. In this study, we develop a new SWIR-VIS scheme that accounts for NDVI SWIR and 15 scattering angle, but also accounts for the UP. 3 Data and Study Region 18 We explore whether MDT surface parameterization can be modified to account for urban 19 surface properties. We start with the research S-MDT version (Levy and Pinker, 2007;and 20 http://darktarget.gsfc.nasa.gov), and develop a replacement SWIR-VIS parameterization that 21 includes dependence on UP. To develop the replacement SWIR-VIS parameterization, we 22 rely on two datasets. These are the MODIS Land Cover Type, and the MODIS Land Surface 23 Reflectance Product. The MDT algorithm with the replacement SWIR-VIS parameterization, 24 we denote as the or C6U version. 25 To test the C6U retrieval, we compare the retrieved AOD with the sunphotometer 26 measurements taken at standard AERONET sites, as well as AERONET's DRAGONs during Atmos. Meas. Tech. Discuss., doi:10.5194/amt-2015-375, 2016 Manuscript under review for journal Atmos.  Figure 1 shows the map of UP for the CONUS region. This map of UP has been 2 used in the C6U algorithm to define and apply the new surface scheme for each MODIS 10km 3 AOD retrieval. The surface parameterization in the MDT algorithm was based on performing atmospheric 7 correction of MODIS-measured reflectance near AERONET sites, which led to formulation of 8 the empirical SWIR-VIS relationships (e.g., Levy et al., 2007b;2013). Rather than repeating 9 this exercise, we rely on the MODIS-derived land spectral surface reflectance product because 10 it is available and has been validated with similar atmospheric correction exercises. (1) Reduce the known bias between MDT retrievals and sunphotometer measurements in 2 urban areas;

3
(2) No degradation of retrieval quality in non-urban areas;

4
(3) Must be an operational way to identify when to apply the new scheme and when not;

5
(4) It must slip into the structure of the standard operational algorithm without requiring 6 extensive additional modifications to the existing algorithm; 7 The first step to meeting these requirements is to develop surface parameterizations 8 characteristics of urban surfaces. To do this we use the MODIS land product as described in 9 Section 3. The main advantage of using this product is its spatial coverage.  (shown as gray color), but bin the data to help visualize the differences from regime to 5 regime. We note that the slopes of the regression between the blue and red wavelengths are 6 not strongly dependent on differences in the UP or the NDVI SWIR , as long as UP > 20%.  Table 1. This provides a quantitative indicator of 13 uncertainties in surface reflectance of visible bands in the new surface scheme.

14
The regression statistics provide the information needed to relate the surface reflectance (ρ s ) 15 at 0.65 μm and 0.47 μm to that at 2.12 µm. We take into account the effects of geometry by As shown in Table 1, the new slopes of SWIR-VIS over urban areas are significantly higher 1 than those assumed in the C6 algorithm for global application, which is consistent with other 2 studies (i.e. Levy et al., 2007b, Castanho et al., 2007, Min et al., 2010. The red-SWIR ratios 3 are also higher for less vegetated areas (NDVI SWIR <0.2) than those with more vegetation 4 areas (NDVI SWIR > 0.2).

5
The new surface scheme is expected to provide better surface reflectance estimates over urban 6 areas than the existing operational scheme (i.e. C6). Figure

13
UP is the only new parameter added to the surface scheme. This is a globally available 14 parameter calculated from a standard MODIS Land product, and is updated annually. The

15
C6U algorithm inputs UP and uses this parameter to decide whether to apply the old or the 16 new surface parameterization. Dependence on UP assures that there will be no changes to the 17 retrieved aerosol products for non-urban surfaces (UP < 20%).

18
We have developed this surface scheme using data from the continental U.S., expecting the 19 scheme to be optimized for this region. Global implementation may be challenging and we 20 discuss these challenges in Section 6.
Expected error (Remer et al., 2005) for AOD retrieved over land, as given by the MODIS 18 science team through validation exercises is presented in equation 6.
where the AOD AERO is AOD from AERONET, and AOD MODIS is the retrieved value from 22 MODIS (either C6 or C6U version). Note the AOD is reported at 0.55 µm, which means that 23 interpolation (quadratic on log-log scale) is performed for the AERONET data.

24
In addition, we have computed linear regression statistics, including correlation coefficient 25 (R), regression coefficients (slope and intercepts) and number of data points (N). as compared to SP AODs. The bias is significantly reduced (to -0.01) when the C6U retrieval 8 is applied (right panel).

9
As discussed in section 4, aerosol retrieval from passive satellite measurements is sensitive to 10 underlying surface reflectance, and the relative uncertainty becomes greater for lower AOD 11 conditions, when the surface reflectance dominates the signal. By improving the overall 12 urban surface reflectance parameterization, the bias is reduced. However, the small negative 13 bias indicates that there is still uncertainty in estimating visible reflectance. We note that most error envelope (EE%) has also increased from 63 % for C6 to 85% for the C6U retrievals.

19
Most of the statistical parameters in comparing MODIS and SP AODs demonstrate 20 improvement in C6U AODs as compared to C6 retrievals.

11
In order to analyze and validate MODIS AODs with lower quality flags, we have grouped the 12 data into several different ways to represent different retrieval conditions. For the third category (UP>20%) this includes only retrievals where C6U retrieval would be 19 applied and different from C6, while the 2 nd category includes retrievals that may be suburbs 20 or small towns, and the 1 st category (ALL) includes everything.

21
Each surface category in Table 2 is further broken down by QAF level. Note that the case of 22 ALL and QAF=3 represents the data in Fig. 4, which are the MDT recommendations of 23 "science quality" data. The correlation for C6U (R=0.86) is marginally higher than for C6 24 (R=0.85), but the bias is reduced to 0.006 from 0.022. The number of retrievals within 25 expected error (EE%) increased by 4%. For the UP>20% cases (e.g. data shown in Fig. 3), 26 there is a huge reduction in bias (from 0.075 to -0.007), and the EE% increases from 58.6% to 27 85.3%.

28
For the C6 data, there is a clear reduction in regression quality (decreased correlation, 29 increased bias, %EE reduction) as QAF criteria are relaxed from 3 to 2 to 1. This is true for window (EE%) at CCNY almost doubled from 46% in C6 to 83% in C6U, and the bias is 10 reduced to zero from 0.09. The C6U results over CCNY observed a slight negative offset of -11 0.03 mainly due to a negative bias in low aerosol loading conditions (AOD<0.1). In general 12 the slope is larger than one for eastern US sites whereas it less than one for western US sites.

13
In Figure 6 we have evaluated the spatial distribution of AOD over the region covering two surrounding suburban and rural areas on the map in the C6 retrievals (Fig 6a). As noted in C6U AOD retrievals can be very helpful in studying these small scales features over large 10 urban centers and their connection with urban pollution and population. Figure 5c shows the 11 difference between C6 and C6U AODs, which is correlated with UP (Fig 5d) and could be as  AERONET is reduced to almost zero, and EE% has gone up to 93%. This implies that the 10 RMSE and slope have also improved in the C6U retrievals, as compared with the C6 11 retrievals, but the offset (I) in C6U is now slightly negative (-0.02), which was zero in the C6 12 retrievals. A closer look at the scatter plot reveals that C6U has introduced some negative bias 13 at very low AOD cases (AOD<0.15), which is consistent with previous analysis over the 14 CONUS region. Again, in order to remove these negative biases; more accurate estimation of encouraging, we decided to run the CONUS-derived C6U algorithm globally and to compare 1 the results with AERONET measurements. Figure 8 presents the difference between C6U 2 AOD with AERONET AOD as a function of UP over the entire global data set, excluding 3 AERONET stations from the CONUS region. The CONUS region has been excluded in order 4 to avoid weighting the results by the formulation data set. The results are surprisingly good.

5
The CONUS-derived C6U has reduced the positive bias over urban surfaces to nearly zero, 6 globally. This analysis is very encouraging, but more in depth analysis for specific 7 stations/regions will be required to better understand the new algorithm before applying it 8 operationally at global scale.

23
Additional testing will be necessary first.

24
As populations flock to urban centers causing the urban landscape around the world to grow  Urban % >20: MODIS pixels with urban percentage larger than 20% selected in collocation.

5
Statistical parameters, number of data points (N) linear correlation coefficient (R), mean bias 6 (Bias), and EE% are reported in this