Interactive comment on “ Development of a new data-processing method for SKYNET sky radiometer observations ” by M . Hashimoto

This paper discusses changes in data processing that improve the retrieval of the single scattering albedo of aerosols measured remotely with the SKYNET sky radiometer. Improvements are compared to AERONET retrievals as a standard because AERONET uses "rigorous calibration routines". The authors demonstrate that the unknown surface albedo, incorrect field of view assignment, and cirrus cloud contamination could explain some of the bias. This led to improved quality control methods that minimize these error sources. They demonstrate that the steps taken to eliminate cirrus contamination may affect retrievals of coarse particles such as dust and steps need to be taken to recognize this situation a priori.


Introduction
Studies in recent years have pointed out that anthropogenic aerosols significantly offset the effect of global warming caused by anthropogenic greenhouse gases.Direct and indirect radiative forcings as large as −0.5 W m −2 and −0.7 W m −2 , respectively, have been estimated (IPCC, 2007), but such estimates still have a large uncertainty because aerosols in the atmosphere have large variations in time and space, and also their physical and chemical properties are complex.Many observational studies of aerosols have been conducted to reduce the uncertainties.Among these, ground-based measurement networks of scanning sunsky multi-wavelength photometers, e.g., NASA AERONET (Holben et al., 1998), can measure and retrieve aerosol optical properties that afford a key to accurately evaluating Direct Aerosol Radiative Forcing (DARF).In particular, the reliable retrieval of the single scattering albedo (SSA) of aerosols (e.g., Dubovik et al., 2002) is a unique function of these networks.In this regard, it is important to recognize that accurate retrieval of the SSA is more difficult than estimation of the value of aerosol optical thickness (AOT) and size distribution (e.g., Loeb and Su, 2010;McComiskey et al., 2008).We refer to the normal optical thickness as optical thickness for normal incidence (Chandrasekhar, 1960).Continued effort toward improving and validating SSA retrieval from networks is, therefore, an important task for us to increase the accuracy of modeling aerosol climate effects.
SKYNET (Nakajima et al., 2007), the focus of our study, is a ground-based network of scanning sun-sky photometers called sky radiometers (Prede Co., Ltd, Tokyo, Japan) with observation sites spread over Asia and other areas.Direct and diffuse solar irradiances are measured with the sky radiometers, and are analyzed to derive the aerosol optical properties, such as AOT, SSA, complex refractive index, and volume size distribution function (SDF) by using SKYRAD.pack,a software package to analyze the sky radiometer data (Nakajima et al., 1996).Aerosol optical properties retrieved from SKYNET have been used to investigate regional and seasonal characteristics of aerosols for climate and environmental studies and to validate satellite remote sensing results (e.g., Higurashi and Nakajima, 2002;Kim et al., 2005;Sohn et al., 2007;Pandithurai et al., 2009;Campanelli et al., 2010;Khatri et al., 2010;Takenaka et al., 2011).
There are several reports that the AOT is obtained with high accuracy compared with that of the standard Langley method and with those from AERONET (Campanelli et al., 2007;Che et al., 2008).However, it is pointed out that the SSA values from SKYNET are 3 % to 8 % larger than those of AERONET at the Beijing site, although the retrieved AOT values from SKYNET agree with those from AERONET within about 1 % (Che et al., 2008).There are several aerosol types with SSA value close to 1.It is, however, claimed that the SKYNET SSA sometimes becomes unnaturally close to unity, as found at the Phimai site in Thailand (H.Tsuruta, personal communication, 2011) and at the Hyderabad site in India (Badarinath et al., 2011).
Figure 1 compares the values of AOT and SSA for SKYNET, at a wavelength of 0.5 µm, and for AERONET, at a wavelength of 0.5 µm, which are interpolated by a polynomial fit in logarithmic space with wavelength, at Pune from April through December 2008.The figure shows that the SKYNET AOT values are in close agreement with those of AERONET, whereas the SSA values from SKYNET varied widely, with a tendency to become larger than those from AERONET.Additionally, some SSA data are unnaturally close to unity, as reported in previous studies.It is difficult to explain this difference by the difference in wavelengths.
It should be noted that the working instruments located at the observation sites, operational systems, and analysis algorithms are somewhat different between two networks.AERONET Cimel sun/sky radiometer scans the sky at both sides of the sun to check the symmetry of spectral radiances relative to the sun in the almucantar by the AERONET processing criteria.The criteria are based on the assumption that the aerosols are uniformly distributed over the ground, and is helpful for the screening of cloud contamination data.On the other hand, the SKYNET Prede sky radiometer scans just one side in all the azimuths.For the cloud screening of AOT, AERONET and SKYNET adopt Smirnov et al. (2000) and Khatri and Takamura (2009), respectively.We describe the SKYNET cloud screening in Sect. 3.2.Calibrations are also different between two networks.Each AERONET instrument is checked by means of intercalibration with reference instrument every 12 months at NASA Goddard Space Center, Maryland.The AERONET reference instruments are calibrated at Mauna Loa site in Hawaii for the calibration constants; V 0 , for direct solar irradiance measurement, by using the normal Langley plot method, and the lamp method for the determination of the instrument solid view angle for sky radiance measurement.For the SKYNET sky radiometer, calibrations for direct solar irradiance and sky radiance are carried out at each site by use of an improved version of the Langley plot method, and the solar disk scanning method for obtaining calibration constants and solid view angles, respectively.Furthermore, the inversion calculation methods adopted by these two networks are different.SKYNET adopts Phillips-Twomey method (Nakajima et al., 1996), while AERONET adopts the maximum likelihood method (MLM) combined with Phillips-Twomey method as a smoothing constraint (Dubovik and King, 2000).The AERONET algorithm includes the particle non-sphericity, but the SKYNET one does not in the present version.For the inversion method, we refer to it in Sect.2.3.
To develop the technique to derive more accurate aerosol properties in the atmosphere, it is important to estimate the error included in the retrieval values and processes.The purpose of this study is to perform sensitivity studies of various aspects of the present SKYNET algorithms, though the true validation of the algorithm should be done through comparison with in situ measurements, which should be our next task.In this study, therefore, we investigated possible causes of error in the SSA retrieved from SKYNET through numerical experiments using a radiation transfer code and through data analysis with real observation data.On the basis of the results of the investigation, we suggest an improved method to estimate SSA more accurately.In Sect.4, we evaluate DARF using improved aerosol optical parameters from the Pune and Beijing sites.

Numerical tests
This section describes sensitivity tests to investigate possible causes of error in SSA values from SKYNET.

SKYRAD.
pack is an open source software package released on the OpenCLASTR web page (http://www.ccsr.u-tokyo.ac.jp/ ∼ clastr/) for data analysis of direct and diffuse solar radiation measurements to give the AOT, SSA, complex refractive index, volume size distribution function (SDF), phase function, and so on.The present version of SKYRAD.pack(version 4) is used for operational analysis by the SKYNET data center at Chiba University (http://atmos.cr.chiba-u.ac.jp/) and the European Skynet Radiometers network (ESR, http://www.euroskyrad.net/index.html).
where F 0 is extraterrestrial solar irradiance in W m −2 µm −1 , which is determined by the improved Langley method; E( ), the monochromatic sky irradiance in W m −2 µm −1 measured at a scattering angle ; τ , ω, and P ( ) are the total optical thickness, SSA, and scattering phase function at scattering angle for the column atmospheric air mass, respectively; m 0 is the optical air mass; , the solid view angle (SVA) of the sky radiometer; and q( ), the contribution of multiple scattering.
In the actual analysis, the measured values, F and E( ), are voltage outputs, and F 0 is a parameter to be determined by the improved Langley method in the voltage unit, called a calibration constant.
is determined by the solar disk scanning method.
Equations ( 1) and ( 2) include information on the aerosol microphysical parameters along with molecular scattering, such as vertically integrated aerosol SDF, v(r), as a function of particle radius, r, and complex refractive index, m, as a function of wavelength, λ.The AOT, τ a , and effective SSA, ω a , of aerosols in the atmospheric column can be derived from the following relations: where K e and K are kernel functions of the Fredholm integral equations as a function of size parameter x(= 2π/λ and the complex refractive index.To solve Eqs.(1) through (3) as an inversion problem, these equations are formulated into a matrix formula by discretization of λ, , and ln(r), and by weighting to make the solution stabilized as proposed by Nakajima et al. (1983): where f is an observation vector; K = K( m(λ)), a matrix of kernel coefficients calculated for fixed values of complex refractive index m(λ); and x, a state vector containing values of size distribution v i = v (r i ) with r i equidistant on a logarithmic scale, i.e., ln (r i+1 ) − ln (r i ) = const.SKYRAD.packuses two methods for removing the multiple scattering contribution in Eq. ( 2) and inversion of Eq. ( 4).The present inversion program, called version 4, uses the iterative relaxation method of Nakajima et al. (1983Nakajima et al. ( , 1996) ) and a statistical regularization method (Turchin and Nozik, 1969) to derive an optimal solution by minimizing the following cost function as proposed by Phillips (1962) and Twomey (1963): where B is the second order derivative matrix with respect to the particle size in ln (r), to generate a priori information that force the obtained solution x to be a smooth function of ln (r).The constant γ is a Lagrange multiplier coefficient and is chosen so as to minimize the first term of the right-hand side of Eq. ( 5).The solution of Eq. ( 4) provides smooth retrieval of size distribution v (r) corresponding to the minimum of e 2 defined by Eq. ( 5).However, in such an approach, both the solution v (r) and e 2 depend on the assumed value of the complex refractive index m(λ), i.e., e 2 = e 2 ( m) Correspondingly, if the value of m(λ) is unknown, the minimization of Eq. ( 5) can be performed for a set of different values mk (λ) (k = 1, 2, . . ., N k ), and mk (λ) corresponding to the smallest e 2 ( m) can be considered as a retrieved value of the complex refractive index.The disadvantage of such a retrieval approach is that the retrieved mk (λ) can only be chosen from the predefined set of values mk (λ) (k = 1, 2, ..., N k ).
If the values of complex refractive index m(λ j ) are directly included in the state vector x, the discretized system analogous to Eq. ( 4) becomes non-linear and should be solved by non-linear iterative techniques.In addition, in order to ensure uniqueness of solution, the minimized cost function of Eq. ( 5) should be modified so that it includes constraints on the retrieved complex refractive index.This can be achieved, for example, by using the maximum likelihood method (MLM) as defined by Rodgers (2000).This method is employed by version 5 and will be described in more detail in the next section.

Sensitivity tests for parameters for data processing
To investigate possible causes of error in the SSA retrieved from SKYNET, we analyzed simulated direct solar irradiances and sky radiances using SKYRAD.pack.
Candidates for possible causes of error are as follows: (1) errors in the input geophysical parameters, such as ground surface albedo and initial condition values of complex refractive index; (2) instrumentation errors, such as minimum observable scattering angle and stray light, and calibration constants and solid view angles of the radiometer; (3) errors caused by inversion algorithms (type of inversion algorithm); and (4) errors attributed to the condition of the atmosphere, such as homogeneity in space and time and the total amount of aerosols in the atmosphere, that is, the value of AOT.The concept and methodology of the tests are somewhat consistent with earlier sensitivity studies conducted by Dubovik et al. (2000) for evaluating the accuracy of AERONET retrievals.
We simulated observation data by using the Open-CLASTR software package for radiative transfer code called Rstar-6b (Radiance System for Transfer of Atmospheric Radiation version 6b) using the formulae of Nakajima andTanaka (1986, 1988).We used the AFGL US standard for atmospheric conditions and the rural aerosol model incorporated in Rstar-6b, which is based on the rural model of the WCP report (Deepak and Gerber, 1983).We set the solar zenith angle θ 0 = 60 • and AOT = 0.5 for a wavelength of 0.5 µm.The ground surface albedo A g , calibration constant F 0 , and SVA are given as 0.2, 1, and 2.5 × 10 −4 str, respectively, at wavelengths of 0.4, 0.5, 0.675, 0.84, and 1.02 µm.In the analysis, we assumed an error of ±5 % for F 0 , ±5 % for SVA, ±50 % (±0.1) for A g , and ±0.005 for the initial value of the imaginary part of the refractive index.We also conducted sensitivity tests in which we changed the minimum observable scattering angle of the sky radiometer from 3 • to 2 • , and in which we also increased the diffused intensity at scattering angle 3 • by 5 % at each wavelength.We compared retrieved SSA values with and without the assumed errors to seek possible causes of error in the SSA value that are consistent with the observed errors.
We show the results of the sensitivity tests in Fig. 2. In the cases where we changed the initial value of the imaginary part of the refractive index, added stray light, and changed the minimum observable scattering angle, there are no significant differences in the SSA from the case without errors, with the differences being less than 1 %.On the other hand, in the case where we assumed an error for A g , SVA, and F 0 , the result shows differences in the SSA relative to the case without errors.
Figure 3 shows the difference in SSA at a wavelength of 0.5 µm between cases with and without error, defined as {[SSA (with error) − SSA (no error)]/[SSA (no error)]}.When we assumed errors of −50 % (−0.1), −5 %, and −5 % for A g , SVA, and F 0 , respectively, the differences in SSA were +3.7 %, +3.0 %, and +5.5 % at a wavelength of 0.5 µm.For the AOT at a wavelength of 0.5 µm, there was no difference, defined as {[AOT (error) − AOT (no error)]/[AOT (no errors)]}, when we introduced an error in A g and SVA, but there was a +2.8 % difference when we introduced an error in F 0 .
On the basis of sensitivity tests, it is concluded that an error in the calibration constant (F 0 ) causes an error in both retrieved SSA and AOT.However, according to a reported comparison of calibration constants from SKYNET with those from AERONET, the improved Langley method adopted by SKYNET seems to yield accurate calibration constants (Campanelli et al., 2004).Furthermore, Che et al. (2008) show that retrieved SSAs are different between SKYNET and AERONET by 4 % at a wavelength of 0.67 µm, whereas the retrieved AOTs from both networks differ only by 1.3 % at a wavelength of 0.67 µm.Therefore, it is unlikely that errors in the calibration constants (F 0 ) are a cause of the SSA difference between the two networks.On the other hand, it is likely that the SSA differences are caused by errors in A g and/or in SVA, because these errors do not cause an error in AOT in the retrieval process.The simulation result indicates that an underestimation of A g or SVA leads to an overestimation of SSA and also of the amount of aerosols, i.e., the value of AOT in the atmosphere is one of shows the result without error or change applied to the parameters.The following errors and changes were applied: ±5 % to F 0 , ±5 % to SVA, ±50 % to A g , ±0.005 to the initial value of the imaginary part of the refractive index, an increase in the diffused intensity at scattering angle 3 • by 5 % at each wavelength (stray light), and a decrease in the minimum observable scattering angle of the sky radiometer from 3 Fig. 3. Summary of the errors in retrieved SSA in percentage at a wavelength of 0.5 µm in the sensitivity experiments.The horizontal axis shows parameters that are applied errors or changes, and also shows the error from the true value in retrieved AOT ( AOT/AOT).The longitudinal axis is the error from the true value in retrieved SSA.SVA is the solid view angle; A g , the ground surface albedo; F 0 , the calibration constant; m i , imaginary part of refractive index; Obs angle min , the minimum observable angle.Applied error or changes are in parentheses.
the causes of error in the SSA.Although the value of A g depends on wavelength and ground conditions, the value used in data processing at the SKYNET data center is set to 0.1 for each wavelength.As for the SVA value, it is determined by the disk scan method, called the Sun scanning method in Nakajima et al. (1996).At present, the disk scan in SKYNET varies according to the observation site.For some sites, the observation software is set to perform a disk scan periodically at intervals of a few days (e.g., 1 week or 10 days) at a certain time (e.g., 11:00 a.m.LT).However, for other sites disk scan data are missing for long periods (e.g., more than 1 yr).The stability of the estimated SVA time series indicates that possible errors included in the SVA are within 5-6 % because of the lens degradation and color aberration of the single lens used in the radiometer, which produces an error of ∼ 3 % in the SSA.SKYNET SVA in the standard analysis is calculated by the point source method (Nakajima et al., 1996) by using disk scan data measured at the interval of few days (e.g., 1 week or 10 days) at certain time (e.g, 11:00 a.m. in the morning), and is one month averaged data of SVA of which deviation is within 10 %.However, for some sites, disk scan data are missing for a long time (more than 1 yr).We took standard deviation of the error of determined SVAs at Pune and Beijing, and the standard deviation is within 5-6 %.We therefore put 5 % error in SVA for the test.For the error 3 % in SSA due to SVA error, it was found from the result of the sensitivity test that 5 % error in SVA occurs about 3 % error in SSA.
We also investigated the error in SSA that is due to the total amount of aerosols in the atmosphere, that is, AOT.We used the data simulated by Rstar-6b for this experiment, and set the simulation conditions as follows: the atmospheric condition is the US standard atmosphere, the ground surface albedo A g , calibration constant F 0 , and SVA are given as 0.2, 1, and 2.5 × 10 −4 str, respectively, at a wavelength of 0.5 µm.The solar zenith angle is given at θ 0 = 60 • .We examined the error in retrieved SSA when AOT is 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, and 1.0, for five aerosol models such as dust-like and rural aerosol types in the WCP report (Deepak and Gerber, 1983), and water-soluble, waterinsoluble, and mineral accumulate aerosol types (Hess et al., 1998).As shown in Fig. 4, the averaged error between the retrieval value and true value of the SSA, which is defined as |SSA(retrieval) − SSA(true)| /SSA(true) , is about 1 % when the value of AOT is larger than 0.2 or 0.3, and the error of the SSA and its variation (standard deviation in the figure) became larger with decreasing AOT values.Therefore, it could be said that the condition of low AOT affects the retrieval accuracy of SSA, especially when AOT is less than 0.2.It is reported that the error in retrieved SSA decreases with increasing aerosol optical thickness for AERONET (Dubovik et al., 2000).We found a similar tendency for SKYNET sky radiometer data analysis in this paper.The amount of aerosols, AOT, in the atmosphere easily varies in one day, and in many places the AOT values are lower than 0.2 or 0.3 at a wavelength of 0.5 µm almost every day.Therefore, in the low AOT case, we should note that it is difficult to retrieve an accurate SSA by using the present algorithm.
In conclusion, it is possible that the errors associated with A g and SVA, and the amount of aerosols in the atmosphere, are the causes of error in the SSA.However, the errors in evaluating SSA due to these parameters and the low AOT condition will be less than 5 %, and these errors should cause both underestimation and overestimation of SSA, so it is difficult to explain the reported overestimation for all these sites.

Sensitivity tests for difference in inversion algorithms
SKYRAD.pack version 4 utilizes the Phillips-Twomey method to minimize the cost function given in Eq. ( 5) to determine aerosol parameters, such as the volume size distribution and refractive index, as described in Section 2.1.In this method, the a priori information for stabilizing the illconditioned Fredholm integral equation is the smoothness of the retrieved SDF given by the B 2 matrix.On the other hand, similar a priori constraints can be applied in the framework of a statistical estimation approach, for example, as implemented in AERONET retrieval via a non-linear maximum likelihood method defined in the logarithmic space of the retrieved variables SDF and m.Here, we base our approach on the MLM method as defined by Rodgers (2000).This method is based on the Bayesian theory.
where p is the probability density function and is defined as the Gaussian distribution; and x and f denote state and measurement vectors, respectively.In the MLM method, x is chosen so that the posterior probability p (x|f ) becomes maximum under the condition that a priori information is already given.Organizing this non-linear equation such that p (x|f ) = max, we obtain the following equation in the tangential space to be solved by a Newtonian method: where x k is the solution at the k-th iteration step; f k = f (x k ), an observation modeled using x k ; x a , the a priori value of x; S e , the measurement error covariance matrix; S a , the covariance matrix defined by a priori and state values, S a = (x − x a ) (x − x a ) T ; and U, the Jacobi matrix, ∂f /∂x.Further, it should be noted that in a manner similar to the AERONET approach (Dubovik and King, 2000), we used a logarithmic scale for the volume size distribution and complex refractive index to prevent x from having a negative value.Thus, analogous to the AERONET retrieval approach, the retrieval algorithm used in version 5 allows rigorous retrieval of both the aerosol size distribution and the spectral complex refractive index.At the same time, some differences remain between the version 5 and AERONET retrieval methods.Specifically, AERONET retrieval does not use the MLM of Rodgers (2000) corresponding to Eq. ( 7).Instead, it uses the multi-term LSM (Least Square Method) described in papers by Dubovik and King (2000), Dubovik (2004), and Dubovik et al. (2011).Similarly to Rodgers' method, multi-term LSM relies on a statistical estimation approach; however, it differs by allowing simultaneous use of multiple a priori constraints.For example, it can include both a priori constraints on the smoothness of retrieved functions and a priori estimates of the state vector.Specifically, instead of Eq. ( 7), the solution of the AERONET algorithm is expressed as follows: where is the smoothness matrix.This is used by Dubovik and King (2000) to apply different smoothness constraints on the size distribution (v (r)), on the spectrally dependent real part of the complex refractive index (n (λ)), and on the spectrally dependent imaginary part of the complex refractive index (m (λ)).The state vector retrieved by the AERONET algorithm can be denoted by x T = (x v ; x n ; x m ; x sph ) T , where x v is a vector including N v = 22 values of v (r i ); x n , a vector including N λ values of n λ ϕ ; x m , a vector including N λ values of m λ ϕ ; and x sph , the value describing the fraction of spherical particles.Using this notation, the smoothness matrix in Eq. ( 8) can be denoted as where and N λ × N λ accordingly, defined via the matrices B v , B n , B m of coefficients for estimating the corresponding derivatives.Specifically, Dubovik and King (2000) constrain the third derivatives of v (r), the first derivatives of n (λ), and the second derivatives of m (λ).It should be noted that matrix B v is analogous to the matrix B used in Eq. ( 5) and defined as prescribed by Phillips (1962) and Twomey (1963); however, matrices B n and B k are slightly more complex because they are defined for non-equidistant discretization points λ ϕ (λ ϕ+1 −λ ϕ = const).The explicit definition of such matrices, as well as the definition of the corresponding Lagrange parameters γ in Eq. ( 8), is given by Dubovik and King (2000) and explained in detail by Dubovik et al. (2011).Thus, the solution given by Eq. ( 8) minimizes the following cost function: It is important to note that the AERONET algorithm uses a priori estimates x a only for two retrieved parameters {x} 1 and {x} 22 corresponding to the values of the retrieved size distribution of the smallest size class (r = 0.05 µm) and the largest size class (r = 15 µm), i.e., the matrix S a has all zero elements except for the first and 22nd elements of the diagonal, i.e., {S} 1 1 and {S} 22 22 .This constraint was introduced by Dubovik et al. (2006) to avoid unrealistically increasing tails of size distribution appearing due to the very low sensitivity of sky radiometer observations to very small and very large particles.This constraint has a rather "cosmetic" objective because, practically, it does not affect the minimum value of the cost function and, as was shown by Dubovik et al. (2000), the retrieval errors for size distribution tails are very high.At the same time, constraining the size distribution tails to small values may give the wrong impression of an absence of very large and very small particles.The correct interpretation should state that there are very large and very small particles that make a contribution sufficient to change the values of the cost function (while the volume of those particles, in principle, cannot be negligible compared to the volume of the rest of the particles).
The values of {x a } 22 are given by AOT(0.44)×0.002.The initial guess for the size distribution is chosen as a straight line with values AOT(0.44)× 0.01; for the refractive index, the initial guess is spectrally independent with values n = 1.45 and m = 0.005.
In order to study whether the reported SSA differences can originate from the difference between the inversion algorithms of SKYNET and AERONET, we performed various test simulations with SKYRAD.packversion 4 and a new version 5, which has been developed using the algorithm that is explained above and is closer to the AERONET algorithm.Version 5 uses an a priori SDF of a bimodal log-normal function, with r m1 = 0.1 µm, r m2 = 2.0 µm, S 1 = 0.4, S 2 = 0.8, C 1 = 1.0 × 10 −12 , and C 2 = 1.0 × 10 −12 following reported climate values (Higurashi et al., 2000).For a priori estimates of the real part (n) and the imaginary part (m) of the refractive index, we usually set n = 1.5 and m = 0.005, which are spectrally independent values.One of the key differences between version 4 and version 5 is that version 5 uses a priori estimation, but version 4 does not.
We first performed numerical tests by using Rstar-6b with two aerosol models incorporated in Rstar-6b, i.e., the dustlike and rural aerosol types in the WCP report (Deepak and Gerber, 1983), and three other aerosol models of Hess et al. (1998), i.e., water-soluble, water-insoluble, and mineral accumulate aerosol types, and we set AOT = 0.5 for each case.As shown in Figs. 5 and 6, the differences for the retrieved SSA are found to be less than 0.01, and there were not large differences for the SDF.Among various other tests, we found one case in which a noticeable difference exists between version 4 and version 5, as shown in Fig. 7.In this case, we assumed that a tested true SDF includes a large amount of coarse particles of the dust-like aerosol type with radius greater than 10 µm.The figure shows that version 4 could retrieve the SDF relatively well, including the coarse mode, in comparison with version 5, because the smoothness condition given by Eq. ( 5) allows the retrieved SDF to be distributed beyond a 10 µm radius.On the other hand, version 5 underestimated the coarse mode of the SDF because of the strong SDF constraint condition given by Eq. ( 11) with a small model radius r m2 = 2.0 µm for the coarse mode SDF.The value of the SSA is then underestimated to compensate the reduced light absorption by coarse particles.Although not shown in a figure, we found no significant error in the retrieved AOT because the inversion process did not bring about a large change in the retrieved direct solar radiation.As discussed later, an enhanced coarse mode SDF is possibly required for several dust storm cases (Chepil, 1957;Gillette et al., 1978;and Feng et al., 2007).The test indicates that version 4 can retrieve accurate SSA values in comparison to version 5 in this case.It is possible that version 5 may underestimate the SSA value because of the constraint on the SDF.
We investigated the effect of the large coarse part of the SDF.We calculated relative intensity at the surface  and Gerber, 1983) and for the water soluble, water insoluble, and mineral accumulate aerosol types of Hess et al. (1998).
./!012'(3*-' 7/85&'0%26$9$5' at a wavelength of 0.5 µm by using Rstar-6b, and compared it with and without a cut above 10 µm for the SDF.We used dust-like aerosol model incorporated in Rstar-6b, which is described by a mono-modal lognormal SDF with mode radius 6.0 µm.The AOT setting at 0.5 µm is 0.5.Figure 8 shows the difference between the relative intensity with and without a cut above 10 µm for the SDF { R = [R(cut above 10 µm) − R(no cut above 10 µm)]/R(no cut above 10 µm)}.The minimum scattering angle is 3 degrees.From this result, the lack of a large coarse part in the SDF causes overestimation of sky radiance at all observation angles.The intensity of forward scattering near 0 degree increases with increase in particle size, but the simulation shows that the diffused intensity without over 10 µm particles is larger than that with over 10 µm particles in the region of measured scattering angles (> 3 degree) in the same condition of AOT.Hence, a lack of large particles (> 10 µm radius) causes "overestimation" of radiance.It is likely that version 5 works to decrease the SSA value to dim the sky radiance in the calculation when a tight constraint on the SDF for particles with radius over 10 µm is applied.The constraint that is used in AERONET data processing for suppressing the concentration of particles corresponding to the largest class (r = 15 µm) may have a similar effect.However, there are other differences, e.g., a priori estimates, other constraints, or determination of a covariance value, etc.Therefore, we need further investigation of this finding.We further performed, as shown in Fig. 9, a test retrieval of a cirrus contamination case for the aerosol atmosphere by using Rstar-6b with a combination of two particle models without an enhanced coarse mode incorporated into Rstar-6b, i.e., the dust-like and ice particle types of the WCP report (Deepak and Gerber, 1983) represented aerosol and cirrus particles, respectively.We set AOT = 0.5 for dust-like particles and an optical thickness of 0.1 for ice particles at a wavelength of 0.5 µm.For the other conditions, we chose US-standard atmosphere and set the solar zenith angle at θ 0 = 60 • .The ground surface albedo A g , calibration constant F 0 , and solid view angle (SVA) are given as 0.2, 1, and 2.5×10 −4 str, respectively, at wavelengths of 0.4, 0.5, 0.675, 0.84, and 1.02 µm.In simulation data, we consider the forward scattering light that measured as a direct solar irradiance in the field of view, and SKYNET retrieval algorithm has the process to remove this overestimation of direct irradiance.
In this case, the SDF consists of the first mode due to aerosol particles and a second coarse mode due to cirrus particles, as shown in Fig. 9.
The inversion result shows that version 4 retrieved the SDF, including contaminating cirrus particles larger than 10 µm, but version 5 successfully filtered out the cirrus particles by the constraint of a reduced SDF for particles with radius greater than 10 µm.As a result, the SSA value retrieved by using version 5 became closer to the true value of SSA.This test indicates that cirrus contamination can cause a serious overestimation of SSA from SKYNET as reported, whereas SSA retrieved by using version 5 is robust and without significant error.3 Analysis of observation data

Case studies
To investigate whether enhanced coarse dust particles and/or cirrus contamination can affect the real SSA retrieval, we analyzed SKYNET data at the Pune site (18.616• N/73.800 • E) in India, which is one of the collocated sites of SKYNET and AERONET radiometers.We used the data for 23 October 2008, when cirrus was detected over the Pune site by the CALIPSO lidar, as illustrated in Fig. 10.As shown in Fig. 11, it is found that version 4 retrieved an SDF with an enhanced coarse mode that seems to be a cirrus particle contribution.
On the other hand, version 5 gave an SDF without a large coarse mode as an a priori value.Values of SSA from version 4 were larger than those from version 5 for the reason we proposed in the preceding section.This result is consistent with the result of numerical experiments on cirrus contaminations, so it is likely that, in this case, the SSA retrieved by version 4 was overestimated because of cirrus contamination.In the period of cirrus contamination, AERONET consistently eliminates the data through their symmetry check of the two sides of the almucantar scan (Holben et al., 2006), whereas the SKYNET scan only has one side of the almucantar and thus cannot use such symmetry quality check.
We next analyzed the SKYNET data at the Beijing site (39.586• N/116.229• E) in China on 14 April 2004 at 09:00 UTC. Figure 12 shows NIES lidar data located at 39.97 • N and 116.37 • E. From the figure, the backscattering intensity shows dense blue color and the depolarization ratio is about 0.1 to 0.2 to a height about 3 km, and it suggests nonspherical particles existing in the amosphere.Furthermore, the high ratio between attenuated backscattering coefficients (Int 1064 /Int 532 ) indicates large sized particles.Additionally, the Ångströme exponent from SKYNET was consistently as small as 0.49, indicating the existence of large particles, and coarse mode particles were detected by in situ measurements (Wu et al., 2009).We compared the SSA value and SDF retrieved by SKYRAD.packversion 4 with AERONET Level 2 data in Fig. 13.As shown in the figure, coarse particles with radius around 10 µm exist in the SKYNET results, whereas the SDF from AERONET does not include particles larger than 10 µm.The SSA values from AERONET are lower than  23October 2008, 08:36:00.2to 08:49:28.9UTC), for aerosol types of dust (yellow), polluted continental aerosol (red), polluted dust (brown), clean continental aerosol (green), and smoke (black).The data were taken from the NASA CALIPSO web site (http://www-calipso.larc.nasa.gov/products/lidar/browseimages/ production/).those from SKYNET, consistent with the result of numerical experiments in the previous section.Since cirrus cloud was not detected from the lidar, it is likely that the difference in the two products was caused by the difference in the inversion method.However, we cannot conclude this point, because there are several differences between the AERONET algorithm and SKYRAD.pack,e.g.input parameters such as surface albedo, instrument calibration method and scanning system.Therefore, more studies, including validation of the SDF, will be needed in the future.

Development of a quality control algorithm
In this section, we mention and suggest an approach to the issue of quality control of observation data and cloud screening.
The present standard process of quality control in SKYNET applies a retrieval error between observations and calculated theoretical values by using retrieval values, σ obs , where (τ meas λ i and R meas λ i ) and (τ λ i and R λ i ) are measured and retrieved observation vectors for the AOT and relative sky radiance defined by Eqs. ( 1)-( 3) at measurement wavelengths.Now we set W e = W P = 1/N total , where N i , N j , and N total = N i +N i ×N j indicate the number of measured wavelengths, scattering angles, and their total, respectively.In the present standard retrieval process in SKYNET, we remove the data if the value of σ obs is larger than 0.2.
The method of Khatri and Takamura (2009) is applied to cloud screening as standard procedure by SKYNET.In this process, the data identified as corresponding to clouds are removed.The first test separates cloud-affected and cloud-free periods of observation days by examining global irradiance data.If the standard deviation of the ratio of the observation flux density from irradiances to the theoretical flux density for each 5 min is equal to or larger than 0.02, that time period is kept in the cloud-affected data group; otherwise, it is kept in the clear sky data group.Secondly, the spectral dependency behaviors of aerosols and cloud, which use a spectral variability cloud screening algorithm (Kaufman et al., 2006) as a reference, are used (Eq.13).The 15-min data are used.
τ 870 nm − τ 400 nm (τ 870 nm /τ 400 nm ) > 0.0075 + 0.03τ 675 nm .(13) Here, τ is an observation value, τ is the maximum deviation of neighboring data, and the subscript notation indicates wavelength in nanometers.The numbers "0.0075" and "0.03" correspond to a noise error and the influence of refractive humidity, respectively.This criterion is applied to the data in the cloud-affected data group, and the data are regarded as cloud-affected data if the data fail to meet the criterion.Finally, statistical analysis tests in Eq. ( 14) are performed to remove outliers that pass the first and second checks, in an approach similar to Smirnov et al. (2000).
This process is applied to the triplet data in a minute.The present cloud screening relies heavily on the global flux test and needs global irradiance data but, almost uniformly, the observation sites in SKYNET do not conduct an observation of solar irradiance.Furthermore, cirrus contamination data are difficult to remove as cloud-affected data.
On the basis of the result of numerical experiments and real data analysis, we develop a quality control (QC) algorithm in this subsection to estimate more accurate SSA.
We drop data according to the following three conditions: (C1) the AOT is less than 0.4 at a wavelength of 0.5 µm, using the AERONET Level 2.0 QC algorithm, because the retrieval error in SSA rapidly increases with decreasing AOT (Dubovik et al., 2000).AERONET uses 0.44 µm as the wavelength for this algorithm, but there is no observation at a wavelength of 0.44 µm in SKYNET.Therefore, instead of AOT(0.44 µm), we use AOT(0.5 µm).(C2) We then reject data with a large deviation of the retrieved observation vector from the measured observation vector by using Eq. ( 12).We defined σ obs = 0.07 as a threshold for data rejection, through semi-empirical judgment of measurement errors found in the data analysis of large volume data sets in the past, including instrumental errors and errors in radiative transfer modeling and optical modeling of scatters.(C3) We then pose a condition regarding the magnitude of the coarse mode of the SDF given by Eq. (3c).
C v × v (2.4 µm) < max{v (7.7 µm) , v (11.3 µm) , v (16.5 µm)} , (15) where C v is a threshold coefficient to be determined for optimum rejection of cirrus contamination.This condition is set to warn the system that the retrieved SDF from version 4 includes a large volume of coarse mode particles larger than 10 µm.We set C v = 2 from the analysis of data at the Pune and Beijing sites, which enables rejection of most cirrus contamination cases.There were some dust cases that also have a large coarse mode, but the magnitude of the coarse mode is lower than that of cirrus contamination and can pass through this condition.The value of C v is determined from one year of data of the Pune and Beijing sites.It will be necessary for future work to determine this value after we collect more dust day data and cirrus contamination data.
We applied this new screening algorithm to the data from version 4 at the Pune and Beijing sites.The data used are from April to December 2008 at the Pune site and from February to September 2004 at the Beijing site.In Fig. 14, we show the monthly mean SSA values and standard deviation with error bars, before and after data screening at a wavelength of 0.5 µm by applying the conditions of C1 through C3 and each condition individually to Pune and Beijing.The screening conditions C1 to C3 contribute to the data screening almost equally, but the results for each condition were different from month to month and also depended on the site.The variability of the SSA was reduced by C1 at Beijing, but the results at Pune did not show a similar reduction.By applying condition C3, the monthly mean SSA at Pune became lower than before screening, because the SSA values that were very high or close to unity were removed by condition C3.On the other hand, the results at Beijing did not show a similar reduction in high SSA values.More detailed investigation of the results shows that the contribution of C2 is larger than those of C1 and C3 at both sites.This result indicates that the σ obs value is useful for detecting illconditioned data caused by cirrus contaminations, horizontally and/or temporally inhomogeneous aerosol stratification, and so on.
Figure 15 shows the normalized frequency distribution of SSA values retrieved by the system with and without data screening by all three conditions, and also shows its value from AERONET.It is found from the figure that SSA variability is significantly decreased after data screening.Moreover, the frequency of cases of SSA values close to unity was noticeably reduced.In addition, the SSA values after data screening become closer to the values from AERONET than they were before the data screening.

Discussion
From the analyses in the preceding sections, the retrieved SSA values from SKYNET at the Pune and Beijing sites became closer in agreement with those from AERONET after the data screening.We used Level 2 AERONET retrievals with that cloud screening.Using the results shown in Fig. 15, we calculated the difference in the SSA values between SKYNET and AERONET SSA = N i=1 SSA SKYNET,i − SSA AERONET,i 2 /N at a wavelength of 0.5 µm.For the calculation, we selected only data for which the measurement times of both networks are within 1.2 h of each other.We also linearly interpolated the SSA values from AERONET at a wavelength channel of 0.5 µm by using SSA values at 0.44 µm and 0.67 µm.The differences in the spring SSA at Pune in May and Beijing in April were 0.073 and 0.008, respectively, and the differences in the autumn SSA at Pune in October and Beijing in September were 0.017 and 0.043, respectively.In order to investigate why the SSA difference in the spring at Pune is as large as 0.073, we made a case study of a dust storm phenomenon on 11 May 2008.Such dust storms are frequently observed in the spring at Pune.As shown in Fig. 16, a backward trajectory analysis revealed that the dust particles were transported from the Arabian Peninsula.Figure 17 shows that the SKYNET SDF on 11 May 2008, had more large coarse particles over 10 µm than that of AERONET, similar to the case of Beijing shown in Fig. 13.From the figure, this difference in SDFs causes a difference in the SSA in the dust season at Pune.It should be stressed that such dust cases seldom arise in observations and that most of the SSA corrections were needed to address cirrus contamination that passed through an insufficient SKYNET cloud screening process.
On the basis of the aerosol parameters obtained by using the new data screening algorithm, we calculated monthly mean of Direct Aerosol Radiative Forcing (DARF), as shown in Fig. 18, by using the approximate formulae of Nakajima et al. (2007) for DARF at the top of atmosphere (TOA) and the bottom of the atmosphere (BOA), as a simple test of how the screened data can approach the AERONET results.The figure shows that the DARF values approach those from AERONET, indicating the new data screening algorithm is effective for SKYNET to improve their SSA retrievals to attain enough accuracy for aerosol forcing estimation.In the case of a large difference in SSA ( SSA = 0.073) in May at the Pune site, the difference in DARF between the two networks was 5 W m −2 at BOA and 10 W m −2 at TOA, even though the difference in AOT between the two networks was about 0.01 at a wavelength of 0.7 µm.We need, therefore, more work in the future to identify the cause of the large SSA difference in the dust storm case to enable better DARF estimation.

Conclusions
We found five sources of possible errors to explain SSA overestimation within the estimation process itself or when comparing values with those of AERONET.These sources are (1) an underestimation of SVA; (2) an underestimation of ground surface albedo; (3) the amount of aerosols in the atmosphere, i.e., low AOT condition; (4) cirrus contamination; and (5) the effect due to dust particles larger than 10 µm.In the low AOT condition, the result corresponds to Dubovik et al. (2000).For the first two sources, it was found that errors of ±5 % in the SVA or ±0.1 in the ground surface albedo result in errors of about ±3.0 % in the SSA.However, reported uncertainties in the SSA retrieval are about ±0.05 (present study; Loeb and Su, 2010) and, furthermore, such error sources can produce not only overestimation but also underestimation in SSA, making it difficult to explain why the SSA is consistently overestimated.We then found that the cirrus cloud contamination cases can be screened by three conditions (C1, C2, C3) given in Sect.3.2.This screening algorithm brought SKYNET SSA and DARF values into close agreement with those of AERONET, within less than 5 W m −2 and 10 W m −2 at TOA and BOA, respectively.One exception occurred for some data at the Pune and Beijing sites, when coarse mode dust particles prevailed at the observation sites.It is found that the version 5 of SKYRAD.pack with posing a constraint of suppressed coarse mode particles of radius larger than 10 µm causes the overestimation of SSA in the dust case with significant concentration of over 10 µm particles.We need more studies to conclude this is the cause of the SSA difference between AERONET and SKYNET in the case of dust case, because there are many other differences between two networks.We reserve more careful investigation of this exception for future work.Nonetheless, it will be beneficial for the two networks to establish suitable a priori SDFs for the cases of enhanced coarse particles larger than 10 µm.There are past reports (e.g.Zhang et al., 1998;Mikami et al., 2006;Feng et al., 2007;and Formenti et al., 2011) that show measured SDFs of soil particles with an extended tail for sizes larger than 10 µm, though these cases are not quite so common in dust cases (e.g.Reid et al., 2003Reid et al., , 2008;;and Johnson and Osborne, 2011).We need to accumulate a priori SDF information for dust cases.
Fig. 1.Time series of AOT and SSA at a wavelength of 0.5 µm at the Pune site.The upper and lower panels show AOT and SSA, respectively.The period is from April through December 2008.The horizontal axis is the time of Julian day from 1 January 2008.Daggers and crosses indicate the results of SKYNET retrieved by SKYRAD.packversion 4 and the results of AERONET, respectively.

Fig. 2 .
Fig.2.SSA values in the sensitivity experiments as a function of wavelength at 0.400, 0.500, 0.675, 0.870, and 1.02 µm.and indicate the results with an error or with a change to the parameter applied as described in the graph legends.shows the result without error or change applied to the parameters.The following errors and changes were applied: ±5 % to F 0 , ±5 % to SVA, ±50 % to A g , ±0.005 to the initial value of the imaginary part of the refractive index, an increase in the diffused intensity at scattering angle 3 • by 5 % at each wavelength (stray light), and a decrease in the minimum observable scattering angle of the sky radiometer from 3 • to 2 • .

Fig. 4 .
Fig. 4. The absolute value of the averaged error in SSA at a wavelength of 0.5 µm on percentage for the value of AOT, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, and 1.0.The bar indicates the standard deviation.

Fig. 7 .Fig. 8 .
Fig. 7. Values of SSA (left panel) and SDF (right panel) for the enhanced coarse mode case.Solid, dashed, and dashed-dotted lines represent true values, retrieved values from version 4, and retrieved values from version 5, respectively.

Fig. 12 .
Fig. 12. Lidar observation at Beijing by NIES lidar from 14:00 UTC to 15:00 UTC on 14 April 2004.The left panel, center panel, and right panel show backscattering intensity at wavelength 0.532 µm, depolarization ratio at wavelength 0.532 µm, and the color ratio of scattering intensities at λ = 1.064 µm and 0.532 µm, respectively.The arrows indicate the time that we checked the SKYNET data (14 April 2004, 09:06:17 UTC).

Fig. 14 .
Fig. 14.Monthly mean and standard deviation of SSA at λ = 0.5 µm before (black) and after (white) screening by all conditions of C1, C2, and C3, and after screening by each condition (C1, C2 or C3) individually (each screened result by C1, C2, or C3 are shown by red, blue, and green, respectively).The top panel shows the result at Pune in 2008 and the bottom panel shows the result at Beijing in 2004.

Fig. 15 .
Fig. 15.Normalized frequency distributions of SSA from SKYNET and AERONET in May and October at Pune and in April and September at Beijing.Upper and middle panels are for SKYNET SSA before and after data screening, respectively, and lower panel is the figure for AERONET SSA.

Fig. 18 .
Fig. 18.Values of Direct Aerosol Radiative Forcing (DARF) derived by SKYNET before and after screening, and by AERONET in spring and autumn at Pune (May and October) and at Beijing (April and September) at the top of atmosphere (TOA) and the bottom of atmosphere (BOA).
SSA values derived by SKYRAD.packversion 4 and version 5 algorithms for the dust-like and rural aerosol types of WCP (Deepak