STUDY OF PLANETARY BOUNDARY LAYER AND AEROSOL TRANSPORT USING CEILOMETER AND AIR QUALITY MODELLING IN NEW SOUTH WALES (NSW), AUSTRALIA

: The planetary boundary layer height (PBLH) is one of the key factors in influencing the dispersion of the air pollutants in the troposphere and hence the air pollutant concentration on ground level. For this reason, accurate air pollutant concentration depends on the performance of PBLH prediction. Recently, ceilometer, a lidar instrument to measure cloud base height, has been used by atmospheric scientists and air pollution control authorities to determine the mixing level height (MLH) in improving forecasting and understanding the evolution of aerosol layers above ground at a site. In this study, ceilometer data at an urban (Lidcombe) and a rural (Merriwa) location in the New South Wales, Australia was used to validate the PBLH prediction from two air quality models (CCAM-CTM and WRF-CMAQ) as well as to understand the aerosol transport from sources to receptor point at Merriwa for the three case studies where high PM 10 concentration was detected in each of the three days. The results show that the PBLH prediction by the two air quality models corresponds reasonably well with observed ceilometer data and the cause and source of high PM 10 concentration at Merriwa can be found by using ceilometer MLH data to corroborate with back trajectory analysis of transport of aerosols to the receptor point at Merriwa. Of the three case studies, one had aerosol source from north and north west of Merriwa in remote NSW where windblown dust is the main source, and the other two had sources from south and south east of Merriwa where anthropogenic sources dominate,


Introduction
The planetary boundary layer (PBL) is considered as the lowest layer of the troposphere directly influenced by the surface forcing, such as heat transfer, frictional drag, topography, and others. The boundary layer structure and its diurnal evolution control the dispersion and the resulting concentrations of pollutants. As the planetary boundary layer height (PBLH) evolves during daytime, primary pollutants from emission sources are diluted within a larger volume of air, leading to cleaner air when photochemical production and the advection of polluted plumes have minor contributions. In contrast, after sunset, air quality may deteriorate with the existence of with strong emissions of pollutants. Monitoring of changes in the PBL height with high spatial and temporal resolution are desirable to improve air quality assessment and forecasting [1]. The diurnal evolution of the PBL is complex and typically consists of the convective mixing layer (ML) during the daytime and the residual layer (RL) during nighttime composing the remains of the daytime ML above the near-surface nocturnal stable layer (SL) [2] [3]. Technically, the daytime ML or convective boundary layer (CBL), determined usually by aerosol mixing layer profile analysis, is a subset of the PBL defined and determined by temperature profile inversion. Shi et al. 2020 [4] distinguished 4 different types of PBLH definition. One is based on temperature profile, one on material or aerosol profile, one on turbulent kinetic energy and one on wind shear. The PBLH based on aerosol gradient profile such as from lidar measurement is less that that based on temperature gradient profile. This is also confirmed in a study by Knepp et al. 2017 [5] on the assessment of ML height estimation from single-wavelength CL51 ceilometer profiles in Colorado by comparing to radiosonde data. They showed that sonde-derived boundary layer heights are higher (10-15 % at midday) than CL51 lidar-derived mixed-layer heights. There is an inconsistent usage of the name MLH, PBLH or ABLH (atmospheric boundary layer height) in literature. In this study, the term MLH indicates boundary layer height as measured by lidar such as Vaisala CL51 while PBLH as layer height predicted from the models or measured with temperature profiles such as from radiosonde.
Ground-based lidar instruments such as Micro-Pulse Lidar (MPL) systems or ceilometers have been used to profile the atmosphere and to determine the boundary layer. Recently, ceilometers are recognised as efficient and affordable ground-based instruments for profiling the atmosphere for cloud and aerosols layer observations. As the name implied, the ceilometer was used originally to determine the cloud base using lidar. Cloud base heights as detected by ceilometer was recently compared with satellite cloud type with good accuracy [6]. Since 2008, aerosol layers were included in the investigation of backscattering lidar as ceilometer reflects every particle including, rain droplets, fogs, moisture droplets, aerosols… A ceilometer measures the optical backscatter intensity in the air that depends on the particle concentration in the air. From the backscatter coefficient profile, the PBL can be identified as above PBL there is no backscattered signal. The Vaisala CL31 ceilometer is a popular instrument and is used in many parts of the world.
The PBLH detection is further enhanced in range with the new Vaisala ceilometer CL51. Both the CL31 and CL51 ceilometers use the Münkel and Roininen algorithm [7] [8] to detect the MLH with confidence levels and error bars added [9]. This proprietary algorithm is implemented in BL-View software provided with the instrument to derive the atmospheric MLH and other aerosol and cloud layer information.
These instruments have also been used to detect aerosols transport and dispersion in the boundary layer. Yang et al. 2020 [10] have used Leosphere WindCube Scan 200S Doppler lidars and Vaisala CL31, CL51 ceilometers to detect and studied dust events over Iceland from volcanic ash eruptions. They found that these instruments provide accurate monitoring of the vertical distribution and temporal evolution of aerosols in Iceland. Shang et al. 2021 [11] has used CL51 ceilometer in Kuopio (Finland) to detect aerosols in the lower troposphere, at 2 to 5 km height on 4 to 6 June 2019. These aerosols were longrange transported from biomass burnings in Canada. Illingworth et al. 2019 [1] has reported on recent developments in Europe on the exploitation of existing ground-based profiling instruments such as ceilometers, MPLs to network them together and to be able to send their real time air pollutant data to forecast centres.
Ceilometer backscatter coefficients (BSCs) data can be used in conjunction with AOD data from sun photometer to estimate the profiled aerosol mass concentration as Shang et al. 2021 have performed and compared with MERRA-2 aerosol concentration.
Currently, the Climate and Atmospheric Science (CAS) Branch of the NSW DPIE (New South Wales Department of Planning, Industry & Environment), Australia operates two Vaisala CL51 ceilometers located at Merriwa and Lidcombe sites. The Merriwa ceilometer near the Upper Hunter (northwest of Sydney) is used to detect the aerosol layers including dust transport from western NSW ( Figure 1). The data collected from the two ceilometers were also used to assess and determine their suitability for calculating the mixing layer height time series continuously within the atmospheric planetary boundary layer (PBL). Vaisala has used the proprietary algorithm based on Münkel and Roininen gradient method in the CL31 and CL51 ceilometers to derive the mixing height time series from the aerosol layer profiles as measured by the ceilometers.
Accurate and continuous measurements of MLH are needed for relevant air quality assessments. Ceilometer technology can be used to retrieve the mixing height and to produce information on the vertical mixing and atmospheric structure over a selected area as well as detecting cloud base and aerosol layers above ground. Results can be used to investigate the relationship between the evolution of the daytime mixed-layer height and air pollution under conditions where changes in concentration depend only on the urban boundary-layer growth and air entrainment from the free atmosphere.
In this study, we used the derived MLH data from the two selected ceilometers to compare with the predicted PBLH from the NSW DPIE's forecast and air quality models (CCAM-CTM and WRF-CMAQ). The data collected at the 2 selected ceilometer sites for the 2 periods in February 2021 and April 2021, were used to derive the correspondent MLH. In previous studies, Uzan et al. 2020 [12] have used MLH observed from ceilometers in Israel to validate two meteorological models, the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS) global model and the COnsortium for Small-scale MOdeling (COSMO) regional model. Both of these models calculate the PBL height using the bulk Richardson number (Rib) method. They found that the results are reasonably good when a correction tool was applied to PBLH prediction at sites in complex terrain.
The validation of the NSW DPIE air quality models is important in having confidence in the forecasting ability and in scenario development of air quality model simulation for policy work. The PBLH is an important parameter in predicting air quality concentration and as such comparison with MLH lidar measurements data from the ceilometer is performed. We also use the aerosol profile and MLH as measured by the CL51 ceilometer to study the aerosol transport from sources to the receptor site at Merriwa in some case studies when high PM10 concentration occurred. Merriwa is a remote site in the Upper Hunter where dust sources from western NSW can cause high dust events in the Hunter and Greater Metropolitan Region (GMR) of Sydney.

Ceilometer Vaisala CL51 and monitoring data
The operation principle of the Vaisala CL 51 ceilometer is simple. Pulses of laser at 910 nm wavelength are emitted from transmitter into the atmosphere. The scattered signals from scatterers (particles) are then received by receiver. Signal is then range corrected (level 1) as the signal get weaker as distance larger. Range corrected signal (RCS) is obtained by multiplying the signal with r 2 where r is the distance. Data is processed to give profile data (level 2). The RCS is then used with gradient minima method to detect particle concentration boundaries. This is based on the observation that as layered particles are confined in the boundary layer, above which there will be no lidar reflection. Boundary layer is then determined from gradient method to find local minima and the lowest of these gradient minima marks the top of the mixed layer.
Various gradient methods have been proposed [13] [14] [15]. The method proposed by [7] is used in Vaisala CL31 and CL51 to produce MLH data (level 3). Validation of the method by comparing MLH as derived from ceilometer measurements with observed radiosonde data was performed at Sterling, VA, USA [7], Colorodo, USA [5], Shanghai, China [15] or with AMDAR (Aircraft Meteorological Data Relay) aviation data at Perth, WA, Australia [16].
The files from these processed levels are stored in netCDF format. Most of the time, one needs to look only at the Level 3 data files which can be used with BL-View software as provided with the instrument to determine the MLH using built-in algorithm. The detection method works best on clear sky day. During rain period where lots of echo were received, even with filtering, the method struggles or does not work. The boundary layer scan resolution of CL51 is 10m and reporting interval of 36s and the range for boundary layer fine structuring profile is 4000m. The cloud base range is 13000m as compared to 7500m of the previous Vaisala CL31.
Recent non-propriety algorithm CABAM (Characterising the Atmospheric Boundary layer based on ALC Measurements) has been used by [17] [18] and validated with AMDAR (Aircraft Meteorological Data Relay) measurements and applied to characterise the urban boundary layer over London. This is the first non-proprietary mixed layer height algorithm, specifically designed for the commonly deployed Vaisala CL31 ceilometer [17] [18]. There are other public domain algorithms to determine the cloud height and MLH from lidar measurements. The STRucture of the ATmosphere (STRAT v1.04) algorithm can determine cloud height and MLH from a variety of lidar instrument including the CL51. STRAT package is available as MATLAB and Python code and is a good alternative to BL-View from Vaisala [5] [19]. The other public available algorithm is the UMBC algorithm which was developed independently for estimating MLHs from lidar backscatter profiles using a covariance wavelet technique (CWT) similar to STRAT [5].
The MLH, as derived by Münkel and Roininen method [7] and provided by Vaisala CL51, is extracted from Level 3 used in this study. Meteorological and air quality data at the ceilometer sites are obtained from the DPIE air quality monitoring network which covers most of New South Wales. Figure 1a shows the location of the two Vaisala CL51 ceilometers in New South Wales. Merriwa is located in the Upper Hunter Region of NSW and known as one of the most fertile farming areas in the country. Dust pollution has been recognised as the major types of pollution occurring in the area. An air quality monitoring station located on the Merriwa Scone Road has been established in 2012, as a background site for fine particles PM10 measurements. The geographic coordinate of the station is -32.12 o latitude and 150.46 o longitude.

Study area
Lidcombe ceilometer is co-located at DPIE Lidcombe air quality monitoring station site situated in centre of metropolitan areas of Sydney (-33.88 o latitude, 151.05 o longitude). The location is ideal to observe the aerosol layers and PBL above the urban areas and the other ground-based meteorological and air quality variables such as wind speed, temperature, Ozone, PM10, PM2.5 and other pollutants.
Topography and meteorological conditions influence the evolution of mixing height profiles and the associated atmospheric stability classes. Selected available meteorological parameters such as wind speed and wind direction, surface temperature, and humidity were selected to study their effects on mixing height at Merriwa and Lidcombe. The time series of selected parameters were plotted against the corresponding mixing heights and the obtained results are presented in following section.
Currently DPIE uses two air quality models for forecasting daily air pollution and policy scenario development: the CCAM-CTM (Conformal Cubic Atmospheric Model -Chemical Transport Model) model developed by CSIRO [20] [21] [22] and the State-ofthe-Science WRF-CMAQ model developed by US-EPA. We use the derived PBL profiles from the ceilometer measurements at Merriwa and Lidcombe sites to validate the predictions from each of the air quality models above. The simulation domain configuration as shown in Figure 1b for the WRF-CMAQ run is a three-nested domain with the outer domain (d01) covering much of Eastern Australia at the resolution of 12 × 12 km. The inner domain (d02) is at 4 × 4 km resolution covering most of NSW, while the innermost domain (d03) is at 1 × 1 km resolution and covers the Greater Metropolitan Region (GMR) of Sydney [23]. The initial and boundary conditions and host data for meteorology is from the US National Centre for Environmental Predic- The 2 periods of validation are 12 to 19 February 2021 and 20 April 2021 to 2 May 2021. Prediction of PBLH from the models are compared with CL51 lidar measurements for those periods. In CCAM, the bulk Richardson number (Rib) method is used to estimate the PBLH while for WRF, the PBLH was estimated by using the Mellor-Yamada-Janjic local scheme (MYJ) [25][26] [27] as a configuration option in our study. Two most often used PBL schemes in WRF are the non-local first order closure YSU (Yonsei University) scheme based on the bulk Richardson number and the local second order closure MYJ scheme based on solution of the turbulent kinetic energy (TKE) second-moment budget equation. Tyagi et al. 2018 [28] in their study of comparison of different PBL schemes in WRF with aircraft observation data have shown that the local second order closure TKE schemes perform better than the first-order closure schemes such as YSU scheme [28].
The bulk Richardson number scheme is based on the simple ratio of convective turbulence (temperature) and mechanical turbulence (wind shear) and is represented by the bulk Richardson number formula.
where g is gravitational acceleration, θ is the potential temperature, z is at height z and s is first level height.
The PBLH is defined as the first height when the bulk Richardson number is greater than the critical value of ¼. The Rib is a simplified empirical form of Richardson number which is expressed as The MYJ parameterization scheme, based on the turbulence kinetic energy (TKE) budget equation, determined the boundary layer height as where the TKE decreases to a prescribed small value (0.2 m 2 s -2 ). The prognostic equation for TKE is solved by using diagnostic estimation of potential temperature, water vapor variance, and covariances [28].

Relation of ceilometer PBL with meteorological variables
The influence of selected meteorological conditions on the derived ML profiles from the ceilometer measurements was studied using the data for the period starting 12 February 2021 to 19 February 2021.
The selected meteorological variables from the DPIE Lidcombe monitoring station shows that temperature and solar radiation are strongly correlated with MLH as shown in Figure 3 and Figure 4. This is expected as during daytime, the growth of the ML is mainly driven by thermal convection.
In comparison with relative humidity (RH), the PBL measurements from ceilometer is negatively correlated with RH ( Figure 4c). The result is consistent with a study by Dang et al. 2016 [29] at Lanzhou, a semi-arid area in northwest China where they found the MLH as measured by a Micro-Pulse Lidar instrument (MPL-4, Sigma Space) was negatively correlated with relative humidity. Allabakash and Sanghun 2020 [30], in their study of climatology of PBLH-controlling meteorological parameters, have also shown that RH is negatively correlated to PBLH over Korea. Australia is a dry continent; high temperature is usually associated with low RH. When temperature is low and the RH is high, the sensible heat flux will be reduced and hence the ML growth is slower. The inverse relation between temperature with relative humidity explains the MLH is positively correlated with temperature but negatively correlated with RH.
For Merriwa, the results are similar to those observed at Lidcombe in terms of correlation of MLH with temperature, solar radiation and relative humidity. The derived ML profiles from CL51 lidar were also used to examine the relationship with the ozone profiles. This comparison illustrated in Figure 5 which shows a strong correlation between ozone and ML height. The increase in the MLH is associated with the increase in temperature. It is known that during daytime the temperature increases will result in higher photochemical reaction rate and higher ozone production, even though the MLH increase also reduces the ground ozone level due to vertical mixing. For this analysis period from 12 February 2021 to 19 February 2021, the photochemical rate production is higher than the dilution factor from vertical mixing process. For other pollutants, such as PM10, PM2.5, NO2, there are no relation between these pollutants at ground level and MLH for the austral summer period 12 to 19 February 2021 with correlation R 2 of 0.192, 0.0068 and 0.0048 respectively. In a study of impact of mixing layer height on air quality in winter in Delhi, India, Murphy et al. 2020 [31] analysed CL51 ceilometer MLH data and ground monitoring data. They found ozone is strongly correlated with MLH but PM2.5 has less correlation with MLH. Our results however are more consistent with [32] in their study on relation between MLH from CL51 and ambient pollutant concentration measured on ground station in Berlin. Their results showed the correlations between MLH and concentrations of pollutants (PM10, O3 and NOx) for different locations in Berlin were varied. There was no clear pattern in correlations with PM10 which were quite different for different sites whereas the correlation with NOx seems to depend on the vicinity of emission sources in main roads [32]. Only in the case of ozone, a clear correlation was found. This is not surprising as ozone is a secondary pollutant formed from photochemical process driven by sunlight and occurred over regional scale and hence has longer correlation distance than other pollutants [33].

Ceilometer MLH and model forecast PBL comparisons
Comparison of the MLH as measured at Merriwa and the predicted PBLHs from CCAM-CTM and WRF-CMAQ for the period 12 to 19 February 2021 and 20 April 2021 to 3 May 2021 based on performance metrics of mean bias (MB), normalised mean bias (NMB), root mean square error (RMSE), correlation coefficient (r) and index of agreement (IOA) is summarised in Table 1. Overall, both models underpredicted the PBLH (MB is negative) and the PBLH as predicted by WRF performed slightly better than CCAM at both sites. And model prediction at Merriwa is better that that at Lidcombe. Another comparison analysis was also performed for an autumn period from 20 April 2021 to 3 May 2021. The time series of the predicted PBL profiles by the selected numerical models and the derived MLH from the ceilometers during the period from 20 April 2021 to 2 May 2021, are shown in Figures 5 and 6. There is a reasonable agreement between the models and the observation. At Lidcombe the peak of the MLH occurred at around mid-day from 12:00 to 15:00 hours as expected, usually at the height from 1500 to 2000 metres. Both CCAM-CTM and WRF-CMAQ have slightly underpredicted the MLH peaks. The model predictions of the nocturnal PBL heights were lower than observations highlighting the complex structure of this layer during this period. Shi et al. 2020 [4] also reported the nocturnal boundary layer heights are seriously underestimated by WRF model in a study in Beijing from 26 to 31 Dec 2017. At Merriwa when compared to ceilometer observations for the period between 20 April 2021 to 2 May 2021, both the air quality models CCAM-CTM and WRF-CMAQ predicted reasonably well the PBLH. Similar to Lidcombe site, the nocturnal PBLH predictions were underpredicted. In particular, on the 23 and 24 April 2021, the estimated nocturnal MLH was high at mostly above 500m indicating warmer air above ground. But the prediction of CCAM-CTM PBLH is less than 100m while WRF-CMAQ nocturnal PBLHs varying between 50m to 400m on 23 April 2021 and are reasonably good as compared to observation. Similarly, WRF-CMAQ nocturnal PBLHs on 24 April 2021 varied between 50 to 600m. The underprediction of nocturnal PBLH in the two air quality models can also be due to MLH being assigned erroneously to the residual layer height rather than the minimum height in the stable layer below it in the minimum gradient detection algorithm used in CL51 ceilometer [32].
Compared to Lidcombe, the CL51 MLH measurements for this summer period showed the maximum MLH at about 1500m compared to 1500-2000 m at Lidcombe during daytime convection condition.
Overall, the performance of the WRF-CMAQ and CCAM-CTM, as currently configured and used for forecasting, in PBLH prediction when compared to CL51 MLH observation is reasonably good based on the performance indicators presented above. The results for CCAM are consistent with CCAM vertical wind and temperature prediction in the previous study [34]. The difference in PBLH prediction between WR-CMAQ and CCAM-CTM can also be due to the different meteorological host data (NCEP FNL and ACCESS-R) used and in the boundary height calculation schemes, Mellor-Yamada-Janjic Scheme (MYJ) in WRF and bulk Richardson number scheme in CCAM.
The bulk Richardson number scheme is also implemented in YSU (Yonsei University) non-local scheme in WRF. And the PBLH calculation in the MYJ scheme in WRF is based on TKE (Turbulent Kinetic Energy) parametrisation while the ceilometer MLH is derived from the aerosol backscatter profile minimum gradient.  Table 1. The chemistry option in their WRF-Chem was based on Regional Acid Deposition Model, version 2 (RADM2) and Modal Aerosol Dynamics Model for Europe/Secondary Organic Aerosol Model (MADE/SORGAM) as the aerosol module.
They also compared the PBLH, as derived from vertical aerosol profile prediction from WRF-Chem using the same wavelet transform on profile data as used in HSRL back scatter data, with HSRL ML height data. Reasonable agreement was also obtained. For example, with Sacramento campaign dataset, the correlation between WRF-Chem and HSRL ML heights across all flights was 0.6 (r value), mean bias difference of 194m and RMS difference (or error) of 586m.
The above comparison results using two different methods of calculating the MLH and PBLH (one using potential temperature and one using aerosol backscatter) is very similar. This result from their study has important implication. It means that PBLHs determined using temperature profile and those using aerosol profile are equivalent. The PBLH prediction in WRF-Chem using MYJ scheme based on turbulent kinetic energy parametrisation gives similar result to that based on aerosol profile peak gradient method, such as the Haar wavelet transform.
Milovac et al. 2015 [36], in their study of comparing different boundary layer schemes and land surface model (LSM) used in WRF in Germany, found that the land surface processes can have an impact not only to the lower CBL but also extend up to the interfacial layer and the lower troposphere. Using 6 different PBL and LSM schemes (including MSY PBL scheme and Noah land surface model as used in our study), the impact of diurnal change in humidity profiles is more significant at the interfacial layer than close to the land surface. They concluded that the representation of land surface processes has a significant impact on the simulation of mixing properties within the CBL.
Other uncertainties include the detection algorithm and MLH gradient method used in processing the lidar echo signals from ceilometers. For example, Bedoya-Velásquez et al. 2021 [37] recently have improved the range corrected signal (RCS) affected by water vapor absorption as it was found that raw ceilometer signal overestimates the water vapor corrected one, mainly below 1 km AGL. Vertical water vapor data if not available, then relative humidity data profiles can be obtained from GDAS (Global Data Assimilation System) database to correct the RCS. And as Geiß et al. 2017 [32] have pointed out, the MLH as determined from layer aerosol profile can give error as the residual layer is misinterpreted as mixing layer during nighttime.
Ceilometer data signal processing and PBL analysis can be enhanced using layer information from MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, version 2) such as vertical aerosol concentration from which lidar optical parameters, such as aerosol extinction and backscattering coefficients, can be calculated.

Case studies at Merriwa
The Merriwa ceilometer has been installed and operational since 16 September 2020. Since then the measurement and its performance in determining the ML height and cloud base height at various hours of the day provided us with an understanding of the evolution of the ML above ground. The following case study analysis of some PM10 event days at Merriwa illustrates the use of ceilometer data to understand and interpretate the source and cause of these high PM10 concentration as measured at Merriwa monitoring station.
The concentrations of PM10 measured at Merriwa monitoring site from 1 July 2020 to 22 September 2021 are shown in Figure 6. The time series of concentrations show no high PM10 events after 12 March 2021 (the date when the Merriwa ceilometer was put in operation). However, there was occasional spikes for a few hours with lower levels. Before March 2021, there were a few particle events with a major dust storm in 20 August 2020 (large spike above 250 μg/m3 in Figure 6). Other smaller particle events happened on 16 October 2020, 10 December 2020 and 19 January 2021.
It is anticipated that dust events could happen in the coming months during spring and summer 2021-2022 and the ceilometer will provide informative measurements to detect the dust plumes. Similar to the study [10] of dust detection, using lidars and ceilometers from volcanic ashes, with case studies of event days where high ground-level PM10 concentrations occurred, in the present analysis, we also analyse the 3 particle events on 16 October 2020, 10 December 2020 and 19 January 2021 at Merriwa to determine where the sources are and the causes of these events at the receptor. These 3 case studies show the influence of meteorology on the transport of emitted aerosols from various sources to the receptor site and the complexity of meteorological processes in understanding the pattern of aerosol profile at the receptor site. The CL51 ceilometer at Merriwa has the profile measurements for those event days which are analysed in detail as described below with corroboration from satellite data and back trajectory analysis.

Event day 16 October 2020
The 16 October 2020 ceilometer measurement is shown in Figure 7. The peak PM10 concentration was detected at 18:00 AEST. This day was cloudy as detected by satellite MODIS Terra/Aqua above Merriwa. The thick foggy water droplet from 5am to 9am was close to the ground on the morning of 17 October 2021 as seen in Figure 7. On 16 October 2020, the MLH as measured by the ceilometer dropped from ~2000m to ~500m at 18:00 AEST. This corresponds to the high PM10 concentration measured at Merriwa monitoring station. Back trajectory analysis for 48 hours was performed to determine the source of this high PM10 concentration detected at Merriwa. Figure 8 shows the source of this high PM10 concentration as from the north west and north of Merriwa where satellite data showing high AOD in this area. The HYSPLIT back trajectory analysis was run at 10m, 20m and 50 m above ground level (AGL) with GDAS (Global Data Assimilation System) 1 o resolution data from NOAA. For the 10 December 2020, high PM10 occurred at 17:00 AEST. Similar to the above analysis, back trajectory analysis was performed and ceilometer data on this day is shown in Figure 9. The MLH increased after 17:00 AEST and only decreased about 3 hours later. This shows that the vertical profile of atmospheric structure had no influence on the high peak ground level concentration at Merriwa. The advection transport of aerosols was the main cause. This is confirmed by HYSPLIT back trajectory analysis. The source of air parcels at 1500m above ground level at Merriwa came from ground level (10m) west of Newcastle 48 hours before as shown by the 'red' trajectory in Figure  10. While air parcels at 10m and 20m AGL at Merriwa came from ground level air parcels north west of Sydney as shown by the 'green' and 'blue' trajectories. There was no change in level heights of these trajectories during the 48 hours transport of aerosols from north west of Sydney passing through the Upper Hunter before reaching Merriwa receptor. In other words, the horizontal advection transport of aerosols, which was emitted along the trajectory path and carried by the wind to Merriwa receptor, caused the high PM10 concentration at 17:00 AEST on 10 December 2020.
(a) (b) Figure 10. (a) HYSPLIT back trajectory analysis for 48 hours at Merriwa on 10 December 2020 at 17:00 AEST showing the sources of air parcels at 10m (green), 20m (blue) and 1500m (red) came from the south east of Merriwa where high columnar AOD (yellow) was detected from MODIS Terra/Aqua satellites sensors (b) Details path trajectories and height AGL for 48 hours.

Event day 19 January 2021
For the high PM10 day of 19 January 2021 at 15:00 AEST, HYSPLIT back trajectory analysis at 10m, 20m and 50m above ground level at Merriwa indicated that the source of the aerosols came from the south east and south. There was no cloud above Merriwa and the rest of NSW on this day, except in the Greater Sydney Metropolitan near the coast.
The ceilometer data at Merriwa on this day showed MLH was reduced from about 1500m to less than 1000m at about 16:00 AEST ( Figure 11). The paths of the back-trajectory analysis at 10, 20 and 50m AGL are at similar levels when passing over the Upper Hunter before reaching Merriwa. While at 1500m AGL at Merriwa, the source of air parcels was from the south west and travelled at height above 1500m as shown in Figure 12. From the 3 case studies above, it has been shown that transport of aerosols to Merriwa receptor can be from western or northern NSW (16 October 2020) but can also be from south east in the Upper Hunter or southern NSW (10 December 2020 and 19 January 2021). Dust is most likely from western and northern NSW while anthropogenic sources from south east in the Hunter region and from southern NSW.
For future study, space-born lidar data from CALIPSO satellite can be used to supplement information on aerosol vertical structure above Merriwa and Lidcombe. Even if there was missing ceilometer data, such as on 20 August 2020, but CALIOP lidar sensor onboard CALIPSO satellite can provide some information on vertical structure above and near the Merriwa region on day when the satellite passed over the region. For the high PM10 events that occurred, ceilometer data and CALIPSO data can be used together to corroborate the measurement results. Li et al. 2021 [38], in their study on using the Lufft CHM15K ceilometer for retrieval aerosol characteristics including AOD, have used CALIPSO level 1 attenuated backscatter coefficient profile to compare and validate with the calibrated ceilometer backscatter coefficient profile Recent usage of ground-based lidar systems such as MPLs and ceilometers attracts more attention on the application of such data to climate model validation such as in CMIP6 (Coupled Model Intercomparison Project Phase 6) for cloud feedbacks in Cloud Feedback Model Intercomparison Project (CFMIP). Kuma et al. 2021 [39] have developed a tool, called the Automatic Lidar and Ceilometer Framework (ALCF), which allows the value of large ALC (Automatic Lidar and Ceilometers) networks worldwide such as ICENET, E-PROFILE and MPLNET to be used for model validation and other scientific studies. This ALCF tool will allow the analysis of lidar data with various methods to extract information on cloud and aerosol structure at many sites.
Besides numerical method based on air quality model to predict PBLH, data-driven methods can be used in future study to predict the PBL based on ceilometer observation at Merriwa and Lidcombe. The methods include different regression techniques (multiple regression, random forest, Ridge regression, decision tree learning, ...) with predictors such as temperature, humidity, time series analysis (ARIMA, SARIMA), ensemble methods (LightGBM, AdaBoost) and Artificial Intelligence (AI) methods (LSTM, Vector output model, Encoder-Decoder method). Recently, AI methods, such as K-means unsupervised algorithm and the AdaBoost supervised algorithm, have been used to determine the MLH based on lidar backscatter profiles [40].

Conclusions
In this study, we used the attenuated backscatter profiles data from the two ceilometers located at Merriwa and Lidcombe to characterise the vertical and horizontal mixing aerosols in the boundary layer. The ceilometers profiles were used to derive the structure of the boundary layer which includes the mixing layer, the nocturnal residual layer and the elevated aerosol layer at these locations. These ceilometers have successfully detected the aerosol dispersion and mixing height profiles above the sites. We have compared these results with the simulated results carried out for the same period, using WRF-CMAQ and CCAM-CTM models. The results show that the mixing height derived boundary layer detection algorithm from the two ceilometers compared satisfactory to those predicted by the numerical air quality models, CCAM-CTM and WRF-CMAQ, used by DPIE for forecasting daily air pollution and development of policy scenarios.
As well as used to verify the PBLH as predicted by air quality models, the ceilometer data also are useful in explaining the nature of aerosol layers change in the vertical atmospheric structure. This dynamic change allows one to explain the transport of aerosol in particular and other pollutants in general above one location as was shown in the trajectory analysis of aerosols at Merriwa in this study.