Vertical profiles of aerosol optical properties and the solar heating rate estimated by combining sky radiometer and lidar measurements

Abstract

The SKYLIDAR algorithm was developed to estimate vertical profiles of aerosol optical properties from sky radiometer (SKYNET) and lidar (AD-Net) measurements. The solar heating rate was also estimated from the SKYLIDAR retrievals. The algorithm consists of two retrieval steps: (1) columnar properties are retrieved from the sky radiometer measurements and the vertically mean depolarization ratio obtained from the lidar measurements and (2) vertical profiles are retrieved from the lidar measurements and the results of the first step. The derived parameters are the vertical profiles of the size distribution, refractive index (real and imaginary parts), extinction coefficient, single-scattering albedo, and asymmetry factor. Sensitivity tests were conducted by applying the SKYLIDAR algorithm to the simulated sky radiometer and lidar data for vertical profiles of three different aerosols, continental average, transported dust, and pollution aerosols. The vertical profiles of the size distribution, extinction coefficient, and asymmetry factor were well estimated in all cases. The vertical profiles of the refractive index and single-scattering albedo of transported dust, but not those of transported pollution aerosol, were well estimated. To demonstrate the performance and validity of the SKYLIDAR algorithm, we applied the SKYLIDAR algorithm to the actual measurements at Tsukuba, Japan. The detailed vertical structures of the aerosol optical properties and solar heating rate of transported dust and smoke were investigated. Examination of the relationship between the solar heating rate and the aerosol optical properties showed that the vertical profile of the asymmetry factor played an important role in creating vertical variation in the solar heating rate. We then compared the columnar optical properties retrieved with the SKYLIDAR algorithm to those produced with the more established scheme SKYRAD.PACK, and the surface solar irradiance calculated from the SKYLIDAR retrievals was compared with pyranometer measurement. The results showed good agreements: the columnar values of the SKYLIDAR retrievals agreed with reliable SKYRAD.PACK retrievals, and the SKYLIDAR retrievals were sufficiently accurate to evaluate the surface solar irradiance.

Introduction

Aerosols significantly affect the Earth's radiation budget by scattering and absorbing incoming solar radiation (direct effect) and by modifying cloud droplet size and number density (indirect effect). In addition, solar heating of absorbing aerosols such as black carbon or dust affects the vertical profile of the temperature and the cloud cover (semidirect effect; Hansen et al., 1997). To better understand these effects, it is essential to investigate the spatial and temporal variability of the microphysical and optical properties of aerosols, in particular the vertical profiles of aerosol optical properties and solar heating of aerosols. The semidirect effect depends on the vertical profiles of aerosols relative to cloud height (Koch and Del Genio, 2010). Under cloudless conditions, the vertical profile of solar heating of aerosols affects the evolution of the atmospheric boundary layer (Yu et al., 2002; Tsunematsu et al., 2006).

To evaluate solar heating of aerosols, vertical profiles of the extinction coefficient (loading), single-scattering albedo (ratio of scattering to scattering $+$ absorption), and the phase function (or asymmetry factor, i.e., the asymmetry of forward and backward scattering) are necessary. The columnar properties of these parameters, but never their vertical profiles, are routinely obtained by observational networks of sun–sky-scanning multi-wavelength photometers, such as AERONET (Holben et al., 1998) and SKYNET (Takamura and Nakajima, 2004; Nakajima et al., 2007). Active remote sensing by conventional elastic backscatter lidar (MIEL) can provide vertical profiles of the extinction coefficient if the value of the extinction-to-backscatter ratio is assumed but not those of the single-scattering albedo and phase function. Particle extinction and backscatter coefficients can be obtained by using Raman lidar (Ansman et al., 1992) and high-spectra-resolution lidar (HSRL; Shipley et al., 1983) measurements without any assumptions being necessary, and several studies have developed methods for estimating the vertical profiles of the aerosol size distribution and the real and imaginary parts of the refractive index from multi-wavelength Raman lidar data (Müller et al., 1999a, b, 2000; Böckmann, 2001; Veselovskii et al., 2002). Then, on the basis of single-scattering theory, the extinction coefficient, single-scattering albedo, and phase function can be calculated from the size distribution and refractive index.

Synergistic algorithms that relate sun–sky-scanning photometer and MIEL data have been developed. The LiRIC (Lidar-Radiometer Inversion Code) algorithm (Chaikovsky et al., 2012; Wagner et al., 2013) uses the size distribution and refractive index of AERONET retrievals and estimates the vertical profiles of the fine and coarse modes of the size distribution from lidar data. Lopatin et al. (2013) developed the GARRLiC (Generalized Aerosol Retrieval from Radiometer and Lidar Combined data) algorithm by extending the LiRIC algorithm. GARRLiC separately estimates the columnar values of the refractive indices of the fine and coarse modes. Cuesta et al. (2008) developed the LidAlm (Lidar and Almucantar) algorithm, which estimates the vertical profiles after decomposing the AERONET size distribution into multiple lognormal modes. These algorithms provide good estimates of the aerosol vertical profiles from the lidar measurements, based on the AERONET retrievals.

SKYNET and AD-Net (Sugimoto et al., 2005) are dense observational networks of the sun–sky-scanning multi-wavelength photometer (sky radiometer, SKYR) and MIEL, respectively, in the East Asian region, which is one of the world's major sources of dust and anthropogenic aerosols. The synergetic algorithms are useful for observing such aerosols. We developed a new algorithm, called SKYLIDAR, to estimate the vertical profiles of aerosol optical properties from the combination of SKYR and MIEL measurements. Similar to the above-mentioned synergetic algorithms, we assumed a bimodal size distribution but we challenged it to estimate the vertical profiles of the refractive index, which are sensitive to lidar measurements and necessary for determining single-scattering albedo and asymmetry factor. In addition, we estimated the solar heating rate from the SKYLIDAR retrievals.

Schematic diagrams of the retrieval procedures: (a) step 1, (b) step 2, and (c) updating parameters $\mathit{x}$.

[Figure omitted. See PDF]

The National Institute for Environment Studies (NIES) is developing the HSRL at 355 and 532 nm for updating AD-Net (Liu et al., 1999, 2002; Nishizawa et al., 2012). For the future SKYNET and AD-Net networks, we designed the SKYLIDAR algorithm so that it could be applied to both MIEL and HSRL measurements, and the algorithm was tested using the simulated HSRL data in this study.

We describe the SKYR, MIEL, and HSRL measurements, the SKYLIDAR algorithm, and the procedure for calculating the solar heating rate in Sect. 2. In Sect. 3, we evaluate the performance of the SKYLIDAR algorithm by sensitivity tests performed with simulated SKYR, MIEL, and HSRL data for three different vertical profiles of aerosols. In Sect. 4, we apply the algorithm to actual SKYR and MIEL measurements obtained at Tsukuba, Japan, estimating the optical properties and the solar heating rate of transported dust and smoke. We also compare the columnar optical properties of the SKYLIDAR retrievals during 2012 and 2013 with SKYRAD.PACK retrievals, and we compare surface solar irradiances calculated from the SKYLIDAR retrievals with those measured by pyranometer. In Sect. 5, we summarize the results.

Method

The SKYLIDAR algorithm uses the solar direct and diffuse radiations measured by the sun photometer SKYR in the SKYNET and the attenuated backscatter coefficient and the total depolarization ratio by the lidar MIEL in the AD-Net. In addition, the algorithm can be applicable to the HSRL, which is being developed by the NIES. The details of these instruments are described in the beginning of this section.

The SKYLIDAR algorithm (Fig. 1) consists of two retrieval steps. In step 1 (Fig. 1a), the columnar microphysical and optical properties of aerosols are estimated from SKYR measurements and the vertically mean depolarization ratio obtained from MIEL measurements. In step 2 (Fig. 1b), the vertical profiles of the microphysical and optical properties are estimated from MIEL (and HSRL) measurements and the columnar properties obtained in step 1.

SKYLIDAR algorithm

Sky radiometer and lidar measurements

The SKYR (Prede Co., Ltd, Tokyo, Japan), deployed in the SKYNET, is a scanning photometer that measures direct solar radiation and the angular distributions of the diffuse radiation in solar almucantar or principal plane geometries at wavelengths of 315, 340, 380, 400, 500, 675, 870, 940, 1020, 1627, and 2200 nm. Aerosol optical thickness is obtained from the direct solar radiation measurement by using a calibration constant determined by the improved Langley method (Nakajima et al., 1996; Tanaka et al., 1986). Diffuse radiance is obtained from the measured diffuse radiation and the field of view (solid view angle), which is determined by scanning the distribution of radiation around the solar disk (Nakajima et al., 1996). Our algorithm uses aerosol optical thickness ${\mathit{\tau }}_{\mathrm{ext}}\left(\mathit{\lambda }\right)$ and the diffuse radiance normalized by the direct solar radiation $I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$ at wavelength $\mathit{\lambda }$ and at scattering angle $\mathrm{\Theta }$ in the solar almucantar geometry (Fig. 2). The normalized diffuse radiance $I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$ is defined as $$\begin{array}{}& & I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)\equiv I\left(\mathit{\theta },,,\mathit{\varphi },,,\mathit{\lambda }\right)={\displaystyle \frac{F_{\mathrm{dif}}\left(\mathit{\theta },,,\mathit{\varphi },,,\mathit{\lambda }\right)}{F_{\mathrm{dir}}\left(\mathit{\lambda }\right)m\mathrm{\Delta }\mathrm{\Omega }\left(\mathit{\lambda }\right)}},& & \cos \mathrm{\Theta }=\cos \mathit{\theta }\cos {\mathit{\theta }}_{\mathrm{0}}+\sin \mathit{\theta }\sin {\mathit{\theta }}_{\mathrm{0}}\cos \left(\mathit{\varphi },-,{\mathit{\varphi }}_{\mathrm{0}}\right),\end{array}$$ where ${\mathit{\theta }}_{\mathrm{0}}$ and ${\mathit{\varphi }}_{\mathrm{0}}$ are the solar zenith and azimuth angles, $\mathit{\theta }$ and $\mathit{\varphi }$ are the observation zenith and azimuth angles, $\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}$ in the solar almucantar plane, $F_{\mathrm{dif}}\left(\mathit{\theta },,,\mathit{\varphi },,,\mathit{\lambda }\right)$ is the diffuse radiation in the solar almucantar geometry, $m=\mathrm{1}/\cos {\mathit{\theta }}_{\mathrm{0}}$ is the optical air mass, $\mathrm{\Delta }\mathrm{\Omega }\left(\mathit{\lambda }\right)$ is solid view angle, and $F_{\mathrm{dir}}\left(\mathit{\lambda }\right)$ is the direct solar radiation. The wavelength $\mathit{\lambda }$ used in the algorithm is 340, 380, 400, 500, 675, 870, and 1020 nm. The scattering angle $\mathrm{\Theta }$ used in the algorithm is 3, 4, 5, 7, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, and 160${}^{\circ }$ in the solar almucantar plane. Note that the maximum value of the scattering angle depends on the solar zenith angle.

Diffuse radiances scanned in the almucantar geometry are suitable for estimating columnar properties because the effect of the aerosol vertical profile to diffuse radiances on the almucantar plane is weak (Torres et al., 2014). Because the SKYR measures both direct and diffuse radiation with the same detector, the calibration constant cancels out by normalizing diffuse radiance to direct radiation, and the calibration uncertainty is neglected in $I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$.

The SKYLIDAR algorithm can be applied to both MIEL and HSRL measurements developed by NIES. We use the attenuated backscatter coefficients for total (particulate $+$ molecular) scattering at 532 and 1064 nm, and the total depolarization ratio at 532 nm in MIEL measurements. These data are routinely calibrated (Shimizu et al., 2004, 2010). The HSRL measurements are the attenuated backscatter coefficients for molecular scattering at 355 and 532 nm.

The lidar signals from near the surface contain errors due to the incomplete overlap between the transmitted laser beams and the receiver field of view. We therefore excluded the original MIEL and HSRL data from below the altitude of 300 m and used data extrapolated linearly from measurements obtained above 300 m in the algorithm.

Observation of the sky radiometer in the solar almucantar geometry.

[Figure omitted. See PDF]

In many lidar retrievals, the attenuated backscatter coefficient is normalized by the attenuated backscatter coefficient at a reference altitude, where aerosol is neglected. However, the reference altitude is usually high and the reference attenuated backscatter coefficient is influenced by the large noise. Therefore, we normalized the attenuated backscatter coefficients for total and molecular scattering by their vertical means from the surface to top of aerosol layer. The top altitude is determined by the empirically determined threshold of the MIEL backscatter signal at 1064 nm. Because the random noise included in the vertically mean attenuated backscatter coefficient is expected to be much smaller than particulate $+$ molecular backscatter, the calibration constant of the attenuated backscatter coefficient cancels out.

Thus, the SKYLIDAR algorithm can be applied to the SKYR data (${\mathit{\tau }}_{\mathrm{ext}}\left(\mathit{\lambda }\right)$, $I\left(\mathrm{\Theta },,,\mathit{\lambda }\right))$ and the following three data sets of MIEL and HSRL:

• (Type 1) ${\mathit{\beta }}_{\mathrm{MIE}}\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ and 1064 nm and $\mathit{\delta }\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ nm;

• (Type 2) ${\mathit{\beta }}_{\mathrm{MIE}}\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ and 1064 nm, $\mathit{\delta }\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ nm, and ${\mathit{\beta }}_{\mathrm{RAY}}\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ nm;

• (Type 3) ${\mathit{\beta }}_{\mathrm{MIE}}\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=355$, 532, and 1064 nm, $\mathit{\delta }\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=532$ nm, and ${\mathit{\beta }}_{\mathrm{RAY}}\left(\mathit{\lambda },,,z\right)$ at $\mathit{\lambda }=355$ and 532 nm,

where $z$ is altitude, $\mathit{\lambda }$ is wavelength, ${\mathit{\beta }}_{\mathrm{MIE}}\left(\mathit{\lambda },,,z\right)$ is the normalized attenuated backscatter coefficient for total scattering and $\mathit{\delta }\left(\mathit{\lambda },,,z\right)$ is the total depolarization ratio, and ${\mathit{\beta }}_{\mathrm{RAY}}\left(\mathit{\lambda },,,z\right)$ is the normalized attenuated backscatter coefficient for molecular scattering. Type 1 data set is only MIEL measurements, and Type 2 and Type 3 data sets include HSRL measurements.

Step 1

The columnar properties of aerosols are estimated from the SKYR measurements (${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{mea}}\left(\mathit{\lambda }\right)$, $I^{\mathrm{mea}}\left(\mathrm{\Theta },,,\mathit{\lambda }\right))$ and the depolarization ratio averaged from the surface to the top of the aerosol layer (${\mathit{\delta }}_{\mathrm{ave}}^{\mathrm{mea}}\left(532,,\mathrm{nm}\right))$. The aerosol parameters to be estimated are the real and imaginary parts of the refractive index at SKYR wavelengths, the volume size distribution, and the volume ratio of non-spherical particles to total particles in the coarse mode.

Our algorithm estimates the aerosol parameters, based on the maximum a posteriori (MAP) scheme. The aerosol parameters for the best fit to all of the measurements and a priori information are obtained by maximizing the posterior probability density function (PDF): $$P\left(\mathit{x},\mathrm{|},{\mathit{y}}^{\mathrm{mea}}\right)=P\left({\mathit{y}}^{\mathrm{mea}},\mathrm{|},\mathit{x}\right)P\left(\mathit{x}\right)/P\left({\mathit{y}}^{\mathrm{mea}}\right),$$ where vector ${\mathit{y}}^{\mathrm{mea}}$ describes the measurements, vector $\mathit{x}$ describes the aerosol parameters to be estimated, $P\left({\mathit{y}}^{\mathrm{mea}},\mathrm{|},\mathit{x}\right)$ is the conditional PDF of ${\mathit{y}}^{\mathrm{mea}}$ given $\mathit{x}$, $P\left(\mathit{x}\right)$ is the prior PDF of $\mathit{x}$, and $P\left({\mathit{y}}^{\mathrm{mea}}\right)$ is the prior PDF of ${\mathit{y}}^{\mathrm{mea}}$. Since $P\left({\mathit{y}}^{\mathrm{mea}}\right)$ does not depend on $\mathit{x}$, the PDF to be maximized is $P\left({\mathit{y}}^{\mathrm{mea}},\mathrm{|},\mathit{x}\right)P\left(\mathit{x}\right)$. Assuming the normal distribution, $P\left({\mathit{y}}^{\mathrm{mea}},\mathrm{|},\mathit{x}\right)P\left(\mathit{x}\right)$ is defined as $$\begin{array}{cc}& & P\left({\mathit{y}}^{\mathrm{mea}},\mathrm{|},\mathit{x}\right)P\left(\mathit{x}\right)\propto \exp [-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}& {\left({\mathit{y}}^{\mathrm{mea}},-,\mathit{y},\left(\mathit{x}\right)\right)}^T{\left({\mathbf{W}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathit{y}}^{\mathrm{mea}},-,\mathit{y},\left(\mathit{x}\right)\right)\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x},-,{\mathit{x}}_{\mathrm{a}}\right)}^T{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left(\mathit{x},-,{\mathit{x}}_{\mathrm{a}}\right)],\end{array}$$ where vector $\mathit{y}\left(\mathit{x}\right)$ comprises the values corresponding to ${\mathit{y}}^{\mathrm{mea}}$ calculated from $\mathit{x}$ by the forward model, and matrix ${\mathbf{W}}^{\mathrm{2}}$ is the covariance matrix of $\mathit{y}$ and is assumed to be diagonal in this study. The diagonal elements of $\mathbf{W}$ are the standard errors in the measurements. The vector ${\mathit{x}}_{\mathrm{a}}$ comprises the a priori value of $\mathit{x}$, and matrix ${\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}$ is the associated covariance matrix. The maximum of Eq. (4) is obtained by minimizing the objective function: $$\begin{array}{cc}& & f\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}^{\mathrm{mea}},-,\mathit{y},\left(\mathit{x}\right)\right)}^T{\left({\mathbf{W}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathit{y}}^{\mathrm{mea}},-,\mathit{y},\left(\mathit{x}\right)\right)& +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{x},-,{\mathit{x}}_{\mathrm{a}}\right)}^T{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left(\mathit{x},-,{\mathit{x}}_{\mathrm{a}}\right).\end{array}$$ We search for the best $\mathit{x}$, which minimizes $f\left(\mathit{x}\right)$, by iterations of the Gauss–Newton method with a line search, ${\mathit{x}}_{i+\mathrm{1}}={\mathit{x}}_i+{\mathit{\alpha }}_j{\mathit{d}}_i$ (“Update $\mathit{x}$” in Fig. 1). This minimization procedure is described in Sects. 2.1.4 and 2.1.5.

In step 1, the vector ${\mathit{y}}^{\mathrm{mea}}$ consists of ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{mea}}\left(\mathit{\lambda }\right)$ at SKYR wavelengths, $I^{\mathrm{mea}}\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$ at SKYR wavelengths and scattering angles, and ${\mathit{\delta }}_{\mathrm{ave}}^{\mathrm{mea}}$ at 532 nm. Vector $\mathit{x}$ describes the aerosol parameters and consists of the real and imaginary parts of the refractive index at SKYR wavelengths ($n\left(\mathit{\lambda }\right)$ and $k\left(\mathit{\lambda }\right))$, and the parameters for the size distribution ($C_{\mathrm{1}}$, $C_{\mathrm{2}}$, $r_{m,\mathrm{1}}$, $r_{m,\mathrm{2}}$, $s_{\mathrm{1}}$, $s_{\mathrm{2}}$, and $\mathit{\epsilon })$. We assumed the following bi-modal lognormal size distribution: $$\begin{array}{cc}& & {\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{dln}r}}=\sum _{i=\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{\mathrm{d}{V}_i\left(r\right)}{\mathrm{dln}r}}& =\sum _{i=\mathrm{1}}^{\mathrm{2}}{\displaystyle \frac{C_i}{\sqrt{\mathrm{2}\mathit{\pi }}{s}_i}}\exp \left[-,{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},{\left({\displaystyle \frac{\ln r-\ln r_{m,i}}{s_i}}\right)}^{\mathrm{2}}\right],\end{array}$$ where $C_i$, $r_{m,i}$, and $s_i$ are volume, radius, and width, respectively, of the fine ($i=\mathrm{1}$) and coarse ($i=\mathrm{2}$) modes. $\mathit{\epsilon }$ is the volume ratio of non-spherical particles to total particles in the coarse mode.

We constructed the forward models $\mathit{y}\left(\mathit{x}\right)$ to calculate ${\mathit{\tau }}_{\mathrm{ext}}\left(\mathit{\lambda }\right)$, $I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$, and ${\mathit{\delta }}_{\mathrm{ave}}\left(532,,\mathrm{nm}\right)$ from the above-mentioned aerosol parameters. The optical properties of aerosols were calculated by a method similar to that of Dubovik et al. (2006) as follows: $$\begin{array}{cc}& & {\mathit{\tau }}_{\mathrm{ext}/\mathrm{sca}}\left(\mathit{\lambda }\right)=\sum _k{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{1}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{\mathrm{ext}/\mathrm{sca}}^{\mathrm{S}}\left(\mathit{\lambda },,,n,,,k,,,r_k\right)& +\sum _k\left(\mathrm{1},-,\mathit{\epsilon }\right){\displaystyle \frac{\mathrm{d}{V}_{\mathrm{2}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{\mathrm{ext}/\mathrm{sca}}^{\mathrm{S}}\left(\mathit{\lambda },,,n,,,k,,,r_k\right)\\ & +\sum _k\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{2}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{\mathrm{ext}/\mathrm{sca}}^{\mathrm{NS}}\left(\mathit{\lambda },,,n,,,k,,,r_k\right),\end{array}$$ $$\begin{array}{cc}& & {\mathit{\tau }}_{\mathrm{sca}}\left(\mathit{\lambda }\right)P_{ii}\left(\mathrm{\Theta },,,\mathit{\lambda }\right)=\sum _k{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{1}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{ii}^{\mathrm{S}}\left(\mathrm{\Theta },,,\mathit{\lambda },,,n,,,k,,,r_k\right)& +\sum _k\left(\mathrm{1},-,\mathit{\epsilon }\right){\displaystyle \frac{\mathrm{d}{V}_{\mathrm{2}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{ii}^{\mathrm{S}}\left(\mathrm{\Theta },,,\mathit{\lambda },,,n,,,k,,,r_k\right)\\ & +\sum _k\mathit{\epsilon }{\displaystyle \frac{\mathrm{d}{V}_{\mathrm{2}}\left(r_k\right)}{\mathrm{dln}r}}{K}_{ii}^{\mathrm{NS}}\left(\mathrm{\Theta },,,\mathit{\lambda },,,n,,,k,,,r_k\right),\end{array}$$ where ${\mathit{\tau }}_{\mathrm{ext}/\mathrm{sca}}\left(\mathit{\lambda }\right)$ denotes the optical thickness for extinction and scattering, and ${\mathit{\tau }}_{\mathrm{sca}}\left(\mathit{\lambda }\right)P_{ii}\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$ denotes the directional scattering corresponding to the scattering matrix elements $P_{ii}\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$. $K_{\mathrm{\cdots }}^{\mathrm{S}}$ and $K_{\mathrm{\cdots }}^{\mathrm{NS}}$ are the kernels of extinction and scattering properties for spherical and non-spherical particles, respectively. The kernel for spherical particles was constructed by using Mie theory. The kernel for non-spherical particles was constructed by using the data table of Dubovik et al. (2006), which assumes randomly oriented polydisperse spheroids with a fixed aspect ratio distribution for mineral dust.

$I\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$ is computed by the radiative transfer code in the SKYRAD.PACK version 4.2 (Nakajima et al., 1996). Although only one atmospheric layer is considered by SKYRAD.PACK, we assumed that the atmosphere consists of two layers. The bottom layer includes aerosols, and its top altitude $z_{max}$ is determined by lidar measurement. The upper layer is aerosol free. The Rayleigh scattering is calculated by the method of Bucholtz (1995). The vertical ozone profile is approximated by the formula of Green (1964), and the ozone absorption coefficient is adopted from the LOWTRAN 7 database (Kneizys et al., 1988).

${\mathit{\delta }}_{\mathrm{ave}}\left(\mathit{\lambda }\right)$ is calculated as $$\begin{array}{cc}& & {\mathit{\delta }}_{\mathrm{ave}}\left(\mathit{\lambda }\right)=\left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda }\right){\mathit{\delta }}_{\mathrm{m}}\left(\mathit{\lambda }\right)}{\mathrm{1}+{\mathit{\delta }}_{\mathrm{m}}\left(\mathit{\lambda }\right)}},+,{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{p}}\left(\mathit{\lambda }\right){\mathit{\delta }}_{\mathrm{p}}\left(\mathit{\lambda }\right)}{\mathrm{1}+{\mathit{\delta }}_{\mathrm{p}}\left(\mathit{\lambda }\right)}}\right)/& \left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda }\right)}{\mathrm{1}+{\mathit{\delta }}_{\mathrm{m}}\left(\mathit{\lambda }\right)}},+,{\displaystyle \frac{{\mathit{\beta }}_{\mathrm{p}}\left(\mathit{\lambda }\right)}{\mathrm{1}+{\mathit{\delta }}_{\mathrm{p}}\left(\mathit{\lambda }\right)}}\right),\\ & & {\mathit{\beta }}_{\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)={\mathit{\tau }}_{\mathrm{sca},\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)P_{11,\mathrm{p}/\mathrm{m}}\left({180}^{\circ },,,\mathit{\lambda }\right),& & {\mathit{\delta }}_{\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}-P_{22,\mathrm{p}/\mathrm{m}}\left({180}^{\circ },,,\mathit{\lambda }\right)/P_{11,\mathrm{p}/\mathrm{m}}\left({180}^{\circ },,,\mathit{\lambda }\right)}{\mathrm{1}+P_{22,\mathrm{p}/\mathrm{m}}\left({180}^{\circ },,,\mathit{\lambda }\right)/P_{11,\mathrm{p}/\mathrm{m}}\left({180}^{\circ },,,\mathit{\lambda }\right)}},\end{array}$$ where ${\mathit{\beta }}_{\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)$, ${\mathit{\delta }}_{\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)$, and ${\mathit{\tau }}_{\mathrm{sca},\mathrm{p}/\mathrm{m}}\left(\mathit{\lambda }\right)$ are the backscatter coefficient, the depolarization ratio, and the scattering optical thickness, respectively, for particulate (p) and molecular (m) scattering in the bottom aerosol layer.

The first term of Eq. (5) is calculated by using the above-mentioned ${\mathit{y}}^{\mathrm{mea}}$ and $\mathit{y}\left(\mathit{x}\right)$. The values of the diagonal matrix $\mathbf{W}$ are the measurement errors, which were assumed to be 0.01 for ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{mea}}\left(\mathit{\lambda },\ge ,500,,\mathrm{nm}\right)$ and 0.02 for ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{mea}}\left(\mathit{\lambda },<,500,,\mathrm{nm}\right)$, 5 % for $I^{\mathrm{mea}}\left(\mathrm{\Theta },,,\mathit{\lambda }\right)$, and 20 % for ${\mathit{\delta }}_{\mathrm{ave}}^{\mathrm{mea}}\left(532,,\mathrm{nm}\right)$.

To reduce the effects of measurement errors on retrievals, Dubovik and King (2000) restricted the spectral variability of the refractive index by limiting the length of the derivatives of the refractive index with respect to the wavelength. They considered this a priori smoothness constraint as being of the same nature as a measurement and incorporated the smoothness constraint into their retrieval scheme. We adopted this consideration and introduced the smoothness constraint for the real and imaginary parts of the refractive index. The first derivatives of the refractive index with respect to the wavelengths are defined as $$\begin{array}{cc}& & {\mathit{y}}_{\mathrm{a}}\left(\mathit{x}\right)=\left({\displaystyle \frac{\ln n\left({\mathit{\lambda }}_{\mathrm{1}}\right)-\ln n\left({\mathit{\lambda }}_{\mathrm{2}}\right)}{\ln {\mathit{\lambda }}_{\mathrm{1}}-\ln {\mathit{\lambda }}_{\mathrm{2}}}},,,\mathrm{\cdots }\right.& {\displaystyle \frac{\ln n\left({\mathit{\lambda }}_{\mathrm{6}}\right)-\ln n\left({\mathit{\lambda }}_{\mathrm{7}}\right)}{\ln {\mathit{\lambda }}_{\mathrm{6}}-\ln {\mathit{\lambda }}_{\mathrm{7}}}}{\displaystyle \frac{\ln k\left({\mathit{\lambda }}_{\mathrm{1}}\right)-\ln k\left({\mathit{\lambda }}_{\mathrm{2}}\right)}{\ln {\mathit{\lambda }}_{\mathrm{1}}-\ln {\mathit{\lambda }}_{\mathrm{2}}}}\mathrm{\cdots }\\ & \left..{\displaystyle \frac{\ln k\left({\mathit{\lambda }}_{\mathrm{6}}\right)-\ln k\left({\mathit{\lambda }}_{\mathrm{7}}\right)}{\ln {\mathit{\lambda }}_{\mathrm{6}}-\ln {\mathit{\lambda }}_{\mathrm{7}}}}\right).\end{array}$$ ${\mathit{y}}_{\mathrm{a}}\left(\mathbf{x}\right)$ is incorporated into $\mathit{y}\left(\mathit{x}\right)$, and the value of ${\mathit{y}}^{\mathrm{mea}}$ corresponding to ${\mathit{y}}_{\mathrm{a}}\left(\mathbf{x}\right)$ is 0. The values entered in the weight matrix $\mathbf{W}$ were 0.2 for the real part and 1.25 for the imaginary part. These values are used in the AERONET retrieval (Dubovik and King, 2000) for constraining the spectral variability of the refractive index to some practically reasonable ranges.

In step 1, $P\left(\mathit{x}\right)$ is assumed to be uniform, and the second term of Eq. (5) is ignored. The objective function (Eq. 5) is minimized by the procedures described in Sects. 2.1.4 and 2.1.5, and the columnar properties of the real and imaginary parts of the refractive index ($n^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and $k^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right))$ at SKYR wavelengths, the size distribution ($C_{\mathrm{1},\mathrm{2}}^{\mathrm{s}\mathrm{1}}$, $r_{\mathrm{1},\mathrm{2}}^{\mathrm{s}\mathrm{1}}$, and $s_{\mathrm{1},\mathrm{2}}^{\mathrm{s}\mathrm{1}})$, and the volume ratio of non-spherical particles to total particles in the coarse mode (${\mathit{\epsilon }}^{\mathrm{s}\mathrm{1}})$ are optimized. The aerosol optical thickness (${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right))$ and single-scattering albedo (${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right))$ at SKYR wavelengths are calculated with Eqs. (6) to (8). $n^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and $k^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ at MIEL wavelengths are calculated by the linear interpolation and extrapolation in the log–log plot. ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and ${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ at MIEL wavelengths are calculated from those with Eqs. (6) to (8). These step 1 results are input to step 2 (see Fig. 1a and b).

Step 2

The vertical profiles of the refractive index at MIEL wavelengths, the size distribution, and the volume ratio of the non-spherical particles to total particles in the coarse mode are optimized to MIEL (and HSRL) measurements and the columnar properties obtained in step 1 by the same strategy. The final outputs of the extinction coefficients, single-scattering albedo, and asymmetry factor at MIEL wavelengths are calculated from the optimized aerosol parameters.

The columnar properties input from step 1 are the aerosol optical thickness ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and the single-scattering albedo ${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ at MIEL wavelengths. The MIEL and HSRL measurements are three data sets described in Sect. 2.1.1. ${\mathit{y}}^{\mathrm{mea}}$ for Type 1 data set consists of ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and ${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ at MIEL wavelengths and the vertical profiles of ${\mathit{\beta }}_{\mathrm{MIE}}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$ and ${\mathit{\delta }}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$ at MIEL wavelengths. ${\mathit{y}}^{\mathrm{mea}}$ for Type 2 and Type 3 data sets consist of ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and ${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ at MIEL wavelengths and the vertical profiles of ${\mathit{\beta }}_{\mathrm{MIE}}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$ and ${\mathit{\delta }}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$ at MIEL wavelengths and ${\mathit{\beta }}_{\mathrm{RAY}}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$ at HSRL wavelengths. The aerosol parameter $\mathit{x}$ contains the vertical profiles of the real and imaginary parts of the refractive index at MIEL wavelengths ($n\left(\mathit{\lambda },,,z\right)$ and $k\left(\mathit{\lambda },,,z\right))$, and the size distribution parameters ($C_{\mathrm{1}}\left(z\right)$, $C_{\mathrm{2}}\left(z\right)$, and $\mathit{\epsilon }\left(z\right))$. The bi-modal size distribution (Eq. 6) is also used in step 2, but the mode radii and the widths of the fine and coarse modes are fixed by the columnar values obtained in step 1 ($r_{\mathrm{1},\mathrm{2}}^{\mathrm{s}\mathrm{1}}$ and $s_{\mathrm{1},\mathrm{2}}^{\mathrm{s}\mathrm{1}})$.

In the forward model $\mathit{y}\left(\mathit{x}\right)$ of step 2, the aerosol optical properties at each altitude are calculated with Eqs. (6) to (8), but note that the extinction/scattering coefficients are calculated.

The normalized attenuated backscatter coefficients for total and molecular scattering are calculated by the lidar equations: $$\begin{array}{cc}& & {\mathit{\beta }}_{\mathrm{MIE}}\left(\mathit{\lambda },,,z\right)=\left({\mathit{\beta }}_{\mathrm{m}},\left(\mathit{\lambda },,,z\right),+,{\mathit{\beta }}_{\mathrm{p}},\left(\mathit{\lambda },,,z\right)\right)& \exp \left(-,\mathrm{2},\underset{\mathrm{0}}{\overset{z}{\int }},{\mathit{\sigma }}_{\mathrm{ext},\mathrm{p}},\left(\mathit{\lambda },,,z^{\prime }\right),+,{\mathit{\sigma }}_{\mathrm{ext},\mathrm{m}},\left(\mathit{\lambda },,,z^{\prime }\right),\mathrm{d},z^{\prime }\right)/{\mathit{\beta }}_{\mathrm{MIE},\mathrm{ave}}\left(\mathit{\lambda }\right),\end{array}$$ $$\begin{array}{cc}& & {\mathit{\beta }}_{\mathrm{RAY}}\left(\mathit{\lambda },,,z\right)={\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda },,,z\right)& \exp \left(-,\mathrm{2},\underset{\mathrm{0}}{\overset{z}{\int }},{\mathit{\sigma }}_{\mathrm{ext},\mathrm{p}},\left(\mathit{\lambda },,,z^{\prime }\right),+,{\mathit{\sigma }}_{\mathrm{ext},\mathrm{m}},\left(\mathit{\lambda },,,z^{\prime }\right),\mathrm{d},z^{\prime }\right)/{\mathit{\beta }}_{\mathrm{RAY},\mathrm{ave}}\left(\mathit{\lambda }\right),\end{array}$$ where ${\mathit{\sigma }}_{\mathrm{ext},\mathrm{p}/\mathrm{m}}$ are the extinction coefficients for particulate (p) and molecular (m) scattering, and ${\mathit{\beta }}_{\mathrm{MIE},\mathrm{ave}}\left(\mathit{\lambda }\right)$ and ${\mathit{\beta }}_{\mathrm{RAY},\mathrm{ave}}\left(\mathit{\lambda }\right)$ are the vertical means of the calculated attenuated backscatter coefficients. The total depolarization ratio ($\mathit{\delta }\left(\mathit{\lambda },,,z\right))$ at each altitude is calculated with Eqs. (9) to (11).

The first term of Eq. (5) is calculated with the above-mentioned ${\mathit{y}}^{\mathrm{mea}}$ and $\mathit{y}\left(\mathit{x}\right)$. The values of the diagonal matrix $\mathbf{W}$ were assumed to be 0.01 for ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda },\ge ,532,,\mathrm{nm}\right)$, 0.02 for ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda },<,532,,\mathrm{nm}\right)$, 0.05 for ${\mathit{\omega }}_{\mathrm{0}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$, 10 % for ${\mathit{\beta }}_{\mathrm{MIE}}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$, 15 % for ${\mathit{\beta }}_{\mathrm{RAY}}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$, and 20 % for ${\mathit{\delta }}^{\mathrm{mea}}\left(\mathit{\lambda },,,z\right)$.

In step 2, the number of estimated parameters is larger than the number of measurements, so the lidar measurements would be insufficient for retrieving unique solutions of the refractive index even if the columnar properties obtained in step 1 were added to ${\mathit{y}}^{\mathrm{mea}}$. Therefore, we gave a priori constraints to the real and imaginary parts of the refractive index using the step 1 results $k^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$ and ${\mathit{\tau }}_{\mathrm{ext}}^{\mathrm{s}\mathrm{1}}\left(\mathit{\lambda }\right)$. Assuming the lognormal PDF for $P\left(\mathit{x}\right)$, the second term of Eq. (5) is replaced by $$\begin{array}{cc}& & {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\ln ,\mathbf{n},-,\ln ,{\mathbf{n}}_{\mathrm{a}}^{\mathrm{s}\mathrm{1}}\right)}^T{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left(\ln ,\mathbf{n},-,\ln ,{\mathbf{n}}_{\mathrm{a}}^{\mathrm{s}\mathrm{1}}\right)& +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\ln ,\mathbf{k},-,\ln ,{\mathbf{k}}_{\mathrm{a}}^{\mathrm{s}\mathrm{1}}\right)}^T{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left(\ln ,\mathbf{k},-,\ln ,{\mathbf{k}}_{\mathrm{a}}^{\mathrm{s}\mathrm{1}}\right).\end{array}$$ The matrix ${\mathbf{W}}_{\mathrm{a}}$, which determines the strength of the constraints, is assumed to be diagonal and the values are obtained by a method similar to that of Dubovik and King (2000). The possible variability ranges of the refractive index for aerosols are from 1.33 to 1.6 for the real part and from 0.0005 to 0.5 for the imaginary part. We considered these intervals to be 68 % confidence intervals, $[\ln n^{\mathrm{s}\mathrm{1}}-\mathrm{\Delta }n$, $\ln n^{\mathrm{s}\mathrm{1}}+\mathrm{\Delta }n$] and [$\ln k^{\mathrm{s}\mathrm{1}}-\mathrm{\Delta }k$, $\ln k^{\mathrm{s}\mathrm{1}}+\mathrm{\Delta }k$], and determined the weight values as $${\mathbf{W}}_{\mathrm{a}}=\left\{\begin{array}{ll}0.5\left(\ln ,n_{max},-,\ln ,n_{min}\right),& \mathrm{for}\mathrm{real}\mathrm{part}\\ 0.5\left(\ln ,k_{max},-,\ln ,k_{min}\right),& \mathrm{for}\mathrm{imaginary}\mathrm{part}\end{array}\right..$$ The objective function is minimized by the procedures described in Sects. 2.1.4 and 2.1.5, and the vertical profiles of the real and imaginary parts of the refractive index ($n\left(\mathit{\lambda },,,z\right)$, $k\left(\mathit{\lambda },,,z\right))$ at MIEL wavelengths, the size distribution ($C_{\mathrm{1}}\left(z\right)$, $C_{\mathrm{2}}\left(z\right))$, and the volume ratio of non-spherical particles to total particles in the coarse mode ($\mathit{\epsilon }\left(z\right))$ are optimized. Finally, the vertical profiles of the refractive index, size distribution $\left(\frac{\mathrm{d}V\left(r,,,z\right)}{\mathrm{dln}r}\right)$, extinction coefficients (${\mathit{\sigma }}_{\mathrm{ext}}\left(\mathit{\lambda },,,z\right))$, single-scattering albedo (${\mathit{\omega }}_{\mathrm{0}}\left(\mathit{\lambda },,,z\right))$, and asymmetry parameter ($g\left(\mathit{\lambda },,,z\right))$ are produced as output. The wavelengths of these optical properties are 532 and 1064 nm for Type 1 and Type 2 data sets, and 355, 532, 1064 nm for Type 3 data set.

Minimization procedure

In both steps 1 and 2, $\mathit{x}$ was optimized to ${\mathit{y}}^{\mathrm{mea}}$ by finding the minimum of $f\left(\mathit{x}\right)$ in $\mathit{x}$ space. We employed the Gauss–Newton method to search for the minimum. The Gauss–Newton method searches for the minimum by iteratively updating $\mathit{x}$, and it is usually combined with a line search method (Nocedal and Wright, 2006). In this procedure (Fig. 1c), $\mathit{x}$ is updated by ${\mathit{x}}_{i+\mathrm{1}}={\mathit{x}}_i+{\mathit{\alpha }}_j{\mathit{d}}_i$, where the value of $f\left(\mathit{x}\right)$ in the $\mathit{x}$ space decreases in the vector ${\mathit{d}}_i$ direction, and ${\mathit{\alpha }}_j$ is a positive parameter that minimizes $f\left(\mathit{x}\right)$ in direction ${\mathit{d}}_i$. ${\mathit{d}}_i$ is determined by the Gauss–Newton method in the outer loop, and ${\mathit{\alpha }}_j$ is determined by a line search in the inner loop. ${\mathit{d}}_i$ is obtained by solving the normal equation: $$\begin{array}{cc}& & \left[\mathbf{J},{\left({\mathit{x}}_i\right)}^T,{\left({\mathbf{W}}^{\mathrm{2}}\right)}^{-\mathrm{1}},\mathbf{J},\left({\mathit{x}}_i\right),+,{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\right]{\mathit{d}}_i& =\mathbf{J}{\left({\mathit{x}}_i\right)}^T{\left({\mathbf{W}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathit{y}}^{\mathrm{mea}},-,\mathit{y},\left({\mathit{x}}_i\right)\right)+{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\left({\mathit{x}}_i,-,{\mathit{x}}_{\mathrm{a}}\right),\end{array}$$ where $\mathbf{J}\left({\mathit{x}}_i\right)$ is the Jacobi matrix and is calculated as the first derivatives of $\mathit{y}\left({\mathit{x}}_i\right)$ in the near vicinity of ${\mathit{x}}_i$. We solved this normal equation by singular value decomposition (Press et al., 1992). After ${\mathit{d}}_i$ is determined, ${\mathit{\alpha }}_j$ is searched for by the iteration of ${\mathit{\alpha }}_{j+\mathrm{1}}=\mathit{\eta }{\mathit{\alpha }}_j$. The initial value of ${\mathit{\alpha }}_j$ is 1.0, and the value of $\mathit{\eta }$ is set to 0.5. ${\mathit{\alpha }}_j$ is iteratively decreased until the Armijo condition is satisfied: $$f\left({\mathit{x}}_i,+,{\mathit{\alpha }}_j,{\mathit{d}}_i\right)\le f\left({\mathit{x}}_i\right)+\mathit{\gamma }{\mathit{\alpha }}_j\mathrm{\nabla }f{\left({\mathit{x}}_i\right)}^T{\mathit{d}}_i,\mathrm{0}<\mathit{\gamma }<\mathrm{1},$$ where $\mathit{\gamma }$ is an arbitrary constant that we set to 0.001. This line search inhibits unstable oscillation in the Gauss–Newton method by limiting the length of $∥{\mathit{x}}_{i+\mathrm{1}},-,{\mathit{x}}_i∥$; as a result, monotonic and stable convergence is obtained.

Logarithmic transformation

In our minimization problem, the number of elements in $\mathit{x}$ is on the order of 10${}^{\mathrm{2}}$ in step 2. Furthermore, the elements in $\mathit{x}$ and $\mathit{y}$ have different units and values that vary over a wide range of magnitude. Under these conditions, there are too many iterations of the Gauss–Newton method, and convergence cannot be obtained. It is therefore necessary to make the convergence efficient. Dubovik and King (2000) developed a logarithmic transformation technique for $\mathit{y}$ and $\mathit{x}$ by which $f\left(\mathit{x}\right)$ becomes dimensionless, because the term $\mathit{y}\left(\mathit{x}\right)-{\mathit{y}}^{\mathrm{mea}}$ in Eq. (5) is expressed as $\ln \left(\mathit{y},\left(\mathit{x}\right),/,{\mathit{y}}^{\mathrm{mea}}\right)$. This makes it simple to operate simultaneously with $\mathit{y}$, which has different units and values. Furthermore, the logarithmic transformation prevents $\mathit{x}$ from assuming a negative value. We developed more effective transformation techniques. $y$ is transformed as $$Y=\ln \left(y,-,y_{min}\right),$$ where $y_{min}$ is a possible minimum value of $y$. When the aerosol load is small, the attenuated backscatter coefficient and depolarization ratio measurements can have negative values because of the large random noise. Although the logarithmic transformation of Dubovik and King (2000) cannot be applied to negative values, Eq. (19) can be.

$x$ is transformed by $$X=\ln \left({\displaystyle \frac{x-x_{min}}{x_{max}-x}}\right),$$ where $x_{min}$ and $x_{max}$ are minimum and maximum possible values of $x$, respectively. This equation can be inverted as $$x={\displaystyle \frac{x_{min}+x_{max}\exp \left(X\right)}{\mathrm{1}+\exp \left(X\right)}}.$$ The value of $x$ is sometimes limited for physical or numerical reasons. In this study, for rapid computation, we constructed the look-up tables of $K_{\mathrm{\cdots }}^{\mathrm{S}/\mathrm{NS}}$ in Eqs. (7) and (8) with minimum and maximum values of the refractive index, 1.33 and 1.60 for the real part, and 0.0005 and 0.5 for the imaginary part, respectively. If the refractive index exceeds its maximum or minimum values in the iteration of ${\mathit{x}}_{i+\mathrm{1}}={\mathit{x}}_i+{\mathit{\alpha }}_j{\mathit{d}}_i$, then the optical properties cannot be calculated and the retrieval process stops. The volume ratio of the non-spherical particles to total particles in the coarse mode must be between 0 and 1 for a physical reason. When the volume ratio of the non-spherical particles has negative values, the values of the size distribution may be negative. These situations can happen due to the large noise of the lidar measurements. The transformation by Eq. (20) prevents $x$ from exceeding its minimum or maximum limitations. Furthermore, because $X$ and $Y$ are dimensionless, it is simpler to deal with the multiple parameters simultaneously. The transformations of the real and imaginary parts of the refractive index by Eq. (20) are illustrated in Fig. 3. Although the magnitude of each part varies over a different wide range, the magnitudes of the transformed values vary over the same range. However, the value of $x-x_{min}$ or $x_{max}-x$ sometimes become 0 due to a rounding error when the value of $x$ is close to $x_{min}$ or $x_{max}$. In this case, Eq. (20) cannot be calculated. Therefore, Eq. (20) should be computed carefully. For example, when the real part of the refractive index becomes its maximum value of 1.60, we change the value to 1.5999.

Transformation of the real and imaginary parts of the refractive indices by Eq. (20). Minimum and maximum values of the real part are 1.33 and 1.6 and those of the imaginary part are 0.0005 and 0.5.

[Figure omitted. See PDF]

Microphysical and optical properties and vertical profiles of the aerosols used in the sensitivity tests.

Aerosol	Components	Particle shape	Size distribution		Refractive index at 500 nm		Relative weight in total optical thickness at 500 nm	Vertical profile
			Mode radius ($\mathrm{\mu }$m)	Mode width	Real	Imaginary
Continental average	Water-soluble	Sphere	0.18	0.81	1.44	0.0026	0.90	$\exp \left(-,z,/,H\right)$, $H=\mathrm{8}$ km
	Soot	Sphere	0.05	0.69	1.75	0.45	0.07	$\exp \left(-,z,/,H\right)$, $H=\mathrm{4}$ km
	Insoluble	Spheroid	5.98	0.92	1.53	0.008	0.03	$\exp \left(-,z,/,H\right)$, $H=\mathrm{2}$ km
Transported dust	Dust	Spheroid	3.23	0.79	1.53	0.0078	0.25	$\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}\exp \left(-,\frac{\left(z,-,z_{\mathrm{c}}\right)}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)$, $z_{\mathrm{c}}=3.5$ km, $\mathit{\sigma }=0.4$ km
	Water-soluble	Sphere	0.18	0.81	1.44	0.0026	0.67	$\exp \left(-,z,/,H\right)$, $H=\mathrm{8}$ km
	Soot	Sphere	0.05	0.69	1.75	0.45	0.05	$\exp \left(-,z,/,H\right)$, $H=\mathrm{4}$ km
	Insoluble	Spheroid	5.98	0.92	1.53	0.008	0.03	$\exp \left(-,z,/,H\right)$, $H=\mathrm{2}$ km
Transported pollution	Water-soluble	Sphere	0.18	0.81	1.44	0.0026	0.08	$\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}\exp \left(-,\frac{\left(z,-,z_{\mathrm{c}}\right)}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)$, $z_{\mathrm{c}}=3.5$ km, $\mathit{\sigma }=0.4$ km
	Soot	Sphere	0.05	0.69	1.75	0.45	0.03	$\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }\mathit{\sigma }}}\exp \left(-,\frac{\left(z,-,z_{\mathrm{c}}\right)}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right)$, $z_{\mathrm{c}}=3.5$ km, $\mathit{\sigma }=0.4$ km
	Water-soluble	Sphere	0.18	0.81	1.44	0.0026	0.79	$\exp \left(-,z,/,H\right)$, $H=\mathrm{8}$ km
	Soot	Sphere	0.05	0.69	1.75	0.45	0.06	$\exp \left(-,z,/,H\right)$, $H=\mathrm{4}$ km
	Insoluble	Spheroid	5.98	0.92	1.53	0.008	0.03	$\exp \left(-,z,/,H\right)$, $H=\mathrm{2}$ km

By applying the above transformations, the problem of searching for the minimum value of $f\left(\mathit{x}\right)$ in $\mathit{x}$ space becomes a search for the minimum in dimensionless $X$ space. This transformation has a scaling effect. In minimization problems, the convergence rate of an algorithm becomes more rapid when the problem is well scaled (Nocedal and Wright, 2006). An example of poor scaling is the function, $f\left(x_{\mathrm{1}},,,x_{\mathrm{2}}\right)={10}^{10}{x}_{\mathrm{1}}^{\mathrm{2}}+x_{\mathrm{2}}$, which is sensitive to small changes in $x_{\mathrm{1}}$ but not sensitive to $x_{\mathrm{2}}$. If we define a new variable $X_{\mathrm{1}}={10}^{\mathrm{5}}{x}_{\mathrm{1}}$ and minimize the function in terms of $X_{\mathrm{1}}$ and $x_{\mathrm{2}}$, the optimum values of $X_{\mathrm{1}}$ and $x_{\mathrm{2}}$ can be found more rapidly. The condition number, which is a measure of the scaling, is defined as the ratio of the maximum to the minimum singular value of matrix $\left[\mathbf{J},{\left({\mathit{x}}_i\right)}^T,{\left({\mathbf{W}}^{\mathrm{2}}\right)}^{-\mathrm{1}},\mathbf{J},\left({\mathit{x}}_i\right),+,{\left({\mathbf{W}}_{\mathrm{a}}^{\mathrm{2}}\right)}^{-\mathrm{1}}\right]$ of Eq. (17) in this study. If the condition number is close to 1, the problem is well scaled. The condition number for step 2 in this study is on the order of 10${}^{19}$ for the case without any transformations, 10${}^{\mathrm{4}}$ when the logarithmic transformation of Dubovik and King (2000) is applied, and 10${}^{\mathrm{2}}$ when our transformations are applied. Actually, no step 2 results can be obtained unless transformations are applied.

Retrieval results from the simulated data for the continental average aerosol: the real and imaginary parts of the refractive index, the volume concentrations of the fine and coarse modes, the volume ratio of the non-spherical particles to total particles in the coarse mode, and the size distributions at two altitudes. “True” values are shown by the black lines, and the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm; HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green). The red and blue lines completely overlap.

[Figure omitted. See PDF]

Retrieval results from the simulated data for the continental average aerosol: the extinction coefficients, the single-scattering albedo, and the asymmetry factor. “True” values are shown by the black lines, and the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm, and HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green). The red and blue lines completely overlap.

[Figure omitted. See PDF]

Same as Fig. 4 but showing the retrieval results from the simulated data with random errors for the continental average aerosol. Note that the ranges of the horizontal axes in the plots of the volume concentration of the coarse mode and the volume ratio of the non-spherical particles to total particles in the coarse mode are different from those in Fig. 4.

[Figure omitted. See PDF]

Same as Fig. 5 but showing the retrieval results from the simulated data with random errors for the continental average aerosol.

[Figure omitted. See PDF]

Retrieval results from the simulated data for the transported dust: the real and imaginary parts of the refractive index, the volume concentrations of the fine and coarse modes, the volume ratio of the non-spherical particles to total particles in the coarse mode, and the size distributions at two altitudes. “True” values are shown by the black lines, the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm; HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green).

[Figure omitted. See PDF]

Retrieval results from the simulated data for the transported dust: the extinction coefficients, the single-scattering albedo, and the asymmetry factor. “True” values are shown by the black lines, the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm; HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green).

[Figure omitted. See PDF]

Figure 8 but showing the retrieval results from the simulated data with random errors for the transported dust.

[Figure omitted. See PDF]

Figure 9 but showing the retrieval results from the simulated data with random errors for the transported dust.

[Figure omitted. See PDF]

Retrieval results from the simulated data for the transported pollution aerosol: the real and imaginary parts of the refractive index, the volume concentrations of the fine and coarse modes, the volume ratio of the non-spherical particles to total particles in the coarse mode, and the size distributions at two altitudes. “True” values are shown by the black lines, the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm; HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green). The red and blue lines completely overlap.

[Figure omitted. See PDF]

Retrieval results from the simulated data for the transported pollution aerosol: the extinction coefficients, the single-scattering albedo, and the asymmetry factor. “True” values are shown by the black lines, the colored lines indicate retrievals from Type 1 data set (MIEL data at 532 and 1064 nm) (blue), Type 2 data set (MIEL data at 532 and 1064 nm; HSRL data at 532 nm) (red), and Type 3 data set (MIEL data at 532 and 1064 nm; HSRL data at 355 and 532 nm) (green). The red and blue lines completely overlap.

[Figure omitted. See PDF]

Same as Fig. 12 but showing the retrieval results from the simulated data with random errors for the transported pollution aerosol.

[Figure omitted. See PDF]

Same as Fig. 13 but showing the retrieval results from the simulated data with random errors for the transported pollution aerosol.

[Figure omitted. See PDF]

Lidar ratio as a function of the mode radius of the mono-modal lognormal size distribution. The mode width is fixed by 0.5. The imaginary part of the refractive index is 0.01 in the left panels. The real part of the refractive index is 1.50 in the right panels.

[Figure omitted. See PDF]

Retrieval results for the transported dust on 2 April 2012 in Tsukuba, Japan: (a) extinction coefficient at 532 nm, (b) single-scattering albedo at 532 nm, (c) asymmetry factor at 532 nm, (d) size distribution, (e) solar heating rate, and (f) solar heating rate difference from that calculated using the vertical mean of the asymmetry factor. The gray shading shows the attenuated backscatter coefficient at 532 nm by MIEL.

[Figure omitted. See PDF]

Solar heating rate

The vertical profiles of the solar heating rate are calculated from the vertical profiles of the extinction coefficient, single-scattering albedo, and phase function in the solar wavelength region by using the radiative transfer model. However, the wavelengths at which the optical properties obtained by the SKYLIDAR algorithm are 532 and 1064 nm for Type 1 data set. We calculated the refractive index between 532 and 1064 nm by the linear interpolation in a log–log plane. For wavelengths less than 532 nm and greater than 1064 nm, the refractive index at 532 and 1064 nm were used. The extinction coefficient, single-scattering albedo, and phase function in the broadband wavelength regions were calculated from the linearly estimated refractive index and the size distribution. The radiative transfer calculation was performed with our developed code (Asano and Shiobara, 1989; Nishizawa et al., 2004; Kudo et al., 2011). The solar spectrum between 300 and 3000 nm was divided into 54 intervals. Gaseous absorption by water vapor, carbon dioxide, oxygen, and ozone were included in the radiative transfer model.

Sensitivity tests using simulated data

To evaluate the performance of the SKYLIDAR algorithm, we conducted sensitivity tests using simulated SKYR, MIEL, and HSRL data (types 1 to 3 data sets). Three aerosol vertical profile patterns were used for the simulation: (1) continental average, (2) continental average $+$ transported dust in the upper air, and (3) continental average $+$ transported pollution aerosol in the upper air. The microphysical and optical properties and the vertical profiles of the continental average, transported dust, and transported pollution aerosols are summarized in Table 1. The continental average was defined as an external mixture of water-soluble particles, soot particles, and insoluble particles (Hess et al., 1999). The pollution aerosol in this sensitivity test was defined as an external mixture of water-soluble and soot particles. The kernels of the non-spherical particles described in Sect. 2.1.2 were used for calculating the optical properties of the insoluble and dust particles. The sensitivity tests were conducted for aerosol optical thickness of 0.05, 0.1, 0.3, 0.5, 0.8, and 1.2 at 500 nm.

Figures 4 and 5 illustrate the retrieval results from the simulated data for the continental average aerosol with an aerosol optical thickness of 0.05 at 500 nm. The “true” values of the refractive index for the externally mixed aerosols were calculated as the weighted arithmetic mean with the scattering coefficient of each aerosol component. The retrieval results with and without HSRL data were the almost same. The retrieval errors of the refractive index were small except for the real part at 355 nm. Although the volume concentrations of the fine and coarse modes were overestimated, the size distributions at two altitudes were reproduced. However, the coarse mode of the size distribution at lower altitude was overestimated. The volume ratio of the non-spherical particles to total particles in the coarse mode was estimated well except for the result for Type 3 data set. The vertical profiles of the extinction coefficients, the single-scattering albedo, and the asymmetry factor were reproduced. There were the bias errors in the retrieved volume concentrations, but their influences to the extinction coefficients were small because the extinction coefficient depends on not the volume but the cross sectional area. The bias errors of the volume concentrations decreased in the case of the larger aerosol optical thickness (not shown).

We also conducted sensitivity tests using simulated data with random errors to investigate the performance of the algorithm under more realistic conditions. The random error distribution was uniform, and the minimum and maximum values of the random errors were $\pm$2 % for direct solar radiation, $\pm$3 % for diffuse sky radiances, $\pm$5 % for the attenuated backscatter coefficient for total scattering, $\pm$10 % for the attenuated backscatter coefficient for molecular scattering, and $\pm$15 % for the total depolarization ratio. Figures 6 and 7 illustrate the retrieval results from the simulated data for the continental average aerosol, but the simulated data include random errors. The influences of the random errors to the refractive index were inhibited by the a priori constraints to the refractive index in the step 2 retrieval. The volume concentrations of the fine and coarse modes, and the volume ratio of the non-spherical particles were significantly influenced by the random errors. However, the retrieval errors of the extinction coefficients were small by the reason mentioned in the previous paragraph. The retrieval results of the single-scattering albedo and asymmetry factor exhibited large oscillations, so their vertical profiles were not clear. These oscillations are caused by the retrieval errors in the volume concentrations of the fine and coarse modes because the single-scattering albedo and asymmetry factor depend on not only the refractive index but also the shape of the size distribution.

Figures 8 and 9 present the retrieval results from the simulated data for the transported dust with an aerosol optical thickness of 0.5 at 500 nm. The vertical profiles of the refractive index, except the real part at 355 and 532 nm, were well estimated. The vertical profiles of the volume concentrations, the size distributions, and the volume ratio of the non-spherical particles were also reproduced. As a result, the vertical profiles of the extinction coefficients, the single-scattering albedo, and the asymmetry factor were reproduced well. Figures 10 and 11 show the retrieval results from the simulated data with random errors. There were small oscillations in the vertical profiles of all the retrievals but the results were almost same as those in Figs. 8 and 9.

Figures 12 and 13 show the retrieval results from the simulated data for the transported pollution aerosol with an aerosol optical thickness of 0.3 at 500 nm. The vertical profiles of the volume concentrations, the size distributions, and the real part of the refractive index were estimated well, but the vertical profile of the imaginary part of the refractive index was not. The volume ratio of the non-spherical particles was underestimated at the upper altitude where the volume of the coarse mode was very small. In this test, the large values of the imaginary part of the refractive index and the small values of the single-scattering albedo at upper altitudes were important characteristics that were not reproduced even when HSRL data was used in the retrieval. Figures 14 and 15 are the retrieval results from the simulated data with random errors. All the vertical profiles had small oscillations but were almost the same as those in Figs. 12 and 13.

Overall, in these tests our algorithm estimated well the vertical profiles of the size distribution. The vertical profiles of the refractive index were also reproduced, but the real parts of the refractive index of the transported dust and the imaginary part of the refractive index of the transported pollution aerosol were not. The vertical profile of the extinction coefficient was reproduced well because the size distribution is a most important factor to determine the value of the extinction coefficient. The vertical profile of the asymmetry factor was also reproduced because the asymmetry factor depends on the shape of the size distribution firstly and the real part of the refractive index secondly. The vertical profile of the single-scattering albedo was reproduced in the case of transported dust because the shape of the size distribution and the imaginary part of the refractive index were reproduced. However, the vertical profiles of the single-scattering albedo and the imaginary part of the refractive index were not in the case of transported pollution aerosol. These characteristics were consistently observed regardless of the aerosol optical thickness value of the simulated data.

To investigate the difference of the retrievals of the refractive index between the transported dust and pollution aerosol, we calculated the dependencies of the extinction-to-backscatter ratio on the refractive index and the mode radius (Fig. 16). The mono-modal lognormal size distribution (Eq. 6) with a mode width of 0.5 was used in the calculation. The optical properties of the randomly oriented polydisperse spheroids described in Sect. 2.1.2 were used in the case of the non-spherical particle. The extinction-to-backscatter ratios at all the wavelengths were widely changed by the real part of the refractive index in the results of the coarse mode (the mode radius from 1 to 5 $\mathrm{\mu }$m), but such large variations were not seen in the results of the fine mode (the mode radius from 0.1 to 0.2 $\mathrm{\mu }$m). The extinction-to-backscatter ratios depended on the imaginary part of the refractive index in the results of both the fine and coarse modes. However, for the fine mode, the dependencies were small in the limited range of the imaginary part from 0.005 to 0.02. This is a range of our defined transported pollution aerosol (Fig. 12). Consequently, the extinction-to-backscatter ratio is sensitive to the real and imaginary parts of the refractive index in the case that the coarse mode is dominant but is not in the case of the fine mode. These features were seen in the results of both the spherical and non-spherical particles.

The retrieval results obtained with and without HSRL data did not differ. The advantage of HSRL data is that particle backscatter and extinction coefficients are obtained separately. Our algorithm in step 1 also estimates the optical thickness and the other parameters related to the backscatter coefficient separately. The retrievals of step 1 are the columnar values, but manipulating the MIEL data with the step 1 retrievals would have an effect similar to the addition of HSRL data. In this regard, our algorithm cannot utilize HSRL data; thus further development of the algorithm is necessary.

The influence of the random error was small when the optical thickness was large. However, the random errors had a large influence when the aerosol optical thickness was small (Figs. 6 and 7). In the methods of Chaikovsky et al. (2012) and Lopatin et al. (2013), the smoothness constraints for the vertical profiles of the volumes in the fine and coarse modes are introduced to overcome the random noise of the lidar measurements. Since this smoothness constraint is expected to decrease the large oscillations of the single-scattering albedo and the asymmetry factor in the optically thin case, we will introduce the smoothness constraint in the future.

Two-day backward trajectories of air masses observed on 2 April 2012. The ochre region and orange areas in the upper panel indicate desert region and urban regions, respectively.

[Figure omitted. See PDF]

Retrieval results for the transported smoke in 8 May 2013 in Tsukuba, Japan: (a) extinction coefficient at 532 nm, (b) single-scattering albedo at 532 nm, (c) asymmetry factor at 532 nm, (d) size distribution, (e) solar heating rate, and (f) solar heating rate difference from that calculated using the vertical mean of the asymmetry factor. The gray shading is the attenuated backscatter coefficient at 532 nm by MIEL.

[Figure omitted. See PDF]

Two-day backward trajectory of air mass observed on 8 May 2013. The upper panel also shows fire activity from 1 to 9 May 2013 (color scale), based on data of MODIS active-fire products (NEO, 2016).

[Figure omitted. See PDF]

Comparisons of the columnar optical properties estimated by SKYLIDAR algorithm (this study) and SKYRAD.PACK during 2012 and 2013 in Tsukuba, Japan. The colors in the single-scattering albedo and asymmetry factor panels show the aerosol optical thickness at 532 nm. In normalized volume size distribution panel, red indicates the SKYLIDAR result, and blue indicates the SKYRAD.PACK result.

[Figure omitted. See PDF]

Application to observational data

We applied our developed SKYLIDAR algorithm to actual SKYR and MIEL measurements obtained during 2012 and 2013 at Tsukuba (140.12${}^{\circ }$ E, 36.05${}^{\circ }$ N), Japan.

Application of the SKYLIDAR algorithm to actual measurements

The retrieval results of the aerosol optical properties and the solar heating rate for transported dust observed in 2 April 2012 are shown in Fig. 17. High values of the extinction coefficients were observed in layer 1 (1–2 km), layer 2 (2–3 km), and layer 3 (3–5 km) (Fig. 17a). The high extinction coefficient in layer 1 was attributed to aerosols remaining in the residual layer from the preceding day. The aerosols in layers 2 and 3 were identified as transported pollution aerosol and dust from China, respectively, because the backward trajectories (Fig. 18) indicated that the layer 2 aerosol had been transported from an urban region, and the layer 3 aerosol had been transported from a desert region. These backward trajectories were calculated with our trajectory model, which was developed following Katsumoto et al. (2002). The temporal and spatial distributions of three-dimensional winds were linearly interpolated from the US National Centers for Environmental Prediction 6-hourly reanalysis data set (Kalnay et al., 1996). The time integration was conducted by the fourth-order Runge–Kutta method. The desert and urban regions shown in Fig. 18 were determined by using data from the Land Cover Type Climate Modeling Grid product (LP DAAC, 2013).

The aerosols in layer 3 had a large asymmetry factor value of more than 0.7 (Fig. 17c), and the coarse mode of the size distribution was dominant (Fig. 17d). Therefore, we interpreted the aerosol in layer 3 as pure dust. The asymmetry factor of the aerosols in layer 2 was about 0.65 or higher (Fig. 17c), and both fine and coarse modes of the size distribution had large values (Fig. 17d). Thus, we interpreted the aerosols in layer 2 as a mixture of the dust and pollution aerosols. The single-scattering albedo was largest in layer 2, and it was mostly constant at about 0.97 above the boundary layer (Fig. 17b); smaller values of the single-scattering albedo, ranging from 0.7 to 0.95, were observed in the boundary layer. Uchiyama et al. (2014) reported that the single-scattering albedo at 550 nm near the surface at Tsukuba is typically from 0.7 to 0.95, based on scattering and absorption coefficients measured by nephelometer and particle soot/absorption photometer. Thus, our estimates of the single-scattering albedo range in the boundary layer were typical.

The solar heating rate calculated from the estimated aerosol optical properties is displayed in Fig. 17e. Although the extinction coefficient was largest in layer 2, the solar heating rate was largest in layer 3. To investigate how the vertical profiles of the single-scattering albedo and asymmetry factor influenced the solar heating rate, we calculated the solar heating rate from the vertical means of the single-scattering albedo and asymmetry factor and then calculated the difference between the solar heating rate calculated from estimated optical properties and that calculated from the vertical means. The difference in the case of the solar heating rate calculated from the vertical mean of the asymmetry factor showed that the solar heating rate decreased in layer 3 and increased in layer 2 (Fig. 17f). Theoretically, a large asymmetry factor increases the downward solar flux and the solar heating rate. Because the asymmetry factor of the pure dust in layer 3 was very large, the solar heating rate was large throughout layer 3 (Fig. 17e), whereas the impact of the single-scattering albedo on the vertical profile of the solar heating rate was small (not shown). Therefore, the vertical profile of the asymmetry factor played an important role in creating vertical variation in the solar heating rate.

We next estimated the aerosol optical properties and the solar heating rate in the case of transported smoke in 8 May 2013 (Fig. 19). A high aerosol load was observed at altitudes from 4 to 6 km (Fig. 19a). In this layer, the fine mode dominated the size distribution (Fig. 19d), and the asymmetry factor was about 0.62 (Fig. 19c). The backward trajectory (Fig. 20) indicated that this aerosol had been transported from the region southeast of Lake Baikal, where forest fires had been observed in early May 2013 (Fig. 20). We interpreted the high volume of fine mode particles as transported smoke from the forest fires. The single-scattering albedo of the transported smoke was estimated to be about 0.96, which is larger than the typical value for smoke (0.88 to 0.94 at 440 and 0.84 to 0.94 at 670 nm; Dubovik et al., 2002). This overestimation may be attributed to the inability of our algorithm to reproduce the vertical profiles of the single-scattering albedo when the fine mode is dominant, as shown by our sensitivity tests (Sect. 3). In the boundary layer, the single-scattering albedo was low and the asymmetry factor was high (Fig. 19b and c). These values may reflect locally emitted soil particles; the land surface at Tsukuba was dry during 10 days before 8 May 2013 (total rainfall was only 0.5 mm), so in many agriculture and urban development areas the ground was bare.

The solar heating rate (Fig. 19e) was consistent with the vertical profile of the extinction coefficients. The difference between the solar heating rate calculated using estimated optical properties and that calculated using the vertical mean of the asymmetry factor (Fig. 19e) showed that the solar heating rate increased in the smoke layer, because the asymmetry factor of the smoke aerosol was smaller than the vertical mean value (Fig. 19c). The influence of the single-scattering albedo on the vertical profiles of the solar heating rate was small (not shown).

The SKYLIDAR algorithm showed the detailed vertical structures for the transported dust and smoke and the relationship of the aerosol vertical profiles to the solar heating rate. Our results suggest that the vertical variation of the asymmetry factor plays an important role in creating vertical variation in the solar heating rate.

Comparisons of columnar properties and surface solar irradiance

To validate the estimated vertical profiles of the aerosol optical properties, direct measurements by airplane or balloon are necessary, but such measurements are not obtained easily. Therefore, we compared the columnar optical properties of the SKYLIDAR retrievals with SKYRAD.PACK retrievals (version 4.2), which have been evaluated by Che et al. (2008) and Estellés et al. (2012). Because the direct observation of the solar heating rate is also difficult, we compared the surface solar irradiance calculated from the SKYLIDAR retrievals with that measured by pyranometer.

We compared the aerosol optical thickness, the single-scattering albedo, the asymmetry factor, and the normalized volume size distribution in the column for 2012–2013 at Tsukuba between SKYLIDAR and SKYRAD.PACK (Fig. 21). The SKYRAD.PACK estimates the columnar values of the real and imaginary parts of the refractive index, the size distribution, the optical thickness, the single-scattering albedo, and asymmetry factor from the SKYR measurements. The most different point from step 1 of the SKYLIDAR algorithm is the assumptions of the size distribution. In the SKYRAD.PACK, only the spherical particle is assumed, and the size distribution consists of the 20 discrete bins. In the SKYLIDAR algorithm, the non-spherical particle is considered, and the bi-modal size distribution is used.

The aerosol optical thickness at 532 and 1064 nm in the SKYLIDAR retrievals agreed well overall with those of the SKYRAD.PACK retrievals. Although slightly underestimated, the SKYLIDAR single-scattering albedo at 532 nm agreed well with the SKYRAD.PACK retrieval for aerosol optical thickness of more than 0.2. The single-scattering albedo estimated from the SKYR measurements by SKYRAD.PACK, however, is larger than that of the AERONET retrievals (Che et al., 2008). Thus, the SKYLIDAR results may be close to the AERONET retrievals. Similarly, the asymmetry factor also agreed with those estimated by SKYRAD.PACK for aerosol optical thickness of more than 0.2. Comparison of the 2-year mean of the normalized volume size distribution showed agreement with respect to the fine mode but not the coarse mode; the difference in the coarse mode was due to the assumption of a bi-modal size distribution by our algorithm. In the SKYRAD.PACK retrieval, a second coarse mode was observed at a radius of 10 $\mathrm{\mu }$m. This second coarse mode was reported by Che et al. (2008) and Estellés et al. (2012), who indicated that it is not observed in AERONET retrievals; therefore, it may be attributed to the difference in the retrieval algorithm. The assumption of a bi-modal size distribution is not ideal, but it is a simple way to prevent the occurrence of an unrealistic second coarse mode.

Surface solar irradiance calculated from the SKYLIDAR retrievals, and also that calculated from the SKYRAD.PACK retrievals, was compared with that measured by the pyranometer during the 2 years (Fig. 22). The surface irradiances calculated from both the SKYLIDAR and SKYRAD.PACK retrievals agreed very well with the measurements, with a very small error of around 10 W m${}^{-\mathrm{2}}$ (about 1.7 %). The SKYLIDAR result was slightly better than SKYRAD.PACK result. The remaining part of the error is attributed to the error in the input water vapor content and the assumed optical properties at wavelengths less than 532 nm and greater than 1064 nm (described in Sect. 2.2). In this study, we used the water vapor content observed by a sonde launched at 09:00 JST near our observation site.

The estimated vertical profiles of the aerosol optical properties and the solar heating rate were not validated against direct measurements in this study, but the columnar properties of the SKYLIDAR retrievals agreed well with the SKYRAD.PACK retrievals, and the surface solar irradiance calculated from the SKYLIDAR retrievals were sufficiently accurate to explain the measured surface solar irradiance.

Summary

We developed the SKYLIDAR algorithm for estimating the vertical profiles of aerosol optical properties from the sun photometer SKYR and the lidar MIEL (and HSRL) measurements. The algorithm consists of two retrieval steps. The columnar properties are first estimated from the SKYR measurements and the vertically mean depolarization ratio obtained from the MIEL measurements. Then, the vertical profiles are estimated from the MIEL (and HSRL) measurements and the columnar properties determined in the first step. The finally derived parameters are the vertical profiles of the size distribution, refractive index, volume ratio of the non-spherical particles to total particles in the coarse mode, extinction coefficient, single-scattering albedo, and asymmetry factor. In addition, we estimated the vertical profile of the solar heating rate from the SKYLIDAR retrievals.

To evaluate the performance of the algorithm, we conducted the sensitivity tests using simulated SKYR and MIEL (and HSRL) data of the vertical profiles of three different aerosols. The vertical profiles of the size distribution, the extinction coefficient, and the asymmetry factor were well estimated in all the tests. The refractive index and the single-scattering albedo in the case of dust were well estimated but not in the cases of pollution aerosol. This is because the extinction-to-backscatter ratio is sensitive to the refractive index (both the real and imaginary parts) in the case that the coarse mode is dominant but is not in the case of that the fine mode is dominant.

We then applied the SKYLIDAR algorithm to actual SKYR and MIEL measurements obtained in 2012 and 2013 in Tsukuba, Japan. Our algorithm showed the detailed vertical structures for transported dust and smoke. In addition, the vertical profiles of the solar heating rate were estimated from the SKYLIDAR retrievals, and the relationship of the aerosol optical properties to the solar heating rate was investigated. The results suggest that the vertical profile of the asymmetry factor plays an important role in creating vertical variations of the solar heating rate.

To evaluate the validity of the SKYLIDAR retrievals and the solar heating rate, we compared the columnar properties between SKYLIDAR and SKYRAD.PACK retrievals during 2012 and 2013, and we compared the surface solar irradiances calculated from the SKYLIDAR retrievals with those measured by pyranometer. The columnar properties of the SKYLIDAR retrievals agreed well with those of the SKYRAD.PACK retrievals when the aerosol optical thickness at 532 nm was more than 0.2. The calculated surface solar irradiances also agreed well with the pyranometer measurements; the mean error was only 1.7 %, despite the assumption that extended the aerosol optical properties at 532 and 1064 nm to the broadband wavelength regions. The columnar properties of the SKYLIDAR retrievals agreed with those estimated by the widely used method of SKYRAD.PACK, and the SKYLIDAR retrievals were sufficiently accurate to evaluate surface solar irradiance.

We have a plan to analyze the long-term measurements of SKYR (SKYNET) and MIEL (AD-Net) in the East Asian region using the SKYLIDAR algorithm. The results of the aerosol vertical profile and the heating rate would reveal the characteristics of the locally emitted and transported aerosols, and their influences to the temperature profiles. We focused on the measurements of SKYNET and AD-Net, but the SKYLIDAR algorithm can be applied to an another data set similar to SKYR and MIEL measurements. This flexibility is expected to be useful for investigating the temporal and spatial distribution of aerosols at different observational sites. In addition, the minimization procedure with our developed logarithm transformations, which works well for hundreds of estimated parameters, is useful for the various remote sensing. Using this minimization techniques, we are now developing the synergetic methods to estimate the aerosol vertical profile from the space-borne lidar and imager data, a combination of CALIOP and MODIS, and a combination of ATLID and MSI/EarthCARE (Illingworth et al., 2015).

Comparisons of the surface solar irradiance calculated from the SKYLIDAR (left) and SKYRAD.PACK (right) retrievals with the pyranometer measurements during 2012 and 2013 in Tsukuba, Japan.

[Figure omitted. See PDF]

Data availability

The lidar data are available from the AD-Net (http://www-lidar.nies.go.jp/). The sky radiometer data are available from the SKYNET (http://www-lidar.nies.go.jp/skynet/), but the sky radiometer data in Tsukuba, Japan, are available on request by contacting the first author of the paper.

Acknowledgements

This work was supported by the Japan Society for the Promotion of Science KAKENHI grant nos. 24510026, 15H01728, and 15H02808. The authors are grateful to the OpenCLASTR project for allowing us to use the SKYRAD.PACK (sky radiometer analysis package) in this research. NCEP reanalysis data were provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from its website at http://www.esrl.noaa.gov/psd/. The MODIS MCD12C1 product was retrieved from the online data pool, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota (https://lpdaac.usgs.gov/data_access/data_pool). Edited by: P. Stammes Reviewed by: two anonymous referees