Interactive comment on “ A new method of measuring aerosol optical properties from digital twilight photographs ” by

The paper is addressing an interesting subject concerning a method to infer aerosol optical properties from digital twilight photographs. The paper is interesting and could be accepted after minor revision that I am suggesting here. My reaction to the paper is that it would be good to compare the method proposed to a much simpler method described in Zerefos et al., 2007, 2014 where digital pictures were also used to calibrate the method before, during and after the overpass of a Sahara dust event over an experimental site. The method is simple and relies on calculating red-to-green ratios measured either in digital photographs, in situ and painted by professional colorists. Therefore I would like to see a short discussion on the methods. It appears that large volcanic eruptions have left a high red-to-green ratio either in digital pictures or in paintings by great masters in the past. With this minor correction I think the paper can be


Introduction
Twilight sky, one of the most beautiful sights in our daily life, varies from day to day.The color of the sky under clearsky conditions gradates from orange-red and salmon red near the horizon to blue-gray and blue at zenith.This gradation results from Rayleigh scattering caused by molecules (Minnaert, 1993).Twilight sky has been studied for more than a half century, elucidating the effects of ozone and stratospheric aerosol on twilight sky colors.As an example, ozone has the Chappuis band, an absorption band with a peak wavelength of about 550 to 600 nm.Hulburt (1953) simulated twi-light sky brightness at zenith with a radiative transfer model (RTM) based on the Rayleigh scattering theory and taking into account ozone absorption.Hulburt's results showed that the contribution of ozone absorption is more important than that of Rayleigh scattering in providing the blue color at zenith. Lee et al. (2011) measured twilight sky in Antarctica with a digital still camera to investigate ozone effects on twilight sky at zenith and the anti-solar horizon, with the results suggesting that scattering light by aerosols is as important as ozone absorption in explaining the color of twilight.These results are consistent with simulation results by Adams et al. (1974).
In 1963, Mt.Agung in Bali, Indonesia, erupted and injected a massive amount of volcanic ash into the stratosphere.After the eruption, stratospheric aerosol effects on twilight sky have been investigated by many researchers.Volz (1964Volz ( , 1965) ) suggested that a persistent purple light in the twilight term seen during 1963 was caused by an increase of stratospheric aerosol as a result of the eruption, using twilight sky measurements and geometrical calculations to support this.Dave and Mateer (1968) demonstrated that stratospheric aerosol was more important than ozone as a mechanism for creating purple light, showing this by simulations with RTM taking into account a spherical-shell atmosphere, dust, and ozone effects.
The general radiative transfer process for radiance at a large viewing zenith angle (VZA) in the twilight sky is more complicated.Because of this, the effects on twilight sky color near the horizon of atmospheric components have not been investigated.To properly do so, an RTM must incorporate a spherical-shell atmosphere, multiple scattering, and refraction effects.Shaw (1981) investigated the dependence of Published by Copernicus Publications on behalf of the European Geosciences Union.stratospheric aerosol vertical profiles on twilight sky at VZA of 70 • by using a single-scattering analytical RTM with a spherical-shell atmosphere.Recently, an RTM for calculating spectral radiation according to a complicated radiative transfer process has been developed (Iwabuchi and Suzuki, 2009); this model achieves realistic simulation of twilight sky near the horizon.Saito et al. (2013) showed that tropospheric aerosol optical thickness (AOT) was sensitive to twilight sky brightness near the horizon by analyzing photographic data and RTM simulations.Those results suggest that multiply scattered light is dominant in twilight sky color and brightness near the horizon.
Some researchers have tried to retrieve aerosol information from twilight sky data.Wu and Lu (1988) developed a method to infer stratospheric aerosol vertical profiles from twilight sky brightness and degree of polarization measurements at a wavelength of 700 nm at zenith.Mateshvili et al. (2013) developed a method that uses optimal estimation to infer the upper tropospheric and stratospheric aerosol vertical profiles from spectral at zenith during twilight, providing stratospheric AOT both pre-and post-eruption of the volcano Nabro.Zerefos et al. (2014) tried to retrieve AOT using a color ratio (R/G), given by digital counts in the red, green, and blue (R, G, and B) channels, from digital photographs and digitalized paintings drawn during the time around sunset in the past hundreds years.The results agree well with other ground-based observations.The aerosol information obtained from these twilight measurements was consistent with that from other measurements, showing that twilight sky measurement is a useful tool for gathering information on aerosols.
It is well known that tropospheric and stratospheric aerosol contributes to climate change.The optical and radiative properties of aerosol can vary strongly according to relative amount, chemical components, and particle sizes.This results in high uncertainty of climate change prediction in terms of radiation budget (IPCC, 2013).Aerosol measurement is carried out with satellite-and ground-based instruments such as the MODerate-resolution Imaging Spectroradiometer (MODIS) devices on-board the Terra and Aqua satellites and the AErosol RObotic NETwork (AERONET) system.This provide on solar and sky measurements (Remer et al., 2005;Holben et al., 1998).However, most of the instruments useful for aerosol measurement require passive remote-sensing techniques, meaning that these instruments cannot infer aerosol properties in regions without direct sunlight.In addition, satellite-based instruments with passive sensors cannot infer aerosol properties in deserts and snowcovered regions because of high surface reflectivity.Particularly in the polar region during winter, aerosol measurements and their data are very limited because of the length of the polar twilight (sunless condition in all day).In the region, a few measurements are carried out by active remote-sensing techniques, in situ sampling, and difficult aerosol measurements that use direct light from the moon and stars as a sub-stitute for direct sunlight transmittance methods (Herber et al., 2002;Berkoff et al., 2011).
Recently, digital cameras have begun being used as measurement instruments (e.g., Lee and Andres, 2003;Lee et al., 2011).Ehrlich et al. (2012) provided a method for constructing a bi-directional reflectance distribution function for cloud from data obtained by a digital single-lens reflex (DSLR) camera installed on an airplane.Kataoka et al. (2013) developed a method for estimating the altitude of visible aurorae from data obtained by a pair of DSLR cameras.Hioki and Iwabuchi (2014) demonstrated a method to find the particle radius of pollen and its column number density from a solar corona image taken with a digital still camera.A two-dimensional brightness distribution can be obtained by digital camera, an ability that is likely to be useful in many scientific fields.The aim of the present study is to infer aerosol optical properties, such as tropospheric and stratospheric AOT and particle size information, from twilight photographs taken with a DSLR camera.This paper is organized into six sections.Section 2 presents the photographic measurement and camera characteristics.Section 3 shows the forward model and the tests to determine the sensitivity of aerosol optical properties to twilight sky.In Sect.4, the method of aerosol retrieval and the error analysis are explained.Section 5 discusses the retrieval results and presents some case studies.The results are compared with those inferred from other ground-and satellite-based instruments.Section 6 summarizes the study and discusses a few points of interest from the retrieval results.

Measurements from a digital camera
The instrument for photographic observations in twilight was a Nikon D7000 DSLR camera equipped with a Nikkor equisolid angle-type full-frame fisheye lens (AF DX Fisheye-Nikkor, 10.5 mm f/2.8GED).The camera configuration simulated an International Organization for Standardization (ISO) film speed of 200 and had an aperture of f2.8.These values were set low enough to reduce complementary metal oxide semiconductor (CMOS) noise and to gain high exposure, even in the dark conditions of twilight.The camera was fixed on a tripod.Table 1 summarizes the measurement setup and the camera characteristics.Photographic observations of twilight sky have been carried out at the Graduate School of Science, Tohoku University (38.26 • N, 140.84 • E; 208 m elevation), Sendai, Japan, since January 2013.Various instruments for aerosol measurements are collocated at the site.These include a Mie scattering lidar system operated by the National Institute for Environmental Studies, a sky radiometer operated by SKYNET and a sun photometer (PREDE PGS-100) originally managed by Tohoku University (Shimizu et al., 2010;Takamura and Nakajima, 2004).The twilight observation criteria are (1) the absence of clouds at the zenith during the twilight, measured by the Mie scat-  tering lidar, and (2) the absence of clouds over a 200 km area on the solar side of the observation site, measured by infrared imagery from the Multifunctional Transport Satellite (MTSAT), a geostationary meteorological satellite.Table 2 describes the VZAs, solar zenith angles (SZAs), and relative azimuth angles (RAAs).The RAAs are the differential angle between the viewing azimuth and the solar azimuth for data used in aerosol retrievals derived from photographic observations.The measurements are carried out for an SZA ranging from 90 to 96 • in intervals of 1 • .The camera exposure time is changed according to SZA so as to be suitable for imaging.To obtain the color and brightness distribution of a twilight sky from a RAW image file, we first extract the digital numbers for the three spectral channels (RGB) in groups of 12 × 6 pixels at corresponding angles (see Table 2 for details).After that, the averaged digital counts for each channel are used to eliminate pixel-level CMOS noise.The RAW format data provide a wider dynamic range than the 8bit JPG format does, which is sufficient to provide brightness information.To read the Nikon D7000 manufacturer-specific RAW format (for the Nikon D7000 the exact format is NEF), we employed the useful open-source tool DCRAW.DCRAW has been described in Ehrlich et al. (2012).Figure 1 shows twilight photographs taken under typical atmospheric conditions.
The color matching functions (CMFs) of the camera for each RGB channels were determined at the National Institute of Polar Research (NIPR) to quantify the spectral radiation from the digital numbers.The details of the method for determining the CMF for a specific camera are given in Saito et al. (2015).This paper briefly describes the methods.The experiment setup is described in detail in Sigernes et al. (2008).The camera and spectroradiometer detect narrowband (i) light in visible wavelengths of 380-740 nm at intervals of 5 nm, which are separated by a grating monochromator, obtaining camera RAW counts u k,i and the corresponding spectral radiance I i (λ).The camera spectral chan- nels and wavelengths are defined as k and λ, respectively.To enable comparison between the values, the values of u k,i are normalized by the photon flux density convert radianceproportional RAW counts u k,i and u k,i , which can be defined as where the function S k (λ) is the relative CMF.The optimalestimation-based inverse calculation is applied to Eq. (1), yielding determining S k (λ) with some associated uncertainty.In this experiment, the same camera configurations are used for all photographic observations, and an exposure time of 8 s is used.The CMFs for the Nikon D200, as determined by Sigarnes et al. (2008), are used as prior information.The vignetting characteristics of the fisheye lens have been investigated by Saito et al. (2015).We tried to eliminate the effects of vignetting from the image digital counts by an experiment using integrated sphere and Halogen light source that enables to generate homogeneous light.In the experiment, the intensity of integrated light at an edge is usually lower by 1 % than that at the center so that the vignetting correction parameter has such biases near the edge of photographs, which does not significantly influence the aerosol retrievals.Therefore, we assume that the vignetting effect is totally eliminated for our purposes.Figure 2 shows the relative CMFs of the Nikon D7000.The center wavelength of each channel (i.e., the most sensitive wavelength) is shifted to a shorter wavelength, by about 10 nm for R channel and 20 nm for the B channel, relative to those of the Nikon D200.Moreover, the shape of the spectral sensitivity on the G channel of the D7000 is broader in shape than that of the D200.
For photographic observation in twilight, exposure times are optimized to sky brightness levels, which make the digital counts on the pixel used for the retrievals in a dynamic range of the camera.For example, we used the exposure time of 1/2048, 1/512, and 1/4 s in the SZAs of 90, 93, and 96 • , respectively, in clean atmospheric conditions.The RAW counts u k are converted to radiance-proportional RAW counts u k , which can be rewritten as cI k , where I k is channel-averaged spectral radiance weighted by the camera spectral response of the channel k.The parameter c is the absolute calibration factor that can be defined as independent of wavelength, because the wavelength-dependence of CMOS detector is considered in the CMFs.Finally, the color ratios (R/G, B/G) are obtained by dividing each of cI R and cI B by cI G , and the normalized brightness G n is obtained by dividing the value of cI G by the value of cI G measured at VZA of 70 • and RAA of 0 • for all VZAs, SZAs, and RAAs in the twilight sky.This process results in an 840-dimensional measurement for each photographic observation.For prior knowledge, this method requires only the relative CMFs of the camera; absolute calibration is unnecessary because the unknown parameter c can be excluded from the data analysis.Absolute calibration can be used to estimate the parameter c, but this is beyond the scope of the current paper and is left for future work.

Model setup for twilight simulations
For twilight sky simulations, we use JACOSPAR, a spherical-shell atmosphere RTM that accounts for refraction and multiple scattering (Iwabuchi and Suzuki, 2009).The RTM calculates singly scattered light by a semi-analytical approach and multiply scattered light by a backward Monte Carlo (BMC) method.The RTM also considers 300 wavelength bands in the range from 370 to 740 nm and models gaseous absorption by the correlated k-distribution method used in the RSTAR RTM (Nakajima and Tanaka, 1986;Sekiguchi and Nakajima, 2008).The earth's surface is assumed to be a Lambertian surface with an albedo of 0.1.In the model atmosphere, we assume the presence of 70 atmospheric layers, from ground to an altitude of 120 km, with geometrically thin layers, particularly in the lower part of the atmosphere.For simulating their scattering properties at each wavelength band, aerosols are assumed to be spherical particles and the effects are calculated according to the Lorenz-Mie theory (Bohren and Huffman, 2008).The normalized vertical profile of stratospheric aerosol extinction coefficient b st (z) at a wavelength of 550 nm is assumed to follow a Gamma distribution b(z), defined as where (ξ ) denotes the Gamma function, and ξ and ψ are parameters that depend on the shape and scale, respectively.The normalized vertical profiles of the tropospheric aerosol extinction coefficient b tr (z) are assumed to be the sum of the vertical profiles of background aerosol b bg (z), aerosol in boundary layer b bl (z), and upper tropospheric aerosol b up (z), each of which follows a Gamma distribution.The shape and scale parameters for each aerosol vertical profile are summarized in Table 3.By substitution, the tropospheric aerosol vertical profile b tr (z) can be written as where a bg , a bl , and a up are, respectively, the AOT fraction of background aerosol, boundary layer aerosol, and upper tropospheric aerosol relative to total tropospheric AOT at a wavelength of 550 nm.In this study, the AOT fractions are assumed to be in the proportions a bg : a bl : a up = 1 : 1.2 : 0.4.Figure 3 shows a priori aerosol vertical profile.In general, the extinction coefficient in the middle-upper troposphere shrinks exponentially with increasing altitude, and aerosol is abundant in the boundary layer.The extinction coefficients β(z) can be calculated by multiplying the AOT by b(z).The distribution of aerosol particle volumes dV /dr is assumed to follow a lognormal distribution: Here, r is particle radius, and r mod , s, and V 0 are the mode radius, geometrical standard deviation, and particle volume, respectively, with subscript F indicating fine particle mode and subscript C indicating coarse particle mode.The parameters r mod and s for this aerosol model are taken from the rural aerosol model proposed by Hänel (1976).To simplify for aerosol properties representation, this paper introduces the coarse-fine particle volume ratio ς , defined as From the above model setting, the spectral radianceintegrated RGB signals can be calculated by multiplying the simulated spectral radiance of twilight sky with the camera CMFs.The normalized brightness G n and the color ratios (R/G, B/G) can then be obtained in the way described in Sect. 2.

Sensitivity tests
The sensitivities of the color ratios and normalized brightness in twilight sky to tropospheric AOT τ tr , stratospheric AOT τ st , coarse-fine particle volume ratio ς , and other atmosphere-surface characteristics were tested.For the sensitivity tests, the complex refractive indices of aerosol were assumed to be 1.5 for the real part and 0.01 for the imaginary part.The results are described in Figs. 4, 5, and 6.In the twilight sky, and particularly near the horizon, as τ tr increases, G n and R/G decrease and B/G increases (Fig. 4a,  b).These relations suggest that the twilight sky becomes dark and blue near the horizon, which is consistent with the results of Saito et al. (2013).Figure 4c and d show the sensitivities of the twilight sky to τ st .With large τ st , the value of G n at VZAs of 75-85 • increases and all R/G (B/G) values increase (decrease), so that the reddened twilight glow is represented.These results are consistent with several reports on the very intense twilight glow (the so-called twilight phenomenon) that results from volcanic aerosol injection into the stratosphere by volcanic eruptions (Volz, 1964(Volz, , 1965)).The sensitivity tests show that R/G increases under high ς conditions only near the horizon (Fig. 4e, f). Figure 5 shows the sensitivities of the brightness and the color ratios in twilight glow with VZA of 85 • , SZA of 93 • , and RAA of 0 • .In the R/G-G n cross-section (Fig. 5b), G n depends strongly and inversely on τ tr , decreasing under typical aerosol property condition (ς = 1) as τ tr increases.In contrast, for R/G, sensitivity to τ st predominates and R/G increases with large τ st (Fig. 5b).The sensitivities of twilight glow are sufficient for modeling with τ tr in 0.01-0.5 and τ st in 0.002-0.1,while the color ratios and normalized brightness in this sector are less sensitive to ς (Fig. 5c, d). Figure 6 is the same as Fig. 5 except near the horizon with VZA of 88 • , SZA of 91 • , and RAA of 0 • .The sensitivity of the twilight sky to τ st and ς is weak for τ tr = 0.13 (Fig. 6a, b). Figure 6c and d suggest that the sky is more sensitive to ς when τ tr is smaller.In such cases, R/G and G n increase with increasing ς .As above, twilight sky is sensitive to ς in the range of 0.2 to 5 under the clean atmospheric condition (τ tr <∼ 0.1).To summarize this section, the sensitivity tests showed that the normalized brightness and the color ratios in the twilight sky have sufficient sensitivity to tropospheric and stratospheric AOT and have moderate sensitivity to the coarse-fine particle volume ratio under small AOT conditions.
The sensitivity to other atmospheric and surface characteristics was investigated, with the results showing that the sensitivity to ozone column of twilight sky color is moderate, as previously shown by Hulburt (1953).However, the sensitivities of sky brightness and color to surface albedo, refractive indices, and water vapor column were found to be very weak (not shown).

Model setup for aerosol retrievals
To apply the JACOSPAR RTM to the forward model for aerosol retrieval, we refined the RTM to allow faster computation by semi-analytically calculating singly scattered light and then linearly interpolating multiply scattered light from a lookup table (LUT), which was pre-calculated by the BMC method.Ozone column amount data were taken with the ozone monitoring instruments (OMIs) on board the Aura satellite and considered in the atmospheric model and in the LUT of the RTM.We assumed that other surfaceatmospheric conditions (such as: surface albedo, water vapor, and refractive indices) agreed with climatic values.

The retrieval method
The retrieval method is a maximum a posteriori estimation method based on optimal estimation (Rodgers, 2000).This method derives an optimal solution from measurement data under the constraints of prior information and is presently in use for aerosol and cloud optical property remote sensing (e.g., Watts et al., 2011;Mateshvili et al., 2013;Iwabuchi et al., 2014).A measurement vector y is derived from twilight photographs taken with SZAs of 90-96 • .The color ratios (R/G, B/G) and the normalized brightness (G n ) around solar direction (VZAs of 60-88 • and RAAs of 0-30 • ) were used in the elements of the measurement vector.A state vector x that includes aerosol optical properties as elements, and a set of model parameters p are defined as where  aerosol.The number of measurement vector elements n is set at 840 because the number is defined as the number of measurement elements multiplied by the number of measurement directions (3 of measurement elements and 280 of measurement directions, which result from 10 VZAs, 7 SZAs, and 4 RAAs).We can formally describe the physics as the forward model function f (x): where e is the measurement-model error vector.
The cost function J (x) can be described as where x a is an a priori state vector, S a is an a priori state error covariance matrix, and S y is a measurement-model error covariance matrix.The derivation of S y will be described in Sect.4.2.We summarize the a priori state and its error variance in Table 4.The assumed S a is a diagonal matrix, with the main diagonal taken as a vector of error variances from prior information.The prior state and its geometrical standard deviation for τ tr are assumed to be 0.13 and 1.19, respectively, which is taken from the values in Dubovik et al. (2002) for AOT in urban, continental, and desert regions, which range from 0.04 to 0.43, as found from ground-based observations.The a priori value of τ st is set to be 0.005 with a geometrical standard deviation of 1.15 (Hess et al., 1998).
The a priori value for ς is assumed to be that of a rural  aerosol model with a bimodal distribution of particle volume, as proposed by Hänel (1976).The geometrical standard deviations for the prior information are set to be large value to avoid overfitting of the state vector x, which has components ln τ tr , ln τ st , and ln ς.When the state vector becomes optimal, the cost function should be below the degrees of freedom of the measurement vector (Rodgers, 2000).Therefore, an optimal state x opt gives Because the measurement vector calculated according to the forward model is nonlinearly dependent on the state vector, we apply the Levenberg-Marquardt method (LMM) to solve the inverse problem.The LMM can provide a stable solution for moderately nonlinear problems (Rodgers, 2000).The method is iteratively applied to find an optimal solution for x.The state vector at the i + 1th iteration, x i+1 , is given by where K i is the Jacobian of the ith iterative calculation, and γ represents an adjustment adaptively chosen at each step so as to minimize the cost function.This construction lets the LMM avoid convergence at only local minima more readily than the Gauss-Newton method can.When the cost function is not sufficiently minimized to m, the forward model is unable to represent the measured twilight sky with the available data under the given assumptions and a priori constraint in the retrieval method.This can occur under two conditions: when clouds, fog, and/or haze is present; and when the values of the state vector and model parameters are out of range for the retrieval method, which results in underestimation of the measurement-model error covariance.Simultaneously with obtaining an optimal solution, the uncertainty can be calculated by the following equation: where S x is an error covariance matrix whose diagonal elements represent the error variances of the solutions.

Measurement-model error
Measurement-model error originates from the measurement error, assumed errors in the forward model, and uncertainties about model parameters.The matrix S y combining these three error components can be found as where S y,m , S y,fwd , and S y,p are the covariance matrices of the measurement error, the assumed errors in the forward model, and the uncertainties about the model parameters, respectively.These errors and the assumptions of the forward model are listed in Table 5.The real (m r ) and imaginary (m i ) parts of the refractive indices are assumed to be 1.5 and 0.01, respectively, which are taken from Hayasaka et al. (1992) by simplifying from the annual variation of refractive indices for aerosol optical properties in Sendai, Japan, as obtained by ground-based measurements (1.44-1.57real part and 0.007-0.057imaginary part).The surface albedo α is assumed to be 0.1 because of the annual variation of α from 0.05 to 0.3, which was found by Zhou et al. (2003).The forward model error originates from uncertainties about the particle size distribution and from the vertical profile shape assumed in the forward model, and it also includes a smoothing error resulting from linear interpolation in the LUT.The forward model error is evaluated by adding Gaussian random noise to the parameters for aerosol particle size distribution and vertical profile shape (specifically to the mode radii for the volume particle size distributions and to the scale and shape parameters of the Gamma distributions for the vertical profile) and running the forward model and JACOSPAR RTM for perturbed sets of model-assumed parameters and finally calculating the difference between each model's output.In general, aerosol vertical profiles are highly variable in the troposphere.The amount and vertical profile in the boundary layer varied according to atmospheric condition, and those in upper troposphere are also variable as typically seen in dust transport events.Considering their variability, we determined the possible ranges of the shape and scale parameters in the forward-model-error analysis.The results from JACOSPAR RTM are taken as ground truth data.The dependence of assumption errors in the forward model on the measurement angles is shown in Fig. 7.The root-meansquare (RMS) errors of the color ratios and the normalized brightness are consistently larger with larger SZAs and are The model parameter errors are evaluated by adding Gaussian random noise to the model parameters shown in Table 5 and running the forward model for perturbed sets of model parameters.Aerosol optical properties, such as the refractive indices and the coarse-fine particle volume ratio, are assumed to be vertically homogeneous through the entire atmosphere from the stratosphere to the troposphere.In general, aerosol optical properties in the troposphere and the stratosphere are essentially different.Most of stratospheric aerosol is sulfate or sulfuric acid, while tropospheric aerosol is highly variable.However, aerosol amount in the troposphere is dominant compared to that in the stratosphere under the present normal conditions.The retrieved optical property, namely coarse-fine particle volume ratio, is representative as a tropospheric one.Furthermore, the twilight sky color and brightness are not very sensitive to the refractive indices of aerosol according to a sensitivity test, suggesting that the assumption of constant refractive indices leads to insignificant uncertainties.We show the dependence of the model parameter errors on the measurement angles in Fig. 8.The color ratio errors depend weakly on the value of RAAs and SZAs and are moderately larger with larger values of VZAs.In con-trast, the error in G n rapidly increases with increasing SZAs and is slightly larger with smaller values of RAAs and larger values of VZAs.The RMS errors for R/G and B/G are in the range 3-5 %.The RMS error for G n is below 5 % for SZAs of 90-94 • but in the range 20-80 % for SZAs of 95-96 • .The model parameter errors are the largest errors and are dominant in the model-measurement error.The covariance matrices for errors in the forward model assumptions and the model parameters are taken as constant matrices for the aerosol retrievals in this study.The applied algorithm could be improved by accounting for the effects of the state vector and model parameter values on these matrices.This is left for future work.
The measurement errors include CMOS noise and biases, uncertainties in the inferred camera response functions, camera-orientation errors, and time-lag error while taking the photographs.The dark current also leads relatively small bias in digital counts in the twilight photographs used in this study: less than 0.5 % (0.1 %) for the camera signal corresponding to the twilight sky at VZA of 60 • (88 • ).CMOS biases can arise from exposure time, aperture, and temperature.CMOS noise errors for each angle are determined by the variance of RAW counts among pixels taken from the corresponding angles.The time-lag error is assumed to be normally distributed with standard deviation 5 s, which gives an uncertainty of ∼ 0.015 • for SZA, and the camera-orientation error is assumed to be normally distributed with standard deviation 0.5 • for the VZA and 0.3 • for the RAA.The inclined camera orientation leads to biases in the normalized brightness and color ratios especially in pixels with large RAAs and small VZAs even when we use a spirit level.Uncertainty of orientation angle is approximately 0.8 • .Moreover, a stray light, which may affect the camera digital counts, is caused by multi-reflected light among the lenses and reflection on near-Lambertian surface, namely a plastic part installed in the equipped lens.The dynamic range of twilight sky brightness is large, and bright light might cause a stray light in the dark part of the photograph.We have experimentally investigated the order of magnitude of a stray light effect on digital counts in twilight photographs.It was found that the stray light effect provided small biases, which became significantly smaller towered to larger angle from a bright light: approximately 1, 0.1, and below 0.05 % at angle of 1, 5, and > 15 • , respectively, from the light source.Furthermore, the above errors are estimated by running the forward model with Gaussian random noise perturbation of the camera angles.As a result, the CMOS noise is on the order of 3-5 %, and other errors including lens characteristics and angle uncertainties are on the order of 5 % for the color ratios and 10 % for G n .The distortion of a fisheye lens may lead to uncertainties in VZA and RAA estimation.Currently, several researchers investigated and showed that the recent lenses for consumers had very slight distortion, leading to a difference between ideal and distorted optical path within a pixel (Schneider et al., 2009;Shahriar et al., 2009).Rabaza et al. (2010) showed that the distortion effect for a Nikkor fisheye lens was considered negligible.For this reason, we assumed the lens distortion did not influence our results.

Retrieval-error analysis
The errors in retrieved aerosol properties originate from the measurement-model error and the sensitivity of the state vector of the sky to atmospheric conditions.For the error sources in Table 3, the retrieval error is evaluated by retrieval simulations.First, after assuming a state x and model parameters p noise and adding Gaussian random noise, the forwardmodel-simulated measurements y noise are calculated, and Gaussian random noise is added to simulate measurement error and the forward model error.Second, we try to retrieve the state vector x noise from y noise with this algorithm.Finally, the truth-retrieval difference x noise -x is obtained.Each state was simulated 200 times and the retrieval error was taken as the average error.The retrieval results with filtering (J (x) < 1000) are shown in Fig. 9.The value of τ tr is accurately retrieved, with RMS error of 54.5 % (absolute RMS error of 0.025) in the range of 0.01-0.5, and positive biases are found particularly under small AOT conditions (Fig. 9a).The retrieval uncertainty is small under the clean and moderately turbid atmospheric conditions for τ tr from 0.01 to 0.3, suggesting that this method can retrieve tropospheric AOT with high precision.In the more turbid atmospheric conditions (τ tr > 0.3), the uncertainty in the retrieved τ tr is around ±0.15.The value of τ st is also correctly retrieved with a RMS error of 53.2 % (absolute RMS error of 0.005) and slightly positive bias (13 %) (Fig. 9b).The retrieved value of τ st is overestimated under small tropospheric AOT conditions as a result of the particular characteristics of this retrieval method (not shown).The uncertainties for τ st retrieval are small for large values of τ st and are slightly larger for small values of τ st , suggesting that the value of τ st can be retrieved under the dense τ st conditions caused by volcanic eruption, such as the Mt.Pinatubo eruption in 1991.
From Fig. 9c, it seems difficult to retrieve the value of ς under typical atmospheric conditions and under fine-particledominance conditions.The RMS error is large, and the retrieved values of ς are scattered, centered mostly around the a priori assumption even under clean atmospheric conditions.This occurs because the sensitivity of ς on twilight sky is too low to overcome the measurement-model error.When coarse particles dominate in the atmosphere (ς > 1) under clean atmospheric conditions, the value of ς might be retrieved, and the RMS error and bias are 62.7 and 8.8 %, respectively.In addition, the averaging kernel shows a sensitivity of retrievals to the true states, demonstrating whether the measurement with associate uncertainties contains enough information to retrieve the state vector or not (Rodgers, 2000).The diagonal elements of the averaging kernel are 0.92, 0.74, and 0.94 (0.94, 0.94, and 0.95) corresponding to tropospheric AOT, stratospheric AOT, and coarse-fine particle volume ratio under the atmospheric condition with tropospheric AOT of ∼ 0.1 (0.3), while the nondiagonal elements are relatively significantly small (below 0.05) except for the sensitivity of coarse-fine particle volume ratio on the true stratospheric AOT state (around 0.3) under small tropospheric AOT condition.Furthermore, the degree of freedom for signal given by trace of the A shows above 2.5 under the general atmospheric conditions.Therefore, it is reasonable to retrieve the three state variables under clean atmospheric condition and at least tropospheric and stratospheric AOT when AOT is large.

Validation
To validate this retrieval method, the retrievals from photographic observations were compared with those from the sky radiometer measurements located at the same site and with satellite-based measurements (Fig. 10).The comparison is limited to selected results with J (x) < 3000.The sky radiometer and an the software tool Skyrad.packdeveloped by Nakajima et al. (1996) were used as the measurement system and tool to infer aerosol optical properties, such as AOTs and distributions of particle volumes, from direct and diffused sunlight.For twilight, sky radiometer measurements are not available because of the absence of direct sunlight.Typically, aerosol optical properties do not vary much within a few hours except when a strong emission source (e.g., forest fire) is present near the site.For this reason, we used the sky radiometer data obtained 1 h before sunset and 1 h after sunrise for comparison in this study.Figure 10a shows the comparison of AOTs (τ tr + τ st ) retrieved from the photographic observations and the sky radiometer measurements.Mostly, the retrieved AOTs from the digital camera were consistent with those from the sky radiometer, with little uncertainty and a correlation coefficient of 0.94.The values of τ st retrieved from the camera data were compared with those from the Canadian optical stratospheric and infrared imaging system (OSIRIS) on board the Swedish satellite Odin (Murtagh et al., 2002;Llewellyn, et al., 2004).OSIRIS measures atmospheric limb radiance of scattered sunlight across visible and near-infrared wavelengths as well as optical thickness, and the vertical profiles of the extinction coefficient at 750 nm in stratospheric aerosol were inferred by the use of the algorithm given by Bourassa et al. (2007).The photographic observations are carried out under the background stratospheric aerosol conditions, in other words, without any effects of volcanic eruptions in the stratosphere.Figure 10b shows the scatter plot of τ st retrieved from the camera and from OSIRIS.Note that we compare stratospheric AOTs at 750 nm from OSIRIS and those at 550 nm from the photographic observations.For stratospheric aerosol particle size, the Ångström coefficient is generally in a range of 0-2, resulting in the stratospheric AOT at 750 nm being smaller than that at 550 nm by 0-26 %.It is shown that average stratospheric AOT from the photographic observations is in the  same order of magnitude but positively biased by 10-50% compared with those from OSIRIS, and no correlation exists for this measurement period.The variability is in the same order of magnitude compared with the RMS error in τ st retrievals from photographic observations according to the retrieval-error tests.The values of ς retrieved from the camera are compared with those calculated from particle volume distributions retrieved from the sky radiometer (Fig. 10c).
The retrieved value of ς tends to be close to the a priori assumption, and thereby no correlation between the camera and the sky radiometer is found.

Case studies
To test the flexibility under various aerosol conditions (e.g., vertical profiles and aerosol species) of this retrieval method, we chose three characteristic cases and analyzed the simula-tion results in detail with respect to the vertical profiles of the aerosol extinction coefficients.Figure 11 shows photographs of the twilight sky taken on (a) 31 January 2013, (b) 23 May 2013, and (c) 3 June 2013.These photographs were taken at an SZA of 94 • .In Fig. 11a, a bright red-orange twilight glow is seen near the horizon, and the twilight sky color gradates from orange to blue.The brightness rapidly decreases as the VZA decreases.In contrast, the color in Fig. 11b seems to be dull across the sky in comparison with Fig. 11a.An orangeto-gray twilight glow is present above a dark-blue segment that appears near the horizon (Saito et al., 2013).In Fig. 11c, the range of the twilight sky color is the smallest among the three cases, and a dark-blue segment can be seen near the horizon, similar to that in Fig. 11b. Figure 12 shows the vertical profiles of aerosol extinction coefficients in the troposphere from the ground to an altitude of 6 km; these are derived from the lidar measurements for each case (Shimizu et al., 2010).In the lidar products, the lidar ratio of 50 sr is assumed, and the inferred extinction coefficients are separated into spherical and nonspherical particles.Table 6 summarizes the results from the ground-based observations (the sun photometer and the sky radiometer, collocated at the photographic observation site) and from the daily means of aerosol optical properties, as calculated by the aerosol transport model SPRINTARS (Takemura et al., 2000(Takemura et al., , 2002(Takemura et al., , 2005)).The sun photometer measures direct solar light in the visible to near infrared wavelengths to infer the AOTs.From the above results, the atmospheric characteristics are clean atmosphere in case (a), dust-pollutant mixed particles transported from eastern China in case (b), and yellow sand transported from inland China in case (c).Additionally, the vertical profiles and aerosol species are vertically inhomogeneous in cases (b) and (c).A comparison of the results retrieved from the photographic observations with the corresponding results from other products shows that the retrieved AOTs are reasonable but slightly lower than the AOTs ob-tained from other products for the three cases with AOTs of 0.056 ± 0.029 (a), 0.091 ± 0.006 (b), and 0.261 ± 0.007 (c) relative to estimates from other ground-based passive instruments.The retrieved values of ς are 2.52 ± 0.34 (a) and 0.99 ± 0.11 (b).The retrieved value of ς is eliminated in the case (c) because of low confidence in the ς retrieval, which was predicted by the retrieval-error analysis (Sect.4.3) for large AOT conditions.The main reason for the under estimation might be the temporal-spatial variability of aerosol optical properties.Specifically, during dust events, aerosol optical properties showed high variability, and so they varied more widely during the time lag between camera observations and other measurements than they did in the other cases, resulting in the difference in AOTs.The AOTs for those three cases are consistent with the total AOTs calculated by SPRINTARS.Thus, the case studies demonstrate that the proposed method can infer AOTs without much dependence on knowing the vertical inhomogeneity for aerosol species and the extinction coefficients.

Conclusions and discussions
An optimal-estimation-based algorithm was developed to infer aerosol optical properties from twilight photographs taken with a digital camera at SZAs of 90-96 • .The following advantages are offered by this method and measurement techniques: (1) tropospheric and stratospheric AOTs can be retrieved by a passive remote-sensing technique, even in the absence of direct sunlight conditions; (2) the retrievals and their uncertainties are simultaneously calculated in the optimal-estimation framework; and 3) a low-cost instrument (a digital camera) can be used.
The measurement vector was set to be the color ratios and the normalized brightness around the solar direction derived from camera RAW data in RGB channels.Ozone column data were taken from OMI/Aura.The forward model includes a spherical-shell atmosphere, refraction, and multiple scattering.Additionally, in the model, singly scattered radiances are calculated semi-analytically, and multiply scattered radiances are calculated by linear interpolation from a pre-calculated LUT.The values of diagonal elements in the a priori state covariance matrix were assumed to be large in or-der to avoid overfitting.The model parameters and their uncertainties were determined by prior measurement and error analysis.The error sources during retrieval were chiefly the model parameter uncertainties, with each of CMOS noise, the assumed distribution of particle volume, and vertical profile shape contributing some error.The retrieval-error analysis showed that the tropospheric and stratospheric AOTs can be inferred for τ tr in the range 0.01 to 0.5 and τ st in the range 0.001-0.1,obtaining a small RMS error (around 50-55 % for τ tr and τ st ).Especially stratospheric AOT retrieval can be outperformed under dense τ st conditions.The retrieved AOT from the digital camera was compared with those from the sky radiometer measurements, and this comparison showed consistency between the results, with a correlation coefficient as large as 0.94.The retrieved stratospheric AOTs from the digital camera were positively biased compared with that from OSIRIS measurements with reasonable variability as shown in the retrieval-error analysis.No correlation was found in the ς comparisons.A case study for three different atmospheric conditions was presented to demonstrate the usefulness of the method.Comparison with products provided by other instruments and results from other model was used to investigate the effects of vertical inhomogeneity for aerosol species and profiles in detail.This investigation showed that AOTs can be retrieved with reasonable accuracy without knowing the details of such properties.
The difficulty in retrieving the coarse-fine particle volume ratio ς can be explained by the following.In cases with a large AOT, the following are true: the sensitivity of ς in twilight sky is weak.In cases with a small AOT, one possible explanation is that polarization effects cause non-negligible differences between the twilight sky brightness and color found by the forward model simulations and those found by measurement.The forward model uses a scalar approximation to model the radiative transfer.Saito et al. (2013) investigated the contribution of multiply scattered light to twilight sky near the horizon with the JACOSPAR RTM, from which the present forward model is descended.The fraction of multiply scattered light is in the range of 20-50 % near the horizon for a VZA of 90-85 • under small AOT conditions (0.05 < τ tr < 0.13).Mishchenko et al. (2006) demonstrated that scalar approximation leads to biases for calculated radiance, with the magnitude of the bias in the range of ± 3-5 % depending on the RAA, relative to the vector RTM outputs for AOTs of 0.1-0.3.This may cause an error in the retrieved brightness and color in twilight near the horizon.Another possible reason is the seeming insensitivity of ς in twilight sky when looking near the horizon (with SZA of ∼ 90 • ), and this sensitivity may be too low to overcome the uncertainties resulting from model errors and camera noise.An investigation of this possibility and improvements to the algorithm to compensate are left for future work.
The retrieval method proposed by this paper can be applied to aerosol measurements in desert, snow-covered, and polar regions, areas where it is difficult to apply the con-ventional solar reflection-based methods.Particularly in the polar region in winter, aerosol measurements are very limited, and the present twilight photometry method should be helpful to increase the availability of aerosol measurements.Furthermore, it is easy to extend the aerosol observations at many points with this method because of the use of a digital camera, which is both inexpensive and easily movable.Digital cameras are widely used in daily life, and many amateur photographers have the opportunity to take pictures of beautiful sights, including twilight sky.In addition to their beauty, twilight photographs have the potential to be used as measurement data for aerosol retrievals.Therefore, the present aerosol retrieval method has unique characteristics, and it may complementarily be used in combination with other measurement techniques.

Figure 1 .
Figure 1.An example of a twilight sky photograph taken on 23 March 2012 in the morning.

Figure 3 .
Figure 3.A prior vertical profile of aerosol extinction coefficient.

Figure 4 .
Figure 4. Sensitivity test results for normalized brightness (G n ) and color ratios (B/G, R/G) of twilight sky: (a, b) sensitivities to tropospheric aerosol optical thickness; (c ,d) sensitivity to stratospheric aerosol optical thickness; (e, f) sensitivity to coarse-fine particle volume ratio.Panels on the right side (b, d, f) show the dependence of normalized brightness variations along with the viewing zenith angles.Panels on the left side (a, c, e) show the dependence of color ratio variations along with the viewing zenith angles on R/G-B/G cross sections.The SZA and RAA are 93 and 0 • , respectively.

Figure 5 .
Figure 5. Sensitivity test results for color ratio values of the twilight glow (VZA of 85 • and SZA of 93 • at the solar azimuth angle).(a, b) show the he dependence of color ratio variations on the R/G-B/G cross section; (c, d) show the dependence of normalized brightness variations on the R/G-G n cross section.Dashed curves in the panels on the right (b, d) are isolines of coarse-fine particle volume ratio, and solid curves are isolines of tropospheric aerosol optical thickness.Dashed curves in the panels on the left (a, c) are isolines of stratospheric aerosol optical thickness, and solid curves are isolines of tropospheric aerosol optical thickness.

Figure 6 .
Figure 6.Sensitivity test results for color ratio values of the twilight sky near the horizon (VZA of 88 • and SZA of 91 • at the solar azimuth angle): (a, b) show the dependence of color ratio variations on the R/G-B/G cross section; (c, d) show the dependence of normalized brightness variations on the R/G-G n cross section.Dashed curves in the panels on the right (b, d) are isolines of tropospheric aerosol optical thickness, and solid curves are isolines of coarse-fine particle volume ratio.Dashed curves in the panels on the left (a, c) are isolines of coarse-fine particle volume ratio, and solid curves are isolines of stratospheric aerosol optical thickness.

Figure 7 .
Figure 7. Forward model errors (absolute values) in (a, b) R/G, (c, d) B/G, and (e, f) G n caused by the model assumptions as functions of SZA and VZA (see text for details).Panels on the left (a, c, e) have RAA =0 • , and those on the right (b, d, f) have RAA = 30 • .

Figure 8 .
Figure 8. Forward model errors (absolute values) in (a, b) R/G, (c, d) B/G, and (e, f) G n caused by uncertainties in the model parameters as functions of SZA and VZA.Panels on the left (a, c, e) have RAA = 0 • , and those on the right have RAA = 30 • .

Figure 9 .
Figure 9.Initial and retrieved (a) tropospheric aerosol optical thickness, (b) stratospheric aerosol optical thickness, and (c) coarse-fine particle volume ratio ς in the retrieval-error analysis.Black and gray marks (circles) show ς retrievals with ς > 1 (ς ≤ 1) under very and moderately clean conditions, respectively.

Figure 10 .
Figure 10.Comparison of aerosol property retrievals from camera measurements and sky radiometer measurements.Error bars show uncertainty of retrievals.Color scales show total retrieval cost J (x).The vertical axis and horizontal axis show retrievals from a camera and those from other instruments, respectively: (a) aerosol optical thickness, from camera and sky radiometer, with correlation coefficient r at 0.94 and 18 trials; (b) stratospheric aerosol optical thickness from camera and OSIRIS on board the Odin satellite; (c) coarse-fine particle volume ratio from camera and sky radiometer.

Figure 11 .
Figure 11.Photographs of twilight sky on the solar side, for the case studies of (a) 31 January, (b) 23 May, and (c) 3 June 2013.In all cases, SZA is uniformly 94 • , and VZA is 90.45 • at the horizon, 91.5 • at the bottoms of the images, and 30 • at the tops of the photographs.
For 31 January 2013 (case a), a very clean atmosphere with an AOT of 0.02-0.03 is simulated by the SPRINTARS and weak backscatter signals are obtained from lidar measurement.Moreover, spherical particles are dominant in this case.In contrast, the SPRINTARS products suggest that sulphate and dust particles were transported from the eastern part of continental China to Japan from 20 May to 26 May 2013.Consistently, the lidar measurements show the presence of nonspherical particles at altitudes of 3-6 km and the presence of spherical/nonspherical mixed particles in the boundary layer for 23 May 2013 (case b).For 3 June 2013 (case c),

Figure 12 .
Figure 12.Vertical profiles of aerosol extinction coefficients for altitudes from 0 to 6 km, as obtained from lidar measurements for the case studies.Data for 0 to 150 m are not available.

Table 2 .
Viewing and solar directions used for aerosol retrievals.

Table 3 .
Parameters for the aerosol vertical profile assumed in the forward model.

Table 4 .
Prior information and its error range for the state vector elements.

Table 5 .
Reference values and variabilities assumed for error analysis.

Table 6 .
Retrieval results and SPRINTARS model simulation results for case studies.
NA: not available