Aerosol characterization using satellite remote sensing of light pollution sources at night

A demanding challenge in atmospheric research is the night-time characterization of aerosols using passive techniques, that is, by extracting information from scattered light that has not been emitted by the observer. Satellite observations of artificial night-time lights have been used to retrieve some basic integral parameters, like the aerosol optical depth. However, a thorough analysis of the scattering processes allows one to obtain substantially more detailed information on aerosol properties. In this Letter we demonstrate a practicable approach for determining the aerosol particle size number distribution function in the air column, based on the measurement of the angular radiance distribution of the scattered light emitted by night-time lights of cities and towns, recorded from low Earth orbit. The method is self-calibrating and does not require the knowledge of the absolute city emissions. The input radiance data are readily available from several spaceborne platforms, like the VIIRS-DNB radiometer onboard the Suomi-NPP satellite.


INTRODUCTION
Determining the aerosol properties at night-time is an essential step for a better understanding of aerosol dynamics, with direct applications to the study of planetary atmospheres, the characterization of potential candidate sites for astrophysical observatories, and light pollution research. Artificial night lights of cities and towns offer a permanent set of light beacons distributed worldwide, whose observation from Earth orbiting platforms has been shown to be useful to determine some basic integral properties of the aerosols contained in the air colum, as e.g. the aerosol optical depth (Choo and Jeong 2016;Johnson et al. 2013;McHardy et al. 2015;Wang et al. 2016;Zhang et al 2008Zhang et al , 2019. However, night-time imagery contains additional useful information, among other the angular dependence of the signal produced by the artificial light scattered by the atmosphere, that can be measured in clear and moonless nights by observing the radiance received at the satellite from directions corresponding to Earth pixels with no light sources, preferably located in the vicinity of urban nuclei with sharp geographical borders. This scattered light appears in the satellite images as a diffuse glow surrounding the city areas, blurring the city ⋆ E-mail: kocifaj@savba.sk limits, and with decreasing radiance for increasing distances (e.g. Sánchez de Miguel et al. 2019). In this Letter we show that the specific way in which the scattered radiance varies with the nadir angle of the observed pixels, as seen from the satellite radiometers, can be used to retrieve more specific key aerosol properties like e.g. their particle size number distribution function, if some additional a-priori information is available. An important feature of this approach is that it is self-calibrating: since it is based on the analysis of the angular dependence of the scattered radiances normalized to the direct radiance received from the city, the knowledge of the absolute emissions from the urban area is not required for its application.

THE MODEL
A light source seen in satellite imagery of an atmosphereless planet would be geometrically sharply confined. However, due to light scattering by atmospheric constituents, the brightness of the artificially lit surface of the Earth at night does not change abruptly when transitioning from bright image pixels to neighbouring pixels at the outer interface of a city or town. Normally the radiance of a surface decreases steeply but in a continuous way, as the angular distance c 2020 The Authors from the city-edge increases, while quickly approaching the background level.
The radiance of the Earth's surface, RS(z), reaching the satellite sensor can be computed as the integral of elementary optical signals produced from light scattering in the atmospheric volumes along the line of sight (see Fig. 1a) where z ≈ z0 + Θ is the zenith angle of a surface area in the city surroundings that has no light sources, while z0 is the zenith angle of an illuminated surface element within the city or town. In Fig. 1a Θ is the angle between the direction of observation and the position of a point s ource of light (or a surface element dS of an artificially lit area). The normalized scattering phase function p(θ)/(4π) shapes with the scattering angle θ cos θ = cos z cos z ′ + sin z sin z ′ cos(∆ϕ) .
The angles z, z ′ and ∆ϕ are defined in Fig. 1a. In Eq. 1 the product of the single scattering albedo,ω, and the volume extinction coefficient, kext, is known as the atmospheric volume scattering coefficient (see e.g. Kokhanovsky 1998). The optical thickness τ (h) is calculated for an atmospheric column that extends from altitude h to the satellite level ξ, i.e.
is the radiant intensity (measured in W sr −1 ) emitted by an elementary area dS of the city, that can be expressed as whereR0(z ′ ) has the units of radiance (W m −2 sr −1 ). The radiant intensity emitted by a city normally depends on the emission angle z ′ . The collective effect of differently oriented and sloped surfaces, however, is that the intensity usually varies slowly with z ′ . This is consistent with the models currently in use (e.g. Fig. 5 in Kolláth and Kránicz 2014). The angular dependences of the remaining terms in the integrand in Eq. 1 are much more pronounced, e.g. p(θ) shows a dramatic decrease rate by one or more orders of magnitude across its definition domain, with a typical strong forward lobe and a weak side scattering. After a considerable manipulation of mathematical expressions (not shown here) we obtain for the angular dependence of the scattered radiance detected by the satellite at wavelength λ. In this equation it is assumed that the satellite passes almost vertically over the city (i.e. at near zero zenith angle), and that the radiant intensity emitted by the city is slowly dependent on z ′ , having an average valueR0dS. Due to atmospheric extinction the radiance of the bright pixels is reduced tõ R 0ξ = e −τ (0)R 0 when measured at the satellite level. The approximations used in deriving the above formula also comprise i) θ ≈ z ′ + Θ ≈ z ′ , ii) z = Θ, iii) the contribution of Rayleigh scattering to the radiance recorded at satellite level for Θ ∼ < 0.09 • is much smaller than the one from aerosols (see Fig. 2), and iv) the Mie theory for spherical homogeneous particles of radii a controls the far-field scattering amplitude through the first element of the scattering matrix, S11, which is also known as the scattering function. The scattered signal from an ensemble of independent particles is proportional to the product of S11(a, z ′ ) and f (a, h), where the latter is the particle size number distribution function, i.e. the number of particles per unit volume having radii between a and a + da at height h above ground level. More accurately, when vegetation and city structures or outdoor luminaire designs are all not efficient enough in blocking light emissions to low elevation angles, an exponential term should be kept in Eq. 4, i.e.

RS(Θ)
The radiance RS(Θ) increases steeply when Θ approaches zero, and we evidence in Fig. 2 that the aerosol contribution is dominant (in comparison with Rayleigh's molecular one) in shaping the radiance patterns for almost all emission angles z ′ and for most typical values of the aerosol asymmetry parameter g. The latter is commonly used in radiative transfer theories to characterize the average cosine of the scattering angle, assuming that the probability of a photon to be scattered per unit solid angle around the direction z ′ is proportional to the phase function p(z ′ ). In the numerical demonstration we have varied the asymmetry parameter in order to model different aerosol types. Values of g approaching unity usually indicate the presence of large particles with a strong forward scattering lobe, while g ≈ 0 corresponds to isotropically scattering media containing particles much smaller than the spectral detection range (e.g. Moosmüller and Ogren 2017). Most typically, g ranges from a few tenths to ≈ 0.9. This is why the numerical demonstration of aerosol contribution to the satellite radiance is performed for g=0.3 ... 0.9.

THE INVERSE SCATTERING PROBLEM
For a satellite located at near zero zenith angle above the city and θ ≈ z ′ we have dh sin 2 z ′ = −ξ sin Θdz ′ , while the mapping from the space of experimental radiance data, Rs(Θ), to the space of particle size distribution functions f0(a) is described by the following integral equation where Θ1 is the smallest angular separation between the bright surface area of a city and the direction of remote sensing, and RS(Θ1) is the radiance scattered towards the satellite from the whole atmospheric column along the line of sight ending in point M (see Fig. 1b). For a vertically stratified atmosphere f (a, h) = f0(a)F (h), where F (h) is commonly approximated by an exponential function F (h) = exp {−h/HA} with HA being the aerosol scale height (e.g. Waquet et al. 2007). In most cases HA ranges from 1 to 3 km. The kernel KS11(Θ1, a) of the integral equation has the following form for the satellite remote sensing of a light pollution source where ∆Θ = ΘN − Θ1 is the angular span between the nearest and farthest bright areas in the city, measured from the point M in the satellite images (cfr. Fig. 1b). ǫ(Θ) is the halfangle subtended by each elementary arc of the city emitting area, as seen from the ground point M , for each value of Θ (see arcs within dS in Fig. 1b). The inverse problem consists in finding the solution vector f0(a) that, for a given KS11 (Θ1, a), produces the best match to the experimentally determined data vector RS(Θ1)/R 0ξ . The retrieval of the particle size number distribution function from Eq. 6 can be difficult because the problem is typically ill-posed until additional (a-priori) information on the aerosols can be applied. For instance, the material composition is one of the important properties that predetermine the range of applicability of the inverse transform. For urban, industrial-like aerosols with refractive index as described in (Kent et al. 1983) we have found that the peak contribution to the scattered signal is due to submicrometer-sized particles at almost all scattering angles θ (see Fig. 3a). On the other hand, the contribution from micrometer-sized particles to the scattered signal is also important if e.g. humidified sea salt is present in the atmosphere (see Fig. 3b). This implies that the method is neither sensitive to very small particles with size parameter x = 2πa/λ ∼ < 1 nor to particles larger than a few micrometers. The traditional solution to the Fredholm integral equation (Eq. 6) is to implement some suitable regularization algorithm, but the problem can be significantly simplified by reducing the number of unknowns. For instance, instead of searching for a N -element solution vector f0(ai) (i = 1..N ), an alternative and more efficient approach consists in using a parametric analytical model for the particle size number distribution function and determining its free parameters by minimizing a cost function, especially if the number of these unknown parameters is low. Ideally e.g. a log-normal distribution with two shaping parameters (the particle modal radius rM and the size-distribution width σ) may allow for transitioning from a complex inversion problem to a simple minimization one.

PROCESSING THE VIIRS-DNB NIGHTTIME LIGHT IMAGERY OF ZARAGOZA CITY
We present in this section a proof-of-concept of the application of this approach to obtain the columnar particle size number distribution functionf0(a) = ∞ 0 f (a, h) dh = HA f0(a) by inverting Eq. 6 with the kernel given in Eq. 7. We used as input data the radiances RS(Θ1) andR 0ξ measured by the VIIRS-DNB radiometer onboard the Suomi-NPP satellite during its pass over the city of Zaragoza (Spain,41.65 N,0.89 W) in the night of July 5th, 2019 at 03:59 CEST (01:59 UTC). This pass took place in a clear and moonless night, with the Moon at 30.14 degrees below the horizon, illuminated at 7.2 per cent, and Sun altitude -20.6 degrees, at the time of the measurements. That night The radiance as a function of the nadir angle α (measured in degrees, as seen from the satellite). The radiance has been taken at two edges for the first and the last contacts with the artificially lit surface area of the city. Red color is for the inner parts of the city.
the satellite passed almost vertically over the the city, at an altitude of 89 degrees over the horizon, corresponding to a zenith angle z0 = 1 • (see Fig. 1a, and ground track in Fig.  4a).
The Suomi-NPP satellite, launched on October 2001, is located in a polar sun-synchronous orbit of 826 x 828 km, with an inclination of 98.7 • , recording the Earth radiance in a panchromatic band (500-900 nm), with 14-bit quantization and a low light imaging detection limit of 5 × 10 −11 W cm −2 sr −1 (Cao and Bai 2014;Elvidge et al. 2013Elvidge et al. , 2017. The radiance data are available as geotiff files in the WGS84 coordinate reference system, with pixel size of 15 arcseconds, measured from the center of the Earth (NOAA 2019). This corresponds to a pixel angular size of about ∆Θ=115 arcseconds in the North-South direction, as seen from the satellite location at low nadir angles. The radiance pixel array analyzed in our study was taken along a North-South path passing through the urban center, which provides theR 0ξ values, and encompassing large unpopulated segments of territory at its extremes, where the recorded RS(Θ1) radiance is mostly due to atmospheric scattering. Zaragoza has a compact urban structure with well defined city limits, subtending, from the satellite, an angle of order 0.58 degrees. The average valueR 0ξ can be obtained from the radiance recorded within the city (see Fig. 4b, red line). The values of RS(Θ1) are measured from the city border, obtaining in this way two sub-arrays of measurements along the path, one before reaching the North border of the city and other starting from the South border (Fig. 4b, full black lines). The dashed black lines in Fig. 4b are logarithmic fits of the raw radiance, removing the small local peaks associated with some residual sparse sources of light.
Two instrumental data arrays can be immediately computed from these measurements: the ratio RS(Θ1)/R 0ξ , and the kernel KS11(Θ1, a) (according to Eq. 7). These arrays are the required inputs for performing the inversion of Eq. 6 in order to obtain the estimatef0(a) . This inversion can be carried out using different approaches. One of them is to perform a conventional linear estimation with an appropriate Tikhonov's regularization parameter. Another useful procedure consists on the choice of a suitable analytic form forf0(a), based on the available a-priori knowledge about the physical properties of the phenomenon under study, and determining the unknown parameters of this trial function by using conventional minimization routines. In this work we resorted to the latter option, by assuming that the particle size number distribution function can be well described by the linear combination of a set of i=1...M log-normal distributions in the form where each i-th distribution is characterized by a small set of parameters (ai, σi, Ni), informing about its position and width on the particle size space, and its relative weight in the linear combination, respectively. Ideally a low number of elementary distributions (M =1 or 2) should be enough for building up a good approximation tof0(a). Some additional a-priori knowledge may be used to improve the estimation. Zaragoza is located in the mid of a semi-arid zone, so that dust-like particles make an important part of the total aerosol content. The characteristic refractive index of dust-like particles, 1.52+0.025i, can be potentially influenced by humidity at night, so we used a reduced index of 1.4+0.002i. With these assumptions we minimized the quadratic differences between the measured radiance ratios, RS(Θ1)/R 0ξ , and those expected from Eq. 8 (substituted into Eqs. 6 and 7) with respect to the unknown parameters ai, σi, Ni (M =2), avoiding too narrow and too wide widths (σi), as well as to small and too large modal radii ai. The computations were made for λ = 650 nm, the weighted average wavelength of the VIIRS-DNB spectral sensitivity band. Incorporating all these constraints we obtain N1/N2 ≈ 125, a1 ≈ 0.08 µm, σ1 ≈ 0.15, a2 ≈ 0.5 µm, and σ2 ≈ 0.3. The resulting particle distribution function is shown in form of the volume density function (Fig. 5b) to be consistent with presentation form of the AERONET data products. Fig. 5a shows that the best match of the theoretical model (Eq. 4) to the experimental data occurs for AOD ≈ 0.3. The RS(Θ1)/R 0ξ ratio as a function of Θ1 becomes steeper for larger AODs, and more flat for lower AODs. In addition, the absolute values of RS(Θ1)/R 0ξ increase with increasing AOD, so one can easily estimate the optimum AOD for the time of measurement. These results are fairly similar to those reported by the AERONET station located in Zaragoza (41.63 N, 0.88 W, height=250 m) a few hours after the Suomi-NPP pass (NASA 2020): the recorded AOD at 05:45 UTC for λ=675 nm was 0.3, and the retrieved particle size volume distribution function clearly shows the same bimodal structure predicted in Fig. 5b, a feature of the aerosol distribution that was maintained along that day. The two peaks of that function at 06:07 UTC, first available data for this day, were centered in aerosol particle radii slightly above 0.1 and 2 µm, with maximum values close to 0.035 and 0.13 µm 3 µm −2 , respectively. The peak locations are similar, and the peak values are of the same order of magnitude but slightly lower than those deduced from the VIIRS-DNB images. Some differences could be expected, since the VIIRS-DNB measurements were obtained at 01:59 UTC, a few hours before dawn.
Let us point out that the actual particle volume distribution function is not expected to be as smooth as that displayed in Fig. 5b. The smoothness of this estimate comes from the fact that it has been obtained using a linear com- Figure 5. Left plot (a): Modelled (black lines) versus measured (green lines) ratios R S (Θ 1 )/R 0ξ for different aerosol optical depths at λ=650nm. The solid green lines correspond to the two experimentally determined data sets, while the dashed line is the weighted average of both. Right plot (b): The particle size distribution function determined by minimization of the differences between the model and the measurements and presented the same way as in AERONET database, i.e. dV /d ln a = a dV /da = (4/3)πa 4f 0 (a), where dV /da, where dV /da is the particle size volume distribution function andf 0 (a) = ∞ 0 f (a, h) dh = H A f 0 (a) is the columnar particle size number distribution.
bination of smooth log-normal elementary distributions, as defined in Eq. 8. A more precise estimate of this function could be obtained by inverting Eq. 6 using a suitable, but numerically more demanding, regularization procedure.

CONCLUSIONS
We demonstrate in this work that the diffuse atmospheric scattered radiance distribution around urban nuclei detected by on-orbit radiometers at night-time contains useful and retrievable information about some key aerosol properties, like the columnar number size distribution function. The method is self-calibrating and does not require the knowledge of the absolute city emissions. This opens the way for more comprehensive studies of the aerosol dynamics at night, improving our knowledge of their characteristics beyond the traditional aerosol optical depth. The results obtained with this approach applied to the radiance data obtained by the VIIRS-DNB radiometer in a particular pass over the city of Zaragoza are presented as a proof-of concept of its feasibility and performance. The DNB band is a panchromatic one, extending from the visible to the NIR. It can be anticipated that the use of multiband radiometry (e.g. RGB calibrated images of night lights acquired from the International Space Station) may provide additional information for a more detailed aerosol characterization.