Capabilities of an Automatic Lidar Ceilometer to Retrieve Aerosol Characteristics within the Planetary Boundary Layer

: Continuous observation and quantitative retrieval of aerosol backscatter coefﬁcients are important in the study of air quality and climate in metropolitan areas such as New York City. are ideal for this application, but aerosol backscatter coefﬁcient retrievals from ceilometers are challenging and require proper calibration. In this study, we calibrate the ceilometer (Lufft CHM15k, 1064 nm) system constant with the molecular backscatter coefﬁcient and evaluate the calibrated proﬁles with other independent methods, including the water-phase cloud method and comparison with the NASA Cloud-Aerosol Lidar and Infrared Pathﬁnder Satellite Observation (CALIPSO) attenuated backscatter coefﬁcient proﬁle. Multiple-day calibration results show a stable system constant with a relative uncertainty of about 7%. We also evaluate the overlap correction for the CHM15k ceilometer (provided by Lufft) with a Vaisala CL-31 ceilometer, and the results show good consistency between two ceilometers’ range-corrected signal (RCS) proﬁles above 200 m. Next, we implement a forward iterative method to retrieve aerosol backscatter coefﬁcients from continuous ceilometer measurements. In the retrieval, the lidar ratio is constrained by the co-located NASA AERONET radiometer aerosol optical depth (AOD) retrieval and agrees with the AERONET lidar-ratio products, derived from aerosol microphysical parameters. The aerosol backscatter coefﬁcient retrievals are validated with co-located elastic-Raman lidar retrievals and indicate a good correlation ( R 2 ≥ 0.95) in the planetary boundary layer (PBL). Furthermore, a case study shows that the ceilometer retrieved aerosol extinction coefﬁcient proﬁles can be used to estimate the AOD of the PBL and the aloft plumes. Finally, simulations of the uncertainty of aerosol backscatter coefﬁcient retrieval show that a calibration error of 10% results in 10–20% of relative error in the aerosol backscatter coefﬁcient retrievals, while relative error caused by a lidar-ratio error of 10% is less than 4% in the PBL.


Introduction
Ceilometers have been widely used to observe cloud base height and mixing layer height (MLH) and are deployed as networks due to unattended operation and easy maintenance [1,2]. Many efforts have been made to optimize or evaluate the algorithm for the MLH estimation [3][4][5][6][7][8][9]. The MLH and cloud base product can be potentially assimilated into the weather research and air quality model [10]. Importantly, continuous ceilometer monitoring can provide range-resolved aerosol distribution in the planetary boundary layer (PBL) and lower troposphere. Vertical distribution and mixing of aerosols and other pollutants in the PBL are critical to characterize air pollution formation and dispersion processes. At the same time, aerosol vertical distribution or range-resolved information is also important to better connect satellite-measured aerosol optical depth (AOD) to the near-surface PM2.5 (particulate matter with diameter equal or less than 2.5 µm) [11,12]. To quantitatively obtain aerosol backscatter coefficients from ceilometers, measured backscatter signal returns need to be calibrated. Since the ceilometer is mainly used to measure cloud base heights and aerosol layers using relative backscatter signals, the absolute calibration is not yet a standard output [3,13,14]. Alternatively, an on-site calibration with the molecular Rayleigh method or the cloud method can be used [15][16][17]. However, independent evaluation or cross validation of the calibration is highly desirable as each method includes a set of assumptions and approximations. For instance, the Rayleigh method assumes that the selected calibration altitude range is free of aerosols while the water-phase cloud method assumes an a priori lidar ratio (S c ) for the cloud, which is defined as the extinction coefficient to backscatter coefficient ratio [18].
The calibrated backscatter coefficient, or attenuated backscatter coefficient, provides the semi-quantitative information for the aerosol loadings as it includes total backscatter coefficient and atmospheric attenuation from both aerosols and molecules. Near the ground level, the ceilometer-measured attenuated backscatter coefficient can be a reasonable proxy to the ground PM2.5 concentration, since the total attenuation is negligible and good temporal correlations between the near-surface attenuated backscatter coefficient and ground PM2.5 concentration have been demonstrated [19,20]. However, when the laser beams propagate a substantial range, the atmospheric transmittance will influence the attenuated backscatter coefficient more significantly. On the other hand, when the aerosol concentration is low, the molecular backscatter coefficient will contribute more to the total attenuated backscatter coefficient. In addition, to obtain the aerosol optical depth, aerosol extinction coefficient profiles are required. Therefore, quantitative retrieval of aerosol backscatter coefficient is necessary to estimate aerosol concentration in the atmosphere.
Conventionally, elastic-scattering lidar systems are used to derive aerosol backscatter or extinction coefficient profiles with a well-established inversion approach: the Fernald-Klett method [21,22]. In this approach, when retrieving each profile, the lidar system constant (C L ) is replaced by a boundary value at a reference height in the upper troposphere where the air is assumed to be aerosol-free or have an assigned small reference value of aerosols. Then, the aerosol backscatter coefficient is retrieved from the boundary value downwards. However, for ceilometers which are eye-safe and low power, the small SNR of the profile at the aerosol-free layer in the upper troposphere can introduce significant uncertainties in the aerosol backscatter coefficient retrieval when using the traditional Fernald-Klett backward method, particularly for the data in the daytime or when low clouds are present [14]. In addition, the a priori lidar ratio of aerosols has to be assumed to solve the ill-posed equation for the single wavelength elastic-scattering lidar.
Some previous efforts were made for retrieving aerosol backscatter coefficient from ceilometer measurements and the results demonstrated good potential and inherent difficulties. Binietoglou et al. proposed a two-step approach based on the Fernald-Klett backward method. In the first step, ceilometer signals are integrated over a long period of time (up to ∼8 h) to improve the SNR and obtain the boundary value at the reference height. In the second step, the boundary value is used to retrieve the aerosol backscatter coefficient with a higher temporal resolution using the Fernald-Klett backward method. They presented a few cases showing the agreement between the ceilometer and Potenza EARLINET Raman lidar (PEARL) retrieved aerosol backscatter coefficients, but further investigation is needed for cases with different aerosol loading and distribution [23]. Wiegner et al. (2014) discussed both forward and backward approaches based on the Fernald-Klett method that derive aerosol backscatter coefficient from ceilometers (Vaisala and Lufft) and their calibration. Their results indicated that a promising relative error in order of 10% may be achieved on retrieving aerosol backscatter coefficients under favorable conditions. However, the water vapor absorption may contribute as much as 20% error for ceilometers emitting wavelengths at 905-910 nm [2,14,16]. On the other hand, Madonna et al. compared the ceilometer attenuated backscatter coefficient to a multiple-wavelength Raman lidar in the INTER-comparison of Aerosol and Cloud Tracking (INTERACT) campaign over a period of 6 months, and the results indicate their differences related to different SNRs and the ambient temperature effects on the stability of ceilometer calibration. Furthermore, they pointed out that the calibration of the ceilometer should be frequently checked and carefully evaluated [24]. Additionally, Jin et al. evaluated ceilometer attenuated backscatter coefficients and calibration constant by using both Rayleigh method and liquid cloud methods. Compared to their lidar calibrated attenuated backscatter coefficient, the CHM15k attenuated backscatter coefficient calibrated by the Rayleigh method is underestimated by 15%, due to the underestimation of the scattering ratio and lidar ratio, while the liquid cloud method underestimates by 9%, caused by the overestimated transmittance below the cloud and uncertainty of the multi-scattering factor [15]. Hopkin et al. showed that the liquid cloud method can be applied to different ceilometers under certain conditions [25].
Overall, the previous work reviewed above indicates the potential of ceilometers for retrieving aerosol backscatter coefficients, but periodic validation of the ceilometer calibration constant and the aerosol backscatter coefficient retrievals are critical to ensure the accuracy. For these reasons, we assess the calibration constant and aerosol backscatter coefficient retrieval of a Lufft CHM15k ceilometer by using a suite of co-located remote sensing observations, including the CCNY elastic-Raman lidar, AERONET sun/sky radiometer and Vaisala CL-31 ceilometer, as well as the NASA spaceborne lidar Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations/Cloud Aerosol Lidar With Orthogonal Polarization (CALIPSO/CALIOP). The cross-validation of the calibration constant is carried out between the Rayleigh and cloud methods and then the calibrated attenuated backscatter coefficient is compared to the CALIPSO data. A forward iterative method is proposed to derive aerosol backscatter coefficients in the daytime and below clouds. In retrieving aerosol backscatter and extinction coefficient profiles, the lidar ratio is constrained by the co-located AERONET AOD and shows reasonable agreement with the independent AERONET retrieved lidar ratio, derived from aerosol microphysical parameters. Then, the retrieved aerosol backscatter coefficients are validated with CCNY elastic-Raman lidar retrieval at the same wavelength, and they show strong correlation (R 2 ≥ 0.95) in the PBL. Furthermore, this retrieval approach is used in a case study to demonstrate that the ceilometer retrieved aerosol extinction coefficient profiles can be used to estimate the AOD of the PBL and of aloft aerosol plumes. Finally, we carry out simulations to estimate the retrieval uncertainty of aerosol backscatter coefficients due to the errors in the calibration constant and lidar ratio. These results show that the ceilometer is able to quantitatively retrieve the aerosol optical properties, which provide valuable information for air quality studies and satellite AOD product evaluation. This paper is organized as follows: Section 2 provides the instrumentation and data processing schemes. Section 3 describes the aerosol backscatter coefficient retrieval procedure using the CHM15k data, including the near-range signal evaluation, the forward iterative inversion method and ceilometer system constant calibration. Section 4 presents the results, cross-evaluation and application of CHM15k ceilometer aerosol backscatter coefficient retrievals. Section 5 discusses the aerosol backscatter coefficient retrieval accuracy influenced by the uncertainties of calibration constant and lidar ratio. Section 6 summarizes the main conclusions.

Instrumentation and Data Processing
In this study, we use a Lufft CHM15k ceilometer carrying out 24 h/7-day measurements on the roof of the Engineering School of City College of New York (CCNY, 40.82 • N, 73.95 • W, altitude: 100 m). A Vaisala CL31 ceilometer, a 3-wavelength elastic Raman lidar and a CIMEL radiometer (NASA AERONET) are also deployed with a distance of 5-10 m to the Lufft CHM15k ceilometer. The Lufft CHM15k ceilometer is equipped with a diode pumped Nd:YAG solid state laser at 1064 nm, with pulse energy of 7-9 µJ and repetition rate of 5-7 kHz. The receiver field of view is 0.45 mrad and a photon counting mode APD is used for backscatter signal detection. Given the biaxial optical design, the CHM15k ceilometer starts to receive useful signals from 120 m and its full overlap range is about 1.5 km [15,26]. The manufacturer provides an overlap function for each instrument upon request and raw data with overlap correction are accessible. The maximum range of CHM15k data is 15 km. All data are stored in NetCDF format on a daily basis with temporal/spatial resolution of 15 s/15 m. The Lufft CHM15k ceilometer provides up to three aerosol layer heights using the gradient method. The PBL-Height (PBLH) can be determined using the lowest aerosol layer heights with an additional quality control/quality assurance process and details are discussed by Li et al. [20].
The Vaisal CL31 ceilometer uses an InGaAs diode laser at a wavelength of 910 ± 10 nm. The pulse energy is 1.2 µJ ± 20% and the repetition rate is 10 kHz. Due to its lower laser power, the maximum detection range is 7.5 km. The spatial/temporal resolution of the Vaisala CL31 ceilometer is 16 s/10 m. Due to the coaxial optical design, its full overlap range is 70 m [27,28].
The CCNY elastic Raman lidar utilizes a Nd:YAG laser, emitting three wavelengths (1064, 532 and 355 nm) at a repetition rate of 30 Hz (Spectra-physics Quanta-Ray PRO-320) and the laser beam divergence is less than 0.5 mrad. The measurement range of the lidar is 0.5-15 km with temporal/spatial resolution of 1 min/3.75 m. Due to the eye safety issue and restrictions of the local airport and Federal Aviation Agency, lidar only collects data in the daytime of weekdays under appropriate weather conditions. More detailed specifications of the CCNY lidar can be seen in [29].
The AERONET CIMEL radiometer provides information to calculate the columnar aerosol optical depth (AOD) between 300 nm and 1640 nm and other aerosol optical properties using the AERONET Version 3 algorithm. There are three quality levels (Levels 1.0, 1.5, 2.0) for each product [30,31]. Because the Level 2.0 data are not available for the study period at the CCNY site, we use Level 1.5 cloud-screened AOD and lidar ratio products for this study.
The data used for aerosol backscatter coefficient retrieval are the Lufft CHM15k ceilometer range-corrected signal (RCS). The RCS below 200 m is replaced by the value at 200 m due to the poor overlap of the CHM15k ceilometer at the near range. All the RCS profiles are averaged every 10 minutes and smoothed along the vertical range with different running mean windows, based on the average observed SNR at different vertical ranges. For z ≤ 1.5 km, a running mean window of 100 m is applied; 1.5 km ≤ z ≤ 3 km: 200 m; z ≥ 3 km: 300 m.

Near Range Signal Evaluation
Prior to aerosol backscatter coefficient retrieval, overlap function correction needs to be applied to the original ceilometer RCS profiles. For the Lufft CHM15k ceilometer, however, some previous studies have found that the overlap function provided by the manufacturer may not be efficient enough to correct the RCS profile [15,24]. To evaluate the manufacturer overlap-corrected CHM15k profiles in the near range, the RCS from a ceilometer with a low full overlap range can be utilized as a reference. At the CCNY site, a Vaisala CL31 ceilometer is deployed 5 m away from the CHM15k ceilometer, carrying out simultaneous measurements. Since the full overlap range of the Vaisala CL31 ceilometer is as low as 70 m, the near-range signal of the Lufft CHM15k can be evaluated by the co-located CL31 ceilometer RCS profile. For comparison, RCS profiles with no cloud are selected and averaged over a 30 min interval. The profiles with multiple aerosol layers are excluded from the comparison because CL31 has a relatively lower SNR and thus may not be sensitive enough to detect the fine aerosol spatial variations. Figure 1a shows the normalized RCS profiles of CHM15k and CL31 ceilometers at 19:00 UTC 27 July 2020. Both profiles are averaged over 18:45-19:15 UTC and are normalized such that they have a unit value around 700 m [14]. Between 200 m and 1600 m, the profiles from two ceilometers generally match, except that the profile of CL31 shows complex "ripple" pattern noise which is discussed by Kotthaus et al. [32]. The signal from Luftt CHM15k below 200 m decreases significantly, indicating the poor overlap. The averaged ratio of the RCS profiles on 27 July 2020 is displayed in Figure 1b. The solid line represents the mean value of the 23 calculated profile ratios from 1:00-23:00 UTC (each ratio is calculated by the 30 min averaged normalized profiles at each hour), and the error bar represents the standard deviation. Figure 1c. shows a similar result to Figure 1b except that the data are selected from 10 to 15 October 2020 (profiles with low clouds and precipitations are excluded). The normalized ratio of CHM15k and CL31 RCS profiles varies within 0.9-1.1 and is relatively stable between 200 m and 700 m. From the RCS profile comparison between two ceilometers, we conclude that the CHM15k signals above 200 m have been properly corrected by the overlap function from the manufacturer but the signals below 200 m cannot be utilized. To simplify the calculation in the retrieval, we assume a uniform distribution below 200 m and replace the RCS below 200 m by its value at 200 m. Given that the molecular and aerosol transmittance below 200 m is nearly 1 for wavelength at 1064 nm [14], this will hardly affect the aerosol backscatter coefficient retrieval above 200 m, but when the near-ground aerosol backscatter coefficients are of interest, the RCS of CHM15k at 100-200 m should be corrected by the normalized ratio of CHM15k and CL31 at 100-200 m.

Aerosol Backscatter Retrieval with Forward Iterative Method
For aerosol backscatter coefficient retrieval, a detailed description of the forward iterative method is given here. Considering the single-scattering lidar equation with both aerosol and molecular components [21]: where P(z) is the returned signal at altitude z and P(z) · z 2 is the RCS. C L is the lidar (ceilometer) system constant whose value is determined by the laser pulse energy and the efficiency of the receiver. β m (z) and β a (z) are the molecular and aerosol backscatter coefficient, respectively. T(z) 2 is the two-way transmittance of the atmosphere, which can be further described by Equation (2): where α m (z ) and α a (z ) are the molecular and aerosol extinction coefficients, respectively. To solve the lidar equation, the aerosol lidar ratio S a is assumed to be constant along each profile, therefore, we have: The value of S a can be estimated based on the a priori knowledge of the aerosol type or constrained by AERONET aerosol optical depth (AOD) when applicable. Another issue is the unknown lidar system constant (C L ), which needs to be calibrated before retrieving the attenuated backscatter coefficients. The details of how to calibrate the ceilometer system constant are described in Section 3.3. After determining the lidar ratio and lidar constant, we obtain: Based on the Rayleigh scattering theory, the molecular backscatter coefficients and the two-way molecular transmittance can be calculated from the pressure and temperature profile and S m = 8π/3 at 1064 nm. In the forward iterative approach, the aerosol backscatter coefficient is retrieved from the near-surface point upwards.
As shown in Figure 2, to begin, since the aerosol transmittance below z 0 , i.e., T a (z 0 ) 2 , is almost 1 for near-infrared, the aerosol backscatter coefficient near the surface β a (z 0 ) is initialized with the attenuated aerosol backscatter coefficient P , then the aerosol extinction coefficient (α a (z 0 )), AOD (τ(z 0 )) and aerosol transmittance (T a (z 0 )) can be determined with the a priori lidar ratio (S a = 40 sr) or the constrained lidar ratio, which will be discussed later. Afterwards, β a (z 0 ) and α a (z 0 ) are calculated again with the updated aerosol transmittance (T a (z 0 )). For the second point (z 0 + δz), the initiated value of α a (z 0 + δz) is set to 0 and then τ(z 0 + δz), T a (z 0 + δz), β a (z 0 + δz) and an updated value of α a (z 0 + δz) is calculated accordingly. This process is iterated until one of the following conditions is satisfied: (1) the relative difference of the present extinction coefficient value and the extinction coefficient value from the previous iteration is less than 0.01%. (2) The iteration number is equal to 30. When the iteration is finished, the outputs τ(z 0 + δz), α a (z 0 + δz) and β a (z 0 + δz) are saved, and the program moves to the next point (z 0 + 2δz) and we follow the same procedure outlined above to calculate τ(z 0 + 2δz), α a (z 0 + 2δz) and β a (z 0 + 2δz). This process is repeated until z = h top (h top = 7.5 km or the cloud base height). To constrain the lidar ratio, we varied its value from 20 sr to 70 sr and match the ceilometer AOD, which is calculated by integrating the extinction coefficient profile from z 0 to z max , to the co-located AERONET retrieved AOD. When no low cloud or aloft aerosol layers appear, z max is set to 4.5 km, as the aerosol loading above the PBL is low, and above 4.5 km the ceilometer signal is often noise dominated (SNR < 1).

Rayleigh Calibration
When using the forward iterative method to retrieve the aerosol backscatter coefficient, it is crucial to calibrate the ceilometer system constant C L . One common approach is the Rayleigh method. To apply this method, we need to find a cloud-free vertical range where no aerosols or very few aerosols appear. Assuming that the aerosol backscatter coefficient in this range is equal to 0 and the aerosol transmittance is close to 1, then the lidar equation in this range can be approximated by Equation (5): where P(z) · z 2 is the RCS and β m (z) · T m (z) 2 is the product of molecular backscatter coefficient and its transmittance. By fitting the ceilometer RCS profile and the reference molecular backscatter coefficient profile with a linear function, the resulting linear regression slope is the ceilometer calibration constant. The reference molecular backscatter coefficient and transmittance are calculated using a radiosonde profile which is measured twice a day (00:00 UTC and 12:00 UTC) at the National Weather Service, Upton, New York (about 90 km from the CCNY site). The major difficulty of this method is that the Rayleigh scattering intensity is proportional to the reciprocal of the wavelength to the fourth power, therefore the molecular backscatter coefficient at 1064 nm is very small. In addition, although at a higher altitude, for example 6-7 km, the air can be generally identified as aerosol-free, the ceilometer SNR is poor for the RCS from that range, due to its inherent low laser power. In this study, we select the RCS profiles from 3-6 km during the nighttime without aerosol plumes or clouds nearby and we average the RCS profiles over a 3-5 h period to improve the SNR. Note that we only select the cases for calibration, in which the correlation R 2 of the ceilometer RCS and the molecular backscatter coefficient profile is greater than 0.9.  Figure 3 is an example of using the Rayleigh method to calibrate the CHM15k system constant. Figure 3a shows the 5 min averaged RCS profiles of 19 December 2018. During the nighttime, PBL-Height (PBLH) is relatively stable, remaining around 700-800 m and no cloud or aerosol plume appears within a 2-6-km range. The AERONET averaged AOD from 13:00-17:00 UTC on this day is 0.015 at 1020 nm and the aerosol transmittance is about 0.97. For calibration, the RCS profiles are further averaged over 0:00-4:00 UTC. The averaged RCS and molecular backscatter coefficients at 0:00 UTC of 3.2-5 km are shown in Figure 3b. The correlation R 2 is about 0.95 and the linear regression slope is 170.97, which represents the ceilometer system constant C L . Figure 3c shows the calibrated attenuated backscatter coefficient profile using C L = 170.97 km 3 ·sr (averaged over 0:00-4:00 UTC, blue solid line) and the molecular backscatter profile coefficient (0:00 UTC, red solid line), demonstrating a good agreement between the attenuated backscatter and molecular backscatter coefficients above 3 km. From December 2018 to February 2020, we selected 12 cases that meet the Rayleigh calibration requirements for calculating the ceilometer system constant and the detailed information is displayed in Table 1. Within these 12 studied cases, we found C L = 161.8 ± 10.78 km 3 ·sr (mean ± std), corresponding to about 7% relative uncertainty. The calibration results are relatively stable and consistent, therefore, we used C L = 162 km 3 ·sr to calibrate the attenuated backscatter coefficients of CHM15k in the retrieval.

Liquid Cloud Calibration
Another typical approach to calibrate the system constant is the liquid stratocumulus cloud method, proposed by O'Connor et al. [18]. Briefly, the path integration of the attenuated backscatter coefficient I(z b , z t ) of the cloud base z b to cloud top z t can be estimated by Equation (6) [15,33]: where T(z 0 , z b ) is the aerosol transmittance from the surface level (z 0 ) to the cloud base height (z b ). η is the multiple scattering factor, which is assumed to be 1; S c and τ c are the lidar ratio and optical thickness of the liquid cloud, respectively. Pinnick et al. [34] found that S c = 18.2 sr at 1.06 µm for a wide range of droplet diameters. When the cloud is optically thick (τ c ≥ 3) and the transmittance below the cloud is close to 1, the system constant can be approximated by: However, because CHM15k uses a photon counting system to collect the return signal, the strong return signal from the low liquid cloud may saturate the detector. In addition, the unknown multiple scattering coefficient can also introduce additional error. Therefore, the cloud sample used for this method needs to be carefully selected to exclude cloud profiles with a saturated signal [25], or a neutral density filter can be utilized to attenuate the return signal [15]. In this study, we applied the liquid cloud calibration to the cloud signals which are temporally stable and located around 1.5-2 km without precipitation. The clouds lower than this range usually result in signal saturation while higher clouds may be mixed-phase or ice clouds.
The liquid cloud calibration is applied to the RCS profiles from 18:00-19:30 UTC 18 November 2018, as shown in Figure 4 (cloud base height is around 1.5 km). The calculated system constant C L = 160 ± 8.2 km 3 ·sr (Figure 4b), which is about 6% lower than the system constant calibrated by the Rayleigh method on 19 December 2018 (as shown in Table 1), shows the consistency between the system constant calibrated by the Rayleigh method and by the liquid cloud method.

Comparison with CALIPSO Attenuated Backscatter Coefficient Profile
With the ceilometer system constant calibrated by the Rayleigh method, the attenuated backscatter coefficients of CHM15k can be obtained. We compared the attenuated backscatter coefficient profiles of CHM15k with the CALIPSO level 1 attenuated backscatter coefficient profile at 1064 nm (version 4) [35]. We only select the nighttime profiles when CALIPSO passes over near the CCNY observation site. The selected CALIPSO attenuated profile is averaged along-track over 50 km and CHM15k ceilometer data are averaged over 30 min, according to the general relationship between aerosol spatial and temporal variation scales [36,37].
The direct comparison between the attenuated backscatter coefficients of CHM15k and CALIPSO at 7:13 UTC 4 July 2019 is shown in Figure 5. The inset of Figure 5 indicates the CALIPSO ground track (red solid line) and the CCNY observation site (green square symbol). The CALIPSO attenuated backscatter coefficient profile is quite noisy, as expected, but both attenuated backscatter coefficient profiles generally show good agreement from 1 km up to 7 km, indicating that the attenuated backscatter coefficients of CHM15k have been properly calibrated. Altitude (km)   Figure 6 shows the aerosol optical property retrievals on 8 October 2020. Figure 6a is the 10 min averaged RCS; Figure 6b displays the retrieved aerosol backscatter coefficients below the clouds; Figure 6c shows the aerosol extinction coefficients; Figure 6d is the comparison of AERONET AOD (red dots) at 1020 nm and AOD from CHM15k retrieval (blue dots), which is the path integration of the aerosol extinction coefficient from the surface to 4.5 km (if there is cloud below 4.5 km, the integration upper bound is set to be the cloud base height). There are some clouds around 1.5 km at 12:00-19:00 UTC, therefore we only utilize the AERONET AOD after 19:00 UTC to constrain the lidar ratio S.

Retrieval Results
As shown in Figure 6d, after 19:00 UTC, the AERONET and ceilometer AODs match very well. The constrained lidar ratio after 19:00 UTC is S = 58.2 ± 5.6 sr. When the AERONET AOD is unavailable or unreliable due to the clouds, the mean lidar ratio S = 58.2 sr is used to retrieve the aerosol backscatter coefficients, extinction coefficients and AOD. The independently retrieved AERONET averaged lidar ratio (version 3 Aerosol Inversion Product) for the same day is 53.8 sr at 1020 nm, which is consistent with our lidar ratio estimation.   Figure 7 b,c are retrieved aerosol backscatter and extinction coefficients, demonstrating that the forward iterative method can retrieve the optical properties of the transported aerosol plumes even above the PBL. Figure 7d displays the good agreement between AERONET AOD and CHM15k retrieved AOD. The constrained lidar ratio S = 46.9 ± 4.1 sr and the averaged lidar ratio from AERONET on this day is 45.6 sr at 1020 nm. The difference between the lidar ratio of 8 October 2020 and 14 October 2020 is mainly due to the different aerosol compositions and size distributions.

Evaluation
We evaluate the aerosol backscatter coefficient retrieval results from CHM15k by comparing with the co-located CCNY lidar aerosol backscatter coefficient retrieval of the same day [29]. The Fernald-Klett backward algorithm is applied to retrieve the lidar aerosol backscatter coefficient and the lidar ratio is 40-45 sr, constrained by AERONET AOD. Figure 8a-f show the aerosol backscatter coefficients retrieved from CHM15k using the forward iterative method and lidar retrieval using the Fernald backward method on 8 and 14 October 2020. The CHM15k and lidar profiles are averaged to match each other temporally and vertically, and the temporal/spatial resolution is 10 min/15 m. On both days, the temporal variations of the CHM15k and lidar aerosol backscatter coefficient retrievals agree very well with each other. On 8 October 2020 (Figure 8a,b), the lidar aerosol backscatter coefficient retrievals below the clouds are smaller than that of the ceilometer, because when using a backward method, the aerosol backscatter coefficient retrieval accuracy below the clouds is degraded by the large uncertainty of the lidar ratio of the clouds. Figure 8c,f are the scatter plots of aerosol backscatter coefficients of the CHM15k and lidar on 8 and 14 October 2020, respectively. To eliminate the influence of the lidar non-overlap effect at the near range and the poor SNR of the ceilometer at a higher altitude, comparison data are selected from 0.7-2.5 km during the cloud-free period. Strong correlations are found on both days between the aerosol backscatter coefficient retrievals of the ceilometer and lidar, with correlation R 2 = 0.97, even though they are calculated independently with different methods. The regression slope is close to one (0.99 for 8 October in Figure 8c and 0.87 for 14 October in Figure 8f). This indicates that the aerosol backscatter coefficient retrieved by the CHM15k ceilometer is consistent with the lidar aerosol backscatter coefficient retrieval. However, a slight difference between the two retrievals can be introduced by the uncertainties of the calibration constant and the lidar ratio. 16 18 20 Time (UTC Hour)

Application
The previous cases demonstrate that the forward iterative method can be utilized to retrieve the aerosol backscatter and extinction coefficient from the CHM15k calibrated attenuated backscatter coefficient. One of its applications is that we can better quantify the aerosol loading within the PBL and of the aloft plumes. Particularly, when aloft aerosol plumes appear, we can retrieve the AOD within the PBL and the AOD of the aloft aerosol plume separately. Figure 9 gives an example of quantifying aerosol backscatter coefficients, extinction coefficients and AODs of aloft plumes that are associated with wildfire smoke transport. Figure 9a,b show the calibrated attenuated backscatter coefficient profiles and aerosol backscatter coefficient retrieval of 23 September 2020, respectively. The heavy dense aerosol plumes can be seen at 1.5-5 km at 0:00-6:00 UTC and some low and moderately dense plume appears below 3 km at 10:00-16:00 UTC. There are some low-level clouds at 20:00-22:00 UTC which totally attenuate the laser beam above the clouds. In Figure 9a,b, we compare the ceilometer attenuated backscatter coefficients and aerosol backscatter coefficient retrievals from the ceilometer and we found that: (1) the aerosol backscatter coefficients of the aloft aerosol plumes are larger than the attenuated backscatter coefficients due to atmosphere attenuation; (2) the retrieved aerosol backscatter coefficients in the lower altitude are smaller than the attenuated backscatter coefficients because of the molecular contribution included in the attenuated backscatter coefficients. Figure 9c illustrates the aerosol extinction coefficients. By integrating aerosol extinction coefficients within different ranges, the segmented PBL-AOD, aloft plume AOD and the total AOD (z max = 7.5 km) are calculated and are shown in Figure 9d. From 0:00 to 15:00 UTC, the major portion of the total AOD is contributed from the aloft aerosol plumes whereas the PBL-AOD is less than 0.02. Starting from 15:00 UTC, the PBLH increases and the aerosols from aloft plume begin to mix into the PBL; the PBL-AOD increases and AOD above the PBL decreases. After 18:00 UTC, most of the AOD is contributed by the aerosols within the PBL. Overall, quantification of the PBL-AOD, AOD of aloft aerosol plumes and columnar AOD is critical to evaluate the impacts of transported aerosol plumes on air quality and it is important for the application of satellite retrieved AOD to air quality studies.

Discussion
As mentioned before, the backscatter coefficient retrieval accuracy depends on the accuracy of different parameters. Herein, we address the relative error of aerosol backscatter coefficient retrieval contributed by the ceilometer system constant and lidar ratio, when using the forward iterative method. We assume an idealized aerosol backscatter coefficient profile, a molecular backscatter coefficient profile and an ideal lidar system with the system constant C L = 3000 km 3 ·sr and lidar ratio S = 40 sr. Then, we obtain the return signal and RCS profile of this lidar using Equation (4). We apply the forward iterative method to the RCS profiles and retrieve the aerosol backscatter coefficients using different sets of parameters (C L , S) and calculate the relative error using Equation (8): First, to analyze the influence of the uncertainty of the lidar system constant on the retrieval accuracy, we set the lidar ratio S = 40 sr and C L to be 2700, 3000, 3150 and 3300 km 3 ·sr and apply the forward iterative method to the simulated RCS profile. In Figure 10a, the magenta solid line shows the idealized aerosol backscatter coefficients β a,true with the most aerosols concentrated in the PBL and aloft aerosol layers around 4-5.5 km; the green solid line is the molecular backscatter coefficient profile and the dashed lines represent the aerosol backscatter coefficient profile retrieved by the forward iterative method using different C L . Figure 10b displays the scattering ratio R, which is defined as the ratio of the total backscatter coefficient to the molecular backscatter coefficient (R = 1 + β a /β m ), and the relative errors of the aerosol backscatter coefficient retrievals with different C L . Within the PBL, where aerosol concentration is larger, indicated by a greater R value, increasing C L by 10% (C L = 3300 km 3 ·sr) will introduce a relative error of ∼11% (purple solid line in Figure 10b) and decreasing by 10% (C L = 2700 km 3 ·sr) will cause a relative error of ∼14% in the aerosol backscatter coefficient retrieval (orange solid line in Figure 10b). Above the PBL where the aerosol concentration is very low, corresponding to the smaller R value, a 10% change can introduce a relative error of more than 35% in the aerosol backscatter coefficient retrieval. For the aloft aerosol layer around 4-5.5 km, increasing by 10%(C L = 3300 km 3 ·sr) will introduce a relative error of ∼18% and decreasing by 10% (C L = 2700 km 3 ·sr) will cause a relative error of ∼24% in the aerosol backscatter coefficient retrieval. Next, the relative error introduced by the uncertainty of the lidar ratio is investigated. We set C L = 3000 km 3 ·sr and vary the lidar ratio S (S = 30, 40 and 50 sr). Figure 10c demonstrates the aerosol backscatter coefficient retrievals with different S and Figure 10d shows the corresponding relative errors. In general, compared to C L , the lidar ratio has less influence on the retrieval accuracy. Increasing the lidar ratio by 25% (S = 50 sr) only introduces a relative error of less than 10% into the aerosol backscatter coefficient retrieval. Figure 11a,b show the relative error of aerosol backscatter coefficient retrieval when we increase by 10% on 8 and 14 October 2020. The relative error is calculated using Equation (8). The blue color indicates that the relative error is below 20% and cyan color means the relative error equals 20-30%. Figure 11c,d display the corresponding scattering ratio for these two days. The region with a warmer color (cyan to red) indicates that the retrieved scattering ratio is larger than 2, that is, the retrieved aerosol backscatter coefficient is larger than the molecular backscatter coefficient. The studied cases are consistent with the simulation results and we conclude that: (1) within the PBL, where the aerosol backscatter coefficient is larger than the molecular backscatter coefficient, an increase of 10% of the system constant will introduce a relative error of 10-20% into the aerosol backscatter coefficient retrieval. (2) A 10% relative uncertainty of the S ratio will only cause a relative error of less than 4% in the aerosol backscatter coefficient within the PBL (not shown).

Conclusions
In this study, we present a forward iterative method for aerosol backscatter coefficient retrieval applied to a Lufft CHM15k ceilometer which can retrieve aerosol backscatter coefficients during the daytime and below the clouds.
The ceilometer system constant (C L ) is calibrated by the Rayleigh method: comparing the ceilometer RCS profiles to the molecular backscatter coefficients at the altitude of 3-6 km (without cloud or aerosol layers). Within 12 studied cases (2018-2020), the calibrated system constant is relatively stable, with ∼7% relative uncertainty. The system constant calibrated by the Rayleigh method is consistent with the result of using a liquid cloud method and the calibrated attenuated backscatter coefficient profile of CHM15k is comparable to the CALIPSO co-located attenuated backscatter coefficient profile at 1064 nm.
For aerosol backscatter coefficient retrieval, the lidar ratio is constrained by the colocated AERONET retrieved AOD. Two cases of CHM15k aerosol backscatter, extinction coefficients and AOD retrieved by the forward iterative method are presented. The retrieved CHM15k aerosol backscatter coefficients are validated with co-located elastic-Raman lidar aerosol backscatter coefficient retrieval at 1064 nm using the Fernald-Klett backward method. The correlation R 2 is greater than 0.95 for both studied cases and the linear regression slopes are close to 1, meaning that the forward iterative method can be applied to retrieve aerosol backscatter coefficients from CHM15k measurements during the daytime with and without the clouds. We also demonstrated that with the retrieved aerosol extinction coefficients, the AOD of PBL and of the aloft plumes can be properly estimated. Finally, we quantitatively analyzed the relative error of aerosol backscatter coefficient retrieval contributed by different parameters. Both simulation and case studies illustrate that the relative error of aerosol backscatter coefficient retrieval using the forward iterative method can mainly be attributed to the uncertainty of the calibration constant. A 10% uncertainty in the calibration constant will cause 10-20% relative error in the aerosol backscatter coefficient within the PBL, while 10% relative uncertainty of the S ratio will only cause a relative error of less than 4% in the aerosol backscatter coefficient within the PBL.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.