An improved method for retrieving nighttime aerosol optical thickness from the VIIRS Day / Night Band

Using Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) data, a method, dubbed the “variance method”, is developed for retrieving nighttime aerosol optical thickness (τ ) values through the examination of the dispersion of radiance values above an artificial light source. Based on the improvement of a previous algorithm, this updated method derives a semi-quantitative indicator of nighttime τ using artificial light sources. Nighttime τ retrievals from the newly developed method are inter-compared with an interpolated value from late afternoon and early morning ground observations from four AErosol RObotic NETwork (AERONET) sites as well as column-integrated τ from one High Spectral Resolution Lidar (HSRL) site at Huntsville, AL, during the NASA Studies of Emissions and Atmospheric Composition, Clouds and Climate Coupling by Regional Surveys (SEACRS) campaign, providing full diel coverage. Sensitivity studies are performed to examine the effects of lunar illumination on VIIRS τ retrievals made via the variance method, revealing that lunar contamination may have a smaller impact than previously thought; however, the small sample size of this study limits the conclusiveness thus far. VIIRS τ retrievals yield a coefficient of determination (r) of 0.60 and a root-meansquared error (RMSE) of 0.18 when compared against straddling daytime-averaged AERONET τ values. Preliminary results suggest that artificial light sources can be used for estimating regional and global nighttime aerosol distributions in the future.


Introduction
Recently, observationally based nighttime aerosol studies have drawn increasing attention from the research community (e.g., Zhang et al., 2008b;Berkoff et al., 2011;Lee and Sohn, 2012;Campbell et al., 2012;Johnson et al., 2013).Berkoff et al. (2011) and Barreto et al. (2013) present surface-based methods to estimate nighttime aerosol optical properties through measuring the attenuated moonlight using a modified version of sun photometers.Based on Berkoff et al. (2011), the NASA AErosol RObotic NETwork (AERONET, Holben et al., 1998) team is developing an operational nighttime aerosol retrieval capability (B. Holben, personal communication, 2015).While mineral dust can be detected and in part quantified from space by infrared means (e.g., Miller, 2003;Hansell et al., 2007;Lee and Sohn, 2012), the recently launched Visible Infrared Imaging Radiometer Suite (VIIRS) instrument, equipped with a calibrated Day/Night Band (DNB), opens a new door for enabling shortwave nighttime aerosol retrievals from spaceborne observations (e.g., Johnson et al., 2013).
The previously mentioned studies are motivated by an acutely felt need for nighttime aerosol observations.In particular, studies have shown notable improvements in aerosol forecasting through the assimilation of satellite aerosol products, mostly from daytime observations (e.g., Zhang et al., 2008aZhang et al., , 2011Zhang et al., , 2014;;Yumimoto et al., 2008;Uno et al., 2008;Benedetti et al., 2009;Schutgens et al., 2010;Sekiyama et al., 2010).To capture the diurnal cycle, the aerosol modeling community requires nighttime satellite aerosol data hav-T.M. McHardy et al.: An improved method for retrieving nighttime aerosol optical thickness ing broad spatial coverage and high temporal resolution to further advance aerosol, visibility, and air quality forecasts (e.g., Zhang et al., 2011Zhang et al., , 2014)).Johnson et al. (2013) were the first to outline a method for retrieving nighttime aerosol optical thickness (AOT, τ ) from the VIIRS DNB.In Johnson et al. (2013), the difference between VIIRS DNB radiances from an artificial light source and a dark background area adjacent to the light source over a cloud-free sky are used to derive nighttime τ .
The algorithm used in the nighttime aerosol retrieval method of Johnson et al. (2013) is challenging to apply on a global scale for several reasons.First, nearly identical VIIRS pixels from different nights are needed from a selected artificial light source and the corresponding nearby background.To implement such a method operationally over a large study domain is difficult, as VIIRS DNB viewing geometries vary from day to day.This makes it difficult to select the same group of pixels daily.Second, the inherent variation within an artificial light source, especially for a large city comprised of hundreds of VIIRS pixels, is not considered.Thus, larger retrieval errors may be expected from complex artificial light sources.Third, the contribution from diffused artificial lights needs to be computed, requiring a priori knowledge about aerosol types and aerosol-absorbing properties.
In this study, modification to the Johnson et al. (2013) algorithm is proposed which considers the reduction of radiance contrast of VIIRS DNB pixels within a given nighttime artificial light source as a proxy for aerosol optical depth.This new method overcomes some of the issues encountered by Johnson et al. (2013) and is easier to implement on a global scale.The outline of this paper is as follows: Sects.2, 3, and 4 describe the data used in this study, highlight the newly developed approach, and explore the sensitivity of radiance dispersion within a given light source to various parameters, respectively.Section 5 shows the intercomparisons of the retrieved nighttime τ from the VIIRS DNB with AERONET and High Spectral Resolution Lidar (HSRL) observations.Considerations of uncertainty are included in Sect.6, and Sect.7 summarizes the salient findings of this research.

Data sets
This study combines VIIRS DNB satellite observations with ground-based sun photometers from four sites in the AERONET network and an HSRL instrument in Huntsville, AL, deployed for the NASA Studies of Emissions and Atmospheric Composition, Clouds and Climate Coupling by Regional Surveys (SEAC 4 RS) campaign.Table 1 lists the locations, study periods, and instrumentation at Alta Floresta, Cape Verde (now officially Cabo Verde), and Grand Forks (the three sites used in Johnson et al., 2013) as well as a fourth site in Huntsville, AL, which is new to this study.Level 1.5 AERONET data were chosen in Johnson et al. (2013) as Level 2.0 data were not available for all desired sites during their respective study periods.We choose to use Level 2.0 AERONET data as they are now available and have included updated comparisons using data from Johnson et al. (2013) to maintain consistency between the two studies.
The nighttime AERONET τ values are estimated by taking the mean of the daytime AERONET retrievals straddling the VIIRS nighttime overpass (the daytime AERONET τ values from the evening before the time of the VIIRS overpass and the morning after the time of the VIIRS overpass, as long as both AERONET retrievals are within 24 h of each other).If, for any reason, there is no available nighttime AERONET τ estimation, the corresponding VIIRS overpass is not included in this study.
Cloud-Aerosol LIdar with Orthogonal Polarization (CALIOP) nighttime τ data are not used as a validation tool in this study due to a lack of collocated CALIOP and VIIRS data pairs over the selected artificial light sources.Also, Omar et al. (2013) found significant τ biases and differences in cloud clearing between CALIOP and AERONET.CALIOP will have a better signal-to-noise ratio for nighttime observations vs. daytime, but Campbell et al. (2012) found that CALIOP τ had similar errors during both day and night.Even with a perfect CALIOP τ , the small spatial coverage combined with the long repeat cycle eliminates CALIOP from the list of potential validation tools for this study.
Still, validating VIIRS nighttime τ retrievals with daytime AERONET data is not ideal, as the diurnal variation of τ is not considered.Thus, nighttime Wisconsin HSRL τ retrievals from Huntsville, Alabama (from 26 June to 26 October 2013), are used as an independent validation data set.Different from CALIOP, HSRL measures both molecular and aerosol backscattering, and thus can be used to derive the ratio of extinction to backscatter (lidar ratio), which can be further used to yield a more reliable τ retrieval (Hair et al., 2008).Besides HSRL measurements, AERONET data are also available for the study period within close proximity to Huntsville, Alabama (UAHuntsville station, 34 • N, 86 • W).Thus, both Level 2.0 AERONET and HSRL data from Huntsville, Alabama (from June-October 2013), are used in the validation efforts.
The HSRL data are processed as a 10 min average at 30 m vertical resolution (i.e., each τ value is the 10 min average of extinction integrated from the surface up to a given altitude, every 30 m from the surface).Because of noise at higher altitudes, rather than simply taking the column-integrated τ value at the highest altitude, HSRL total column τ was estimated by taking the average of all column-integrated τ retrievals between altitudes of 10 and 15 km (Ralph Kuehn, personal communication, 2015).HSRL profiles are visually inspected for the presence of clouds, and VIIRS overpasses that occurred on nights with no available HSRL data are not used.
In order to obtain τ 0.532 , Eq. ( 1) is rearranged such that where, in this case, λ is 0.532 µm; λ 0 is 0.675 and 0.7 µm (the center of the DNB wavelength range) for AERONET and VIIRS retrievals, respectively; and τ λ 0 is the originally retrieved τ at those wavelengths.
The three sites (Alta Floresta, Cape Verde, and Grand Forks) from Johnson et al. (2013) were originally chosen to represent smoke, dust and, urban cases, respectively.It is important to note that they also vary in spatial extent (6-150 pixels).The Huntsville test site covers a significantly larger area (500 pixels) than the other three sites, allowing for detailed testing of the impact of city size upon τ retrievals.Basic information about the four sites examined in this study is found in Table 1.
The VIIRS/DNB Sensor Data Record-SDR (SVDNB) and VIIRS/DNB SDR Geolocation Content Summary (GDNBO) data are used for providing radiance and geolocation values, respectively.The VIIRS Cloud Cover/Layers Height Data Content Summary (VCCLO) data are used for cloud clearing of the selected VIIRS pixels.In addition to the calibrated pixel-level top-of-atmosphere radiance values, the VI-IRS DNB data products include metadata on satellite and lunar geometry, specifically, satellite, solar, and lunar zenith angles; lunar illumination fraction; and quality assurance flags.

Methodology
Following Johnson et al. (2013), the satellite-observed radiance is composed of three major components, as illustrated in Eq. (3).
Here, τ is the total optical thickness, µ is the cosine of the satellite zenith angle, D is the additional diffuse radiance, I p is the path radiance (including aerosol-layer-reflected moonlight), and I s is the surface upward radiance for an artificial light source.The first term from the right-hand side of Eq. ( 3) represents the VIIRS-observed surface radiance through direct attenuation.The second and the third terms from the right-hand side of Eq. ( 3) show the diffuse and path radiances, respectively (Johnson et al., 2013).As mentioned in Johnson et al. (2013), I s can be defined as where F 0 is the moonlight at the top of the atmosphere, µ 0 is the cosine of the lunar zenith angle, r s is the surface reflectance and r is the reflectance from the aerosol layer.T (µ 0 ) is the diffuse transmittance, and I a is the emission from the artificial light source.The first two terms from the numerator represent reflected direct and diffuse downward moonlight, respectively.Here we assume that emissions from artificial light sources (I a ) are stable throughout the respective study periods.Also, because pixels comprising a given city are within close proximity, we assume that total optical thickness, diffuse radiance, path radiance, and the reflected direct (r s µ 0 F 0 e τ/µ 0 ) and diffuse (r s µ 0 F 0 T (µ 0 )) downward moonlight are invariant within a given city.Like in Johnson et al. ( 2013), the r s r term is assumed to be negligible, although this term can be significant for optically thick aerosol plumes.We further assume that the retrieved total optical thicknesses over cloud-free skies are aerosol optical thicknesses, although a detailed analysis is needed to separate AOT from Rayleigh optical thickness (Johnson et al., 2013).
Equation 6 shows that optical thickness can be calculated for a single artificial light source using only the cosine of the satellite zenith angle (µ) and the pixel-to-pixel differences in satellite-observed radiance and upwelling radiance from the artificial light source ( I sat and I a , respectively).C is a correction term applied to account for a dependence of upwelling radiance on satellite viewing angle, which was found to be applicable only at the Grand Forks location as identified from Johnson et al. (2013), for plausible reasons such as stray light contamination.In Johnson et al. (2013), a linear relationship between cosine viewing angle and VIIRS radiance is constructed, and for each observation the C value is estimated by dividing the observed radiance by the radiance value predicted from the linear relationship for the Grand Forks site.We have adopted C values directly from Johnson et al. (2013).C is set to 1 for all other sites.A more detailed explanation of this correction term can be found in Johnson et al. (2013).It was suspected that this dependence on viewing geometry was a product of parallax effect, which was corrected with the release of terrain-corrected geolocation data (released in early 2014).The difference between terrain-corrected and non-terrain-corrected geolocation data was determined to be subpixel in magnitude for Grand Forks.Therefore, the lack of terrain-corrected geolocation for the study periods in this study is not expected to introduce additional error.
As this study intends to examine the radiance contrast between pixels within a given artificial light source, a way to quantify the pixel-to-pixel differences of radiance values within this light source is needed.Although, mathematically, Eq. ( 5) is achieved by taking the spatial derivative of Eq. (3), it is beneficial for this method to look at the spatial derivative in the statistical domain.In other words, a measurement of statistical dispersion is used in place of taking the actual derivative.A sensitivity study comparing optical thickness values calculated using different measures of statistical dispersion (not shown) determined that simple standard deviation produces the best results when compared against AERONET optical thickness values; therefore the operator in Eqs. ( 5) and ( 6) is the "the standard deviation of . . .".The sensitivity of τ in Eq. ( 6) to varying I sat and satellite zenith angle (with a fixed I a ) is shown in Fig. 1.
The method described above, dubbed the "variance method", is based on radiance values within a given artificial light source; therefore, it is necessary to determine which pixels comprise an artificial light source ("city pixels").With potential automation in mind, it was decided that city pixels should be determined algorithmically.As suggested from Johnson et al. (2013), any pixel with a radiance value greater than 1.5 times the mean radiance of a given section of a VI-  6) with varying standard deviations of radiance ( I sat , x axis) and satellite zenith angles, simulating retrievals made using the variation method.I a is held constant at 10 −7 W cm −2 sr −1 .IRS DNB granule enclosing the target city and with a radiance value greater than 0.25 ×10 −8 W cm −2 sr −1 is considered a city pixel.Figure 2 shows an example of artificial light detected using this method for Huntsville, Alabama.Because atmospheric conditions, and therefore VIIRS radiance values, change on a nightly basis, the number of city pixels for a given city may also change from night to night.In order to remain consistent, n brightest city pixels are used to compute the standard deviation of radiance values for each nighttime VIIRS overpass within the study period of a given city, where n is the minimum number of algorithmically determined city pixels among overpasses that passed VIIRS cloud screening and visual inspection for each city.VIIRS cloud screening is performed by removing any pixels that were flagged as containing clouds at any level in the previously mentioned VCCLO product.Any VIIRS granule of an artificial light source that contains a cloudy pixel in the vicinity of the artificial light source (∼ 0.2 • latitude/longitude) is not used in this analysis.
The only values required to calculate τ from Eq. ( 6) over an artificial light source from a single VIIRS nighttime overpass are I sat and I a -the standard deviations of radi-ance values of city pixels from that overpass and the stable (aerosol-, cloud-, and moon-free) standard deviation of radiance values of the city pixels, respectively.The former is straightforward; the latter, however, presents multiple challenges.First and foremost, a nighttime atmosphere totally free of aerosols and lunar illumination is very unlikely, so error is inherent.Due to the lengths of the VIIRS data record and study periods used for this demonstration, the number of VIIRS overpasses that occur on even nearly aerosol-, cloud-, and moon-free nights within the study periods is extremely limited.It was therefore decided to take I a for each city as the mean of the two greatest standard deviations of city pixel radiance values within the study period.The rationale behind this is that, in theory, the night with the lowest aerosol loading and least lunar illumination should have the most variance among city pixels.This is not always the case; however, the I a terms for all locations are taken on nights with τ values below 0.1 (according to the straddling daytime-averaged AERONET τ ).A database of VIIRS radiance values for artificial light sources worldwide is currently being constructed to improve the robustness of this parameter.
Lastly, different from Johnson et al. (2013), only pixels within an artificial light source are needed.It is assumed that the diffuse radiance is spatially and temporally invariant within a given light source, and the derivative is 0. Thus, the newly proposed approach avoids the need for a priori knowledge of aerosol microphysical properties, making it more suitable for a larger-scale analysis.

Sensitivity studies
The contrast of satellite-observed radiances within an artificial light source ( I sat ) is used as the basic information content for the study.It is therefore necessary to examine the impacts of viewing angles and lunar illumination on I sat .For example, for a given aerosol loading, do we expect a change in I sat for a moonless night vs. a night with significant moonlight?Figure 3 shows the variance of VIIRS radiance values of artificial light source pixels plotted against lunar fraction (from the VIIRS GDNBO product).There is a plausible relationship between variance and lunar fraction for the Alta Floresta, Grand Forks, and Cape Verde sites, but the Huntsville site displays no such relationship.Lunar zenith angle (from the VIIRS GDNBO product) was examined in the same manner (not shown), but the data showed no plausible relationship with variance for any location.A variance correction using a linear model was attempted for the three locations showing a potential relationship; however this correction produced results that are worse than those without the correction, and thus it is not used hereafter.
It is anticipated that lunar illumination would introduce a change in I sat .To test this assumption, the I sat values over the Huntsville site are computed for nights with and without lunar illumination that have similar nighttime HSRL τ val- ues.Thus, if the lunar illumination affects I sat under a similar aerosol loading scenario, a relationship between I sat and lunar fraction should be observed.An examination of the difference in the variance of VIIRS radiance values between nights with lunar illumination and nights without lunar illumination is shown in Fig. 4.Each data point consists of a pair (one light, one dark) of retrievals with similar HSRL total-column τ (within 0.008 of each other).Also, to minimize the effects of aerosol loading on the derived relationship, only data pairs that have HSRL τ values less than 0.25 (and greater than 0.16, the lowest HSRL τ among retrieval nights) are chosen.The y axis displays the relative difference in variance for each pair -this is a proxy for the contribution of lunar illumination to the VIIRS τ retrieval.The x axis is the difference in lunar fraction for each pair.Values of relative contribution range from approximately −0.12 (meaning the light night had lower variance) to approximately 0.42.While the sample size is small, this figure displays no clear relationship between lunar fraction and VIIRS radiance dispersion.
To further test the effects of lunar illumination on the VI-IRS retrievals, all retrievals were separated into two groups -those made during which the moon was present and those made without the presence of the moon.Retrievals made in the presence of the moon were determined as all nights having a lunar fraction greater than 15 % and having a lunar zenith angle below 89 • .Table 2 shows the root-mean-squared error (RMSE) of the VIIRS-retrieved τ values against straddling daytime-averaged AERONET τ values separated by the presence of the moon, as well as data counts for each subset separated by location.Retrievals were not divided further (such as in Fig. 6) in order to keep the sample size in each subsection sufficiently large and to prevent locational bias.The RMSE for retrievals made on nights with significant lunar illumination is nearly 0.1 lower than for nights without significant lunar illumination.With this small sample size, lunar illumination cannot be ruled out as a source of retrieval error; however, our tests find no evidence for degraded retrieval performance under moonlit conditions.
Similar to Johnson et al. (2013), the relationship in between I sat and viewing geometry is studied.Figure 5 shows I sat as a function of the cosine of the satellite zenith angle for Huntsville, Alta Floresta, Grand Forks, and Cape Verde.Similar to what is reported by Johnson et al. (2013), no viewing angle dependence is found for I sat for Alta Floresta and Cape Verde.Huntsville also does not display a dependence on viewing geometry.Again, as in the previously mentioned study, a significant positive relationship does exist between I sat and the satellite zenith angle for Grand Forks.The reason for this relationship, as suggested from Johnson et al. (2013), is still not known and is left to future research once a larger analysis sample is available.
A major caveat regarding the previously discussed sensitivity studies that has been mentioned is sample size.While these sensitivity studies appear to be relatively inconclusive, the study period has not been extended to achieve statistical robustness for a few reasons.The first is that the primary goals of this study are to demonstrate the efficacy of the variance method and to compare the results directly with the results presented in Johnson et al. (2013).Second, the short study periods are selected to ensure the relative stability of artificial light sources.For a longer study period, the seasonal variations in artificial light sources will need to be accounted for, which is beyond the scope of this paper and is a subject of our next planned study.6a, circle), Grand Forks (left, Fig. 6a, "x"), and Cape Verde (right,Fig. 6b,asterisk).One-to-one (dotted) and bestfit (solid) lines are also shown.Cape Verde is isolated due to its relative small size (6 pixels).Retrievals made on nights with a lunar fraction less than 15 % or with a lunar zenith angle greater than 89 • (below the horizon) are shown in black.Retrievals made on nights with a lunar fraction greater than 15 % but less than 50 % and a lunar zenith angle less than 89 • are shown in red.Retrievals made on nights with a lunar fraction greater than 50 % and a lunar zenith angle less than 89 • are shown in blue.

Results
As the first step, the VIIRS DNB nocturnally retrieved τ values at 0.7 µm (τ 0.7 ) derived based on Eq. ( 6) are intercompared with the straddling daytime-averaged AERONET τ at 0.675 µm (τ 0.675 ) values. Figure 6a shows VIIRSretrieved τ 0.7 values plotted against the corresponding straddling daytime-averaged AERONET τ 0.675 for Huntsville (square), Alta Floresta (circle), and Grand Forks ("x").Cape Verde (Fig. 6b, asterisk) is isolated, as it is comprised of significantly fewer pixels than the other locations.It is suspected that the cause of the poor performance of the VIIRS I a (the "stable" upwelling radiance from the artificial light source) for each city is taken from the average of two nights within the study period as explained above, there is one VI-IRS τ retrieval for each location that is at or below 0, which is not physically possible, and one VIIRS τ retrieval slightly above 0.These retrievals were not removed for the calculations of r 2 , RMSE, or best-fit lines, and they do negatively impact the results.
One of the main purposes of this study is a comparison between the variance method developed here and the method presented in Johnson et al. (2013, hereafter "background method").Table 3 shows the r 2 , RMSE, and bestfit slope of nighttime τ 0.7 retrievals made using the variance method (Table 3a) and nighttime τ 0.7 retrievals made using the background method (Table 3b, Johnson et al., 2013) compared against the corresponding straddling daytime-averaged AERONET τ 0.675 values.Huntsville data are not included in Table 3b, as they are not reported in Johnson et al. (2013).For Alta Floresta (moderate aerosol loading, moderate size), the variance method performs better than the background method across the board.For Grand Forks (low aerosol loading, large size), the methods perform comparably.We suspect the poor performance of both methods over Grand Forks is due to the very low aerosol loading and relatively high error inherent in these methods.For Cape Verde (high aerosol loading, very small size), the background method performs significantly better than the variance method.Based on these results, it is hypothesized that the background method performs better for smaller or more isolated artificial light sources and that the variance method performs better for larger artificial light sources.This hypothesis will be further investigated once a large-scale analysis has been performed.
We have further evaluated the variance method with the use of HSRL data over Huntsville. Figure 7a shows VIIRSretrieved τ 0.7 values interpolated to 0.532 µm (τ 0.532 ), plotted against the corresponding time-averaged HSRL-retrieved total-column τ 0.532 values for Huntsville.The r 2 of the interpolated VIIRS-retrieved nighttime τ 0.532 values and averaged HSRL-retrieved total-column τ 0.532 values is 0.48, with a RMSE of 0.08 (with a mean HSRL τ of 0.25, for context), indicating that the VIIRS-based nighttime τ retrieval method has some skill.We have also plotted the straddling daytimeaveraged AERONET τ 0.532 (interpolated to 0.532 µm) vs. the averaged HSRL τ 0.532 values as shown in Fig. 7b.Clearly, Fig. 7b suggests there is value to using the average of the "bookend" daytime AERONET τ observations as nighttime ground truth when there is a lack of better options.Figures 6 and 7 also show that, while far from perfect, the method presented in this paper has some skill at capturing variations in nighttime τ .
To further investigate the effects of lunar illumination on the VIIRS τ retrievals, relative errors (VIIRS retrieval, interpolated, minus HSRL total-column, divided by HSRL totalcolumn τ ) for the Huntsville retrievals are plotted as a function of HSRL total-column τ as shown in Fig. 8.The lunar fraction (scaled from 0-100) on the night of the retrieval is indicated by the plot symbol size.This figure shows that there is no apparent relationship between retrieval error and lunar fraction.It is clear from Fig. 8, however, that the relative error of the VIIRS retrievals decreases as HSRL total-column τ increases, as expected.This indicates that VIIRS retrievals  7a) and VIIRS-retrieved τ (bottom, Fig. 7b) as a function of Wisconsin HSRL-retrieved total-column τ for Huntsville, Alabama.One-toone (dotted) and best-fit (solid) lines are also shown.are more valuable for aerosol loading events with τ greater than approximately 0.25.
Because this method performs the worst for Cape Verde (the smallest site by number of city pixels), a sensitivity study is performed for Huntsville, Alta Floresta, and Grand Forks to determine if the number of city pixels used in a retrieval is an important factor in its quality.Cape Verde is not used because reducing the already small number of pixels produces results that are not robust.To achieve this goal, the relationship between the number of artificial city light pixel using in the VIIRS τ retrievals and the I sat values are studied over relatively aerosol-free nights (HSRL and/or AERONET τ less than 0.2; see Table 4).For each location, the n brightest pixels of the original algorithmically determined city pixels (i.e., the data that are used to create Figs. 6 and 7, also included in Table 4) are treated as the baseline of the study.For a given location, I sat values are computed for the identified relatively aerosol-free nights using the baseline-defined number of brightest artificial light pixels.Then, a standard deviation of I sat values ( I sat ) is computed.I sat represents the relative stability of I sat values used in the VI-IRS τ retrievals.Next, we repeat the process by incrementally reducing the number of artificial light pixels used in the analysis.For example, we reduce the number of pixels used for a Huntsville retrieval to 400 (as opposed to the original 509 pixels used) by taking the 400 brightest pixels out of the original 509 pixels.I sat is then computed for each night in Table 4, with the 400 brightest pixels treated as the artificial light source.The result of the exercise is shown in Fig. 9.In Fig. 9, the x axis is the number of brightest artificial light pixels used in computing I sat normalized by the number of pixels used in the actual retrieval for each respective location (see Table 4).The y axis is I sat normalized by the I sat from the baseline cases.It is clear from Fig. 9 that the normalized I sat increases as the number of pixels used in the retrieval decreases.This increased variability translates into a less reliable retrieval.It is hypothesized that this is because simple standard deviation as a quantification of statistical dispersion is more robust as the sample size increases.Therefore, it is suggested that only medium to large cities (on the order of tens of VIIRS pixels in spatial coverage at least) are used for nighttime retrievals via the method outlined in this study; however, an exact size has yet to be determined and is a topic left for future studies.

Uncertainties and limitations
Similar to Johnson et al. (2013), by taking the total derivative of Eq. ( 6  ) as a function of the normalized number of brightest artificial city light pixels (n) for nights that are relatively aerosol-free (HSRL and/or AERONET τ < 0.2) for Huntsville, Alta Floresta, and Grand Forks.For a given location and for a given number of brightest artificial city light pixels taken from among the baseline case pixels, the standard deviations of I sat values are first computed using data from selected relatively aerosol-free nights.Both n and I sat value are then normalized based on the values used in the actual retrievals (n , I sat ) as shown in Figs. 6 and 7 (details are listed in Table 4) for a comparison among different locations.
Thus, the uncertainty in the estimated VIIRS τ is a combined relative uncertainty from the correcting factor C (only applied to Grand Forks), as well as relative errors in the standard deviation of satellite radiance and surface outgoing radiance within an artificial light source, weighted by the cosine the viewing angle.For example, a 10 % uncertainty in the derived I s could introduce an error in VIIRS τ of 0.05 at the viewing angle of 60 • (µ = 0.5), and at nadir the retrieved error increases to 0.1.
Still, there are major issues that need to be explored that could potentially limit implementation of the research presented in this paper.For example, I a , which is the surface artificial light source emission, is assumed to be invariant through the study period.However, I a may indeed, and almost certainly does, change with time and viewing geometry (e.g., Román and Stokes, 2015).Even with a constant satellite zenith angle, viewing direction could present an addi-tional issue, especially when viewing artificial light sources.Thus, to fully investigate this problem, a careful analysis of the spatial and temporal variations of artificial light sources is needed.Also, please note that the uncertainties at the low and high τ ranges could be dominated by different factors.In a high-τ regime, the ignored r s r term can be significant.In a low-τ regime, however, the r s r term can be ignored, while the ignored Rayleigh optical depth can be comparable in magnitude to the retrieved AOT.Other currently unresolved factors related to the artificial light sources that may introduce uncertainty include an unexpected increase in the artificial light emissions (e.g., a large-scale fire event).In this study, city pixels are determined algorithmically using purely radiance values, which means that the exact same locations on the earth were not necessarily chosen for each retrieval.While there may be a way to implement official boundaries into determining which pixels are city pixels, cities often do not simply end at boundaries.These sources of uncertainty are in addition to previously mentioned relationships due to lunar and satellite variables.
Cloud contamination is likely to be another major source of uncertainty.With only limited visible and infrared channels, nighttime clouds, and especially thin cirrus clouds, detections may be problematic.Since this is a case study, all VIIRS granules are visually inspected.However, to fully automate the process, additional cloud screening checks will most likely be needed, which will require the use of other VIIRS channels or data from other sensors.

Conclusions
This paper presents a new method for using radiance values observed by the VIIRS DNB over artificial light sources to estimate nighttime aerosol optical thickness.This method is based on theoretical radiative transfer equations for the retrieval of optical thickness using contrast reduction of radiances from artificial light sources within a close proximity.This study seeks to improve upon the shortcomings inherent in the method presented by Johnson et al. (2013), such as the difficulty of implementing said method over a large area in a timely manner.This study suggests the following: 1. Improved agreements vis-à-vis Johnson et al. (2013) are found between VIIRS DNB nighttime τ and nighttime τ values estimated from AERONET and HSRL measurements.This indicates that the variance method can be used as a way of estimating nighttime aerosol properties.
2. A reasonable agreement is also found between nighttime HSRL τ values and nighttime AERONET τ values that are approximated based on daytime AERONET data.Thus, daytime AERONET τ observations can be used semi-quantitatively as validation for nighttime τ retrievals when other options are not present.
3. The proposed algorithm has its limitations.Expected uncertainties originate from lunar illumination, viewing geometry, and city characteristics.Various tests described in this study did not show any systematic impact of lunar illumination on the retrieval, but the small sample size prevents drawing strong conclusions from this.Also, to implement the algorithm in a regional or global analysis, an automated cloud screening scheme is needed to take the place of the manual cloud screening performed in this study.
4. While it is expected that lunar illumination has an effect on the retrievals presented in this study, the magnitude of contribution to retrieved τ values from lunar contamination is unclear.For example, data from the largest artificial light source (Huntsville), which has coincident nighttime HSRL measurements, display no relationships between lunar fraction and VIIRS-retrieved τ values, relative error, or variance of radiance values.
Lastly, an operational nighttime aerosol product from a passive remote-sensing method is still non-existent.CALIOP does provide nighttime aerosol measurements, but the coverage of the non-scanning CALIOP observations is very limited (a single track roughly 70 m across).This limits the use of CALIOP data in aerosol modeling and forecasts.This study presents a novel technique for deriving nighttime aerosol properties from the VIIRS that, despite important limitations, shows skill in comparison with ground-based measurements.The methods described here will lead to future improvements in developing a reliable nighttime aerosol product for the aerosol modeling community.
5) I sat and I a are the pixel-to-pixel change in satelliteobserved radiance and surface emission for a single artificial www.atmos-meas-tech.net/8etal.: An improved method for retrieving nighttime aerosol optical thickness light source, respectively.Solving Eq. (

Figure 1 .
Figure1.AOT computed by using Eq.(6) with varying standard deviations of radiance ( I sat , x axis) and satellite zenith angles, simulating retrievals made using the variation method.I a is held constant at 10 −7 W cm −2 sr −1 .

Figure 2 .
Figure 2. (a) A Raw VIIRS DNB image of Huntsville from the night of 26 June 2013.(b) The same image from (a) with algorithmically determined city pixels in red.(c) The same as in (b) except that only the 100 brightest city pixels are in red.

Figure 3 .
Figure 3. Variance of VIIRS radiance values vs. average lunar fraction for each night an τ retrieval was made, separated by location.

Figure 4 .
Figure 4.The relative difference in radiance variance for nights with lunar contamination ( I on the y axis label) and without lunar contamination ( I ).Each point consists of a pair (one light, one dark) of retrievals with similar (within 0.008) τ .The x axis indicates the difference in lunar fraction between the paired retrievals.The symbol sizes represent the magnitude of HSRL τ values, which range from 0.15 to 0.25.

Figure 5 .Figure 6 .
Figure 5. Variance of VIIRS radiance values vs. the cosine of average satellite zenith angle for each night that a retrieval is made, separated by location.

Figure 8 .
Figure 8. error of each VIIRS τ retrieval (interpolated to 0.532 µm) with respect to HSRL-retrieved total-column τ , plotted as a function of the HSRL retrieval.The lunar fraction (on a scale of 0-100) is indicated by plot symbol size, as shown by the key in the bottom right corner.
), we have dτ = µ dC C + d I a I a − d I sat I sat .

Table 1 .
Basic information for each of the test sites, including site name, AERONET and HSRL site latitude-longitude locations, and date range of the respective study periods.The "Size of the city" column indicates the number of pixels used in the calculations.

Table 2 .
RMSE of VIIRS-retrieved nighttime τ against straddling daytime-averaged AERONET τ for retrievals occurring on nights with (top) and without (bottom) the moon present.Nights with lunar presence are determined as those having a lunar fraction greater than 15 % and having a lunar zenith angle less than 89 • .The number of retrievals in each category separated by location is also shown.

Table 3 .
Coefficient of determination, root-mean-squared error (RMSE), and best-fit slope of VIIRS estimated nighttime τ compared against estimated nighttime Level 2.0 AERONET τ for the new method presented in this study (a) and the method presented inJohnson et al. (2013) (b).Average AERONET τ for each location is provided for context.

Table 4 .
Dates with HSRL (for Huntsville) or straddling daytime-averaged AERONET (Grand Forks and Alta Floresta) τ less than 0.2, which are used for the computations shown in Fig.9, and the baseline (i.e., used in the actual τ retrieval) n pixels used for Grand Forks, Huntsville, and Alta Floresta.Date format: month/day.