Near-real-time processing of a ceilometer network assisted with sun-photometer data : monitoring a dust outbreak over the Iberian Peninsula

The interest in the use of ceilometers for optical aerosol characterization has increased in the last few years. They operate continuously almost unattended and are also much less expensive than lidars; hence, they can be distributed in dense networks over large areas. However, due to the low signal-to-noise ratio it is not always possible to obtain particle backscatter coefficient profiles, and the vast number of data generated require an automated and unsupervised method that ensures the quality of the profiles inversions. In this work we describe a method that uses aerosol optical depth (AOD) measurements from the AERONET network that it is applied for the calibration and automated quality assurance of inversion of ceilometer profiles. The method is compared with independent inversions obtained by colocated multiwavelength lidar measurements. A difference smaller than 15 % in backscatter is found between both instruments. This method is continuously and automatically applied to the Iberian Ceilometer Network (ICENET) and a case example during an unusually intense dust outbreak affecting the Iberian Peninsula between 20 and 24 February 2016 is shown. Results reveal that it is possible to obtain quantitative optical aerosol properties (particle backscatter coefficient) and discriminate the quality of these retrievals with ceilometers over large areas. This information has a great potential for alert systems and model assimilation and evaluation.


Introduction
Atmospheric aerosol is one of the main responsible factors of climate radiative forcing through multiple processes including aerosol-radiation and aerosol-cloud interactions (IPCC, 2014).The aerosol direct effects depend on the optical properties and spatial and vertical distribution of the aerosol in the atmosphere.In spite of the recent advances on instrumentation that has improved the ability of characterizing key aerosol properties and increase the spatial resolution, the associated uncertainties are still considered to be one of the majors in climate forcing (Boucher et al., 2013).
In this sense, the implementation of observational networks is crucial for spatial characterization of aerosol properties.Ground-level aerosol measurement networks represent key tools in the study of aerosol radiative forcing.These observational networks provide surface measurements Published by Copernicus Publications on behalf of the European Geosciences Union.distributed over large areas, e.g., the Global Atmospheric Watch, GAW (GAW, 2011), and ACTRIS (www.actris.eu)for Europe.In addition, one of the recognized instruments for the retrieval of column-integrated aerosol properties is the robotic sun and sky photometer that is used in the global Aerosol Robotic NETwork (AERONET; Holben et al., 1998;Dubovik et al., 2006).Lidar systems are well-known active remote sensing instruments for the vertically resolved characterization of aerosol optical and microphysical properties (Winker et al., 2003).GAW Atmospheric Lidar Observation Network (GALION) has emerged as an initiative of the GAW aerosol program (GAW, 2008).Its main objective is to provide the vertical component of the aerosol distribution through advanced laser remote sensing in a network of ground-based stations.Among other networks, GALION includes the European Aerosol Research Lidar Network (EARLINET) that provides vertical aerosol profile observations over Europe based on 27 instruments in 16 countries (Pappalardo et al., 2014), the Micro-Pulse Lidar Network, MPLNET (Welton et al., 2001), and the Latin American Lidar Network, LALINET (Guerrero-Rascado et al., 2016).
In order to obtain a larger spatial coverage than groundbased networks, in the last few years some space missions have been promoted focusing on aerosol measurements from satellites, e.g the Lidar in Space Technology Experiment, LITE (McCormick, 1997), and the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation, CALIPSO (Winker et al., 2003).The main disadvantage of measurements from spaceborne platforms is the low temporal resolution, since the measurements are limited to the satellite passes over a region.
The complexity of lidar systems requires staff to be trained in their operation, and the analysis procedures are not fully automated in many stations.In this sense, continuous operation of lidar systems is not feasible for most stations.In addition, economic and operational costs hinder the implementation of dense lidar networks.On the other hand, ceilometers are one-wavelength (near infrared) lidars with simple technical specifications (eye-safe low pulse energy and high pulse repetition frequencies) allowing for unattended and continuous operation.Originally designed for cloud base determination, their performance has been improved in the last few years.Their capabilities have been shown for determining planetary boundary layer (e.g., Wiegner et al., 2006;Münkel et al., 2007;Haeffelin et al., 2012;Pandolfi et al., 2013), detection and forecast of fog (Haeffelin et al., 2016), and recent efforts have been conducted to quantify the aerosol optical information that can be derived from ceilometers (Frey et al., 2010;Heese et al., 2010, Wiegner et al., 2014).
The main advantage of the use of ceilometers for aerosol characterization is, on the one hand, the automatic and much simpler operation compared to lidars and, on the other hand, the possibility of installing them distributed over large areas.Meteorological services such as those in Germany, France, the Netherlands or the United Kingdom are deploying ceilometers networks to cover their national territories with the objective of reaching a spatial density of nearly one device every 100 km (e.g., de Haij and Klein-Baltink, 2007;Flentje et al., 2010).Due to a dense number of instruments and continuous measurements, operative networks need an automated processing and a protocol that ensures the quality of the data.
In this sense, two programs in Europe are dealing with the use of automated lidars and ceilometers for aerosol and cloud properties characterization.The COST Action ES1303 TOPROF (TOwards operational ground based PROFiling with ceilometers, doppler lidars and microwave radiometers for improving weather forecasts) aims in one of its working groups at better characterizing the parameters that can be derived from ceilometer measurements and related uncertainties.At the same time, E-PROFILE, a program of EU-METNET (EUropean METeorological services NETwork), focuses on the harmonization of ceilometer measurements and data provision across Europe.
In this study we present the implementation of procedures to manage a regional ceilometer network for aerosol characterization over the Iberian Peninsula, the Iberian Ceilometer Network (ICENET).An automatic calibration procedure is applied to the ceilometers and this calibration is used to validate the elastic inversion automatically applied to the profiles.This method uses additional aerosol optical depth (AOD) information during the calibration for the quality assurance of the data.
All processes can be performed unattended and in nearreal time with the objective of obtaining reliable vertically resolved aerosol optical properties.This information is especially useful for strong events, such as mineral dust outbreaks, volcanic plumes, severe biomass burning episodes or contamination episodes.Thus, the aerosol information obtained can be potentially used as an alert system for aviation or weather services or to feed models for assimilation and validation in near-real time.
The capabilities of this distributed network are explored by characterizing an unusually intense dust outbreak affecting the Iberian Peninsula on 20-24 February 2016 and a multi-wavelength (MW) Raman lidar is used to validate the retrievals from ceilometers.
The next section describes the Iberian Ceilometer Network and related instrumentation used in this study.Section 3 presents the methodology, including the calibration of ceilometers (Sect.3.1) and the use of the calibration for inversion validation, as well as a validation with an independent lidar system (Sect.3.2).The results are presented in Sect.4, with a description of the dust event (Sect.4.1 and 4.2) and retrievals from ceilometers (Sect.4.3).Finally, conclusions are presented in Sect. 5.

Instrumentation: the Iberian Ceilometer Network
An initiative of the Atmospheric Physics Group of the University of Granada has been the coordination of a network of ceilometers (ICENET) combined with sun photometers for the characterization of atmospheric aerosol with the objective of obtaining reliable vertically resolved aerosol optical properties in near-real time.The first goal is obtaining the total attenuated backscatter for all ceilometers in the network, i.e., to obtain calibrated output from ceilometers, and the second one is applying an inversion algorithm to the ceilometer profiles in order to obtain the particle backscatter coefficient.All sites of this new network have a co-located AERONET CE318 sun-sky photometer (Cimel Electronique) that is used to constrain the ceilometers calibration and inversion retrievals.In addition, the high-performance lidar system MULHACEN, located at the EARLINET Granada station, is used as an independent validation of the inversions.This nested approach combining high-performance systems like those operated in EARLINET and the distributed ceilometer plus sun photometer is an example of synergy among active and passive remote sensing observations in the ACTRIS research infrastructure (www.actris.eu).
Figure 1 shows a map of the ceilometer distribution over the Iberian Peninsula, and Table 1 presents the characteristics of each site.
All sites operate a Jenoptik (now Lufft) CHM15k-Nimbus ceilometer and have a co-located AERONET sun photometer, except Montsec station (MSA), which has the photometer 770 m above the ceilometer and at a horizontal distance of 2 km approximately (Titos et al., 2017).The ceilometer at Murcia (UMH) was not operative during the outbreak studied in this work.
The CHM15k is a ceilometer that operates with a pulsed Nd : YAG laser emitting at 1064 nm.The energy per pulse is 8.4 µJ with a repetition frequency in the range of 5-7 kHz.The laser beam divergence is less than 0.3 mrad and the laser backscattered signal is collected on a telescope with a field of view of 0.45 mrad.The signal is detected by an avalanche photodiode in photon-counting mode.Complete overlap of the telescope and the laser beam is found about 1500 m above the instrument (Heese et al., 2010).According to the overlap function provided by the manufacturer, the overlap is 90 % complete between 555 and 885 m a.g.l.The vertical resolution used is 15 m and the maximum height recorded is 15360 m a.g.l.Ceilometers at Granada (UGR), Tabernas (PSA) and Valladolid (UVA) operate at a temporal resolution of 15 s, while ceilometers at Montsec (MSA) and Badajoz (UEX) operate at a temporal resolution of 1 min.
The process of calibration for ceilometers described on the next section is assisted with AOD data from co-located AERONET stations.All sun photometers near the ceilometers belong to the Iberian network for aerosol measurements (RIMA), a regional network associated with AERONET.This means that all instruments are routinely calibrated following the same protocol and the data are quality-controlled.The sun photometer provides solar extinction measurements at 340, 380, 440, 675, 870, 936 and 1020 nm, allowing for computing the AOD at these wavelengths (except 936 nm).The AOD uncertainty ranges from ±0.01 in the infraredvisible to ±0.02 in the ultraviolet channels (Holben et al., 1998).For comparison with the ceilometers the AOD is extrapolated to 1064 nm by the Ångström law (Ångström, 1964) using the AOD measurements at 870 and 1020 nm.Level 1.5 AERONET data, which are automatically cloudscreened and delivered in near-real time, are used in this analysis.
At UGR station a multi-wavelength Raman lidar system (MULHACEN) is used for validation of the ceilometer inversions.The upgraded LR331-D400 (Raymetrics Inc.) operated at IISTA-CEAMA (Andalusian Institute for Earth System Research) has been part of EARLINET since April 2005.This lidar system is a ground-based, six-wavelength system with a pulsed Nd : YAG laser.The emitted wavelengths are 355, 532 and 1064 nm with output energies per pulse of 60, 65 and 110 mJ, respectively.It has elastic backscatter chan- ).Full overlap is reached around 1220 m a.g.l., although the overlap is complete at 90 % between 520 and 820 m a.g.l.(Navas-Guzmán et al., 2011).Appropriate overlap corrections are derived following the procedure of Wandinger et al. (2002).

Methodology
The principle of measurement for elastic lidars and ceilometers is the same, and retrieval of optical properties in both systems follows the lidar equation (the dependency with the wavelength has been omitted for simplicity since it is always the same in ceilometers): where P (z) is the backscattered power received in the telescope from a distance z, C L is a parameter that depends on the geometry and characteristics of the instrument and universal constants, and the term z 2 accounts for the acceptance solid angle of the receiver optics with the distance to the laser.The backscattered signal collected by the telescope depends on the overlap between the laser beam and the telescope field of view, and the degree of overlap is quantified by O (z), ranging from 0, if there is no overlap, to 1, if overlap is complete.β (z) is the atmospheric backscatter coefficient and T (z) estimates the atmospheric transmittance of the laser signal (squared due to travel back and forth).Also, both properties can be split into contributions of particles and molecules (β (Fernald, 1984).
In Eq. ( 1) the only properties depending on the medium are β(z) and T (z).Thus, the atmospheric attenuated backscatter is defined as (2)

Ceilometer calibration
The ceilometers used in this study provide the rangecorrected signal (RCS(z) = P (z) • z 2 ) as output, using an overlap function determined by the manufacturer and corrected for the number of laser shots.Therefore, the only parameter that needs to be addressed is C L .Wiegner et al. (2014) describe a method to find the C L parameter in ceilometers, commonly referred as ceilometer calibration.This method compares the RCS from the ceilometer in a particle-free region with the molecular attenuated backscatter that can be calculated using Rayleigh theory.The Rayleigh fit compares the gradient with altitude (the slope) of both profiles and looks for a region in the ceilometer profile that has the same trend as the expected molecular profile.In this study, we select regions of 990 m with a difference in gradients below 1 %.Thus, in that region or reference height (z ref ), C L can be calculated: At this reference height, the backscattering is only due to molecules.The transmittance due to molecules (T m ) can be easily determined from Rayleigh theory but the transmittance due to particles (T p ) is unknown.However, if a co-located sun photometer is available, T p can be calculated, using the AOD at 1064 nm: When trying to automate the calibration process, the main problem is that z ref must be a particle-free region and, due to the low signal-to-noise ratio, finding z ref is not always possible.In some cases, the region might be a non-particle-free region that, on average, follows the molecular trend (they have a similar gradient).Also, we might find several regions that meet the criteria, but it is complicated to discriminate automatically which one is the most appropriate.A way to ensure that the z ref selected is a molecular region is by applying the Klett-Fernald inversion algorithm (Klett, 1981(Klett, , 1985;;Fernald et al., 1972;Fernald, 1984) as follows.
First, we need to determine the z ref .Thus, after finding z ref , the Klett-Fernald inversion is applied.Heese et al. (2010) and Wiegner et al. (2012Wiegner et al. ( , 2014) ) showed the capabilities of ceilometers by applying this inversion algorithm to study a few cases.Using the AOD measurements, the lidar ratio (Lr) of the inversion can be adjusted.This can be done by matching the integral of the particle extinction coefficient profile (i.e., particle backscatter coefficient profile multiplied by the Lr) with the AOD.Wiegner et al. (2012) applied this procedure to a Jenoptik CHM15kx ceilometer, obtaining reasonable values for the Lr.
In summary, the calibration process is carried out as follows: 1. First, temporal averaging of the profiles is performed (hourly averages are used for the calibration).The first 300 m of the profile are assigned to the value at 300 m to avoid large overlap correction.
2. Second, for each profile a set of potential z ref is obtained by comparing the profiles of the RCS and β m , which is obtained from a standard atmosphere profile scaled to ground temperature and pressure.The slopes are calculated over a 990 m window.All regions with slope differences below 1 % are selected.
3. For each z ref , and Lr from 20 to 80 sr, Klett-Fernald inversion is applied and the resulting profile for the backscattering coefficient is integrated and multiplied by the Lr and compared to the AOD.The pair z ref and Lr that minimizes the difference between the integral of the particle extinction coefficient profile (i.e., particle backscatter coefficient profile multiplied by the Lr) with the AOD is selected.
4. Finally, C L is calculated using Eq.(3) if the minimum difference calculated in step 3 is below 10 %.
C L calculated with this method uses Eq. ( 4) to calculate T 2 p .In the case of MSA the sun photometer is 770 m above the ceilometer and the T 2 p calculated would not be representative of the entire column.However, MSA is a remote mountain site and the effect of accounting for T 2 p values obtained using the AOD measured from 770 m above the ceilometer can be considered negligible.Thus, the total attenuated backscatter can be calculated by applying the following equation: Calibration values can be used individually or averaged over a period of time.Table 2 shows a mean calibration factor (± standard deviation) calculated using this method for all sites in ICENET for the period 1 May 2014 to 1 May 2016.

Ceilometer inversion
Total attenuated backscatter obtained by applying the calibration factor allows the comparison between ceilometers since the signal is corrected for instrument characteristics.Also, a long time series of the calibration allows determining possible problems or degradation of the systems.However, the total attenuated backscatter is influenced by transmission, so in order to be able to monitor and compare singular events at multiple sites, the backscattering coefficient is more appropriate.Section 3.1 showed that it is possible to apply the Klett-Fernald inversion to ceilometer data, but the challenge is to determine automatically, without human supervision, whether the inversion is correct or not.A common step between the calibration proposed in Sect.3.1 and Klett-Fernald inversion is finding a reference height (Z ref ).At the Z ref selected for the inversion with Rayleigh fit, applying Eq. (3), we can obtain a value that has to be close to the C L calculated for the instrument on a longer period of time.If a simultaneous AOD measurement is available, the calibration process itself provides the inverted backscattering coefficient profile (steps 3 and 4 of the calibration process) and the inversion can be marked as valid or invalid based on the value of the C L compared to a long-term C L value.If no simultaneous AOD measurement is available (e.g., during nighttime or partially cloudy skies), an approximation of the AOD needs to be used in order to apply Eq. ( 4).In this case, an interpolated value or an averaged value can be used.
The next section quantifies the differences between these backscattering coefficient inversions and the inversions calculated independently with a multi-wavelength Raman lidar.

Lidar -ceilometer comparison
During the dust outbreak affecting the Iberian Peninsula between 20 and 24 February 2016, the multi-wavelength lidar operated on 22 February between 07:30 and 14:00 UTC and on 23 February between 08:00 and 13:30 UTC.Elastic inversions using the Klett-Fernald method were applied to 30 min average profiles at 1064 nm using a fixed lidar ratio of 50 sr.Thus, a total of 24 particle backscatter coefficient profiles were obtained.Coherence of the inversion at 1064 nm was checked against the Klett-Fernald and Raman methods at 355 and 532 nm.The resolution of the multi-wavelength lidar (7.5 m) has been downscaled to 15 m for the comparison with the ceilometer.
The ceilometer elastic inversion, using the Klett-Fernald method, was also applied to 30 min average profiles for the same period; a total of 15 profiles were successfully inverted (a reference height was found automatically).The calibration factor at the reference height was calculated using the average AOD for the entire dust event.If negative C L values are discarded, a total of 11 profiles are comparable with lidar inversions.
Each one of the derived C L values at the reference height selected for the inversion is compared with long-term C L calculated for the Granada station ceilometer over a long period of time (see Table 2), classifying the situation according to statistical parameters measuring the agreement between the ceilometer and lidar retrievals.The normalized mean bias (NMB) in particle backscatter of the ceilometer and lidar profile is calculated following Eq.( 6).The center of mass of the profiles is calculated with Eq. ( 7), as is the relative difference between ceilometer and lidar center of mass.Finally, the coefficient of correlation (R) of the profiles is determined.
In Eq. ( 6), β ceil and β lidar are the mean particle backscatter coefficient from ceilometer and lidar, respectively, for the entire retrieved profile, and β in Eq. ( 7) may refer to ceilometer or lidar particle backscatter coefficient depending on the case.
Figure 2 shows for the 11 comparable profiles, the retrieved calibration factors at the reference height on the ceilometer profile versus the NMB (panel a), center of mass relative differences (panel b) and R (panel c).It is evident that ceilometer profiles with a calibration factor closer to the mean calibration factor have inversions closer to the lidar inversions.Figure 2 also shows that, for the statistics NMB and R, the difference between the calibration factor and the mean calibration factor is related, and the farther the profile calibration factor is from the mean value, the worse the mentioned statistics.Thus, it seems feasible to determine the quality of the profiles by selecting an appropriate threshold for this difference.Considering a maximum discrepancy between the particular calibration factor and the long-term calibration factor equal to 1 standard deviation of the mean value of the calibration factor (dotted lines in Fig. 2), we obtain four profiles that have a NMB smaller than 15 %; the center of mass of the profiles is practically the same, with a relative difference smaller than 2 %, and finally R of the profiles is above 0.92.
A sequence of ceilometer and lidar particle backscatter profiles from 23 February 2016 is shown in Fig. 3.The first ceilometer profile (marked in blue) has a calibration factor of 2.57 × 10 11 (m 3 sr) and hence is rejected according to the threshold described above.For this case, the NMB of the ceilometer and lidar profiles is −0.31, the center of mass relative difference is −0.06 and the R is 0.84.The other four ceilometer profiles (marked in red) have calibration factors within the standard deviation of the mean calibration factor.The profiles on 23 February at 09:00 and 09:30 UTC correspond with a decoupled dust layer.Those profiles have a NMB of −0.08 and 0.1, respectively, the center of mass relative difference is −0.01 and −0.03, respectively, and R is 0.95 and 0.97, respectively.The profiles on 23 February at 12:00 and 12:30 UTC show that the previous dust layer is mixed with the boundary layer.In these cases, profiles have a NMB of 0.14 and −0.12, respectively, the center of mass relative difference is 0.006 and −0.01, respectively, and R is 0.99 and 0.93, respectively.

Results
The capabilities of the ceilometer network for aerosol optical properties characterization and the near-real-time processing have been tested with the analysis of the African dust outbreak that affected the Iberian Peninsula on 20 February 2016 and persisted until 24 February 2016.Sorribas et al. (2017) studied the same event and compared it with meteorological parameters, aerosol properties and ozone from historical data sets on a site in southern Spain.They concluded that the event was exceptional because of its unusual intensity, its impact on surface measurements and the month of occurrence.In addition, Titos et al. (2017) also analyzed this event using 250 air quality monitoring stations over Spain to investigate the impact and temporal evolution of the event on surface PM 10 levels.They also investigated aerosol optical properties, including attenuated backscatter from ceilometer during the event at Montsec station (one of the station included in ICENET).They concluded that the impact on surface PM 10 was exceptional and highlighted the complexity of the event.Next section provides a detailed description of the dust mobilization and arrival of the plumes to the Iberian Peninsula.

Description of the dust episode
The evolution of the dust outbreak is illustrated in Fig. 4, where a sequence of the false-color RGB dust images from MSG-SEVIRI is shown.This product makes use of three thermal satellite channels to contrast the brightness temperature signal between surface, cloud and dust (Lensky and Rosenfeld, 2008), in a color scheme in which dust appears in magenta.The presence of dust plumes over the High Plateau located between the Sahara and Tell Atlas in Algeria at 12:00 UTC on 20 February 2016 is shown in Fig. 4a.Dust migrated northwestward and passed over the Alboran Sea from the Algerian-Moroccan border around 14:00 UTC, reaching the southern Iberian Peninsula at 18:00 UTC (Fig. 4b) and continued moving northwestward (Fig. 4c).A second dust plume migrated northwards on 21 February 2016 at 16:00 UTC (Fig. 4d).
SYNOP meteorological observations and aerodrome routine (METAR) and special (SPECI) reports in northern Africa recorded a strong reduction in horizontal visibility (MOR, meteorological optical range), down to 2 km, between 07:00 and 08:00 UTC (20 February) at distant loca-tions situated at the edges of the Great Western Erg in Algeria.At the eastern part of the Western Erg and in the Great Eastern Erg, visibility lowered to less than 5 km between 09:00 and 11:30 UTC.In good accordance with the satellite images, at the Saharan Atlas and the High Plateau area, with heights over 1000 m, visibility less than 2 km was recorded at 10:00 UTC at Mecheria, while at the other stations in that area values went down to 2-3 km at 12:00 UTC.High relative humidity and clouds were found in the westernmost sites, also in agreement with the satellite images.It is remarkable that no significant visibility reduction was reported in the north-facing downslope areas of the Tell Atlas and in the Rif Mountains, close or at the coast of Algeria and Morocco.This indicates that dust was uplifted before passing over the northern slope of the Atlas Mountains and the North African coast.Correspondingly, no station in the southern Iberian Peninsula reported a reduction in visibility when the dust plume reached Spain.
The entrance of dust-laden air masses above the ground level in the Iberian Peninsula is confirmed by a backtrajectory analysis performed with HYSPLIT (Hybrid Single-Particle Lagrangian Integrated Trajectory) model (Draxler and Hess, 1998;Stein et al, 2015) using ERA-Interim data of 0.5 • resolution.The trajectory analysis provides an estimate of the range of heights at which the dustladen air masses passed over the study sites.This is illus-  (Pappalardo et al., 2004;Mattis et al., 2016) and are on the order of 10 −9 (m −1 sr −1 ); therefore, they are not shown.trated in Fig. 5 for two stations: in the south where dust reached first the Iberian Peninsula (Granada) and in the westernmost station, on the cyclonic track of the dust (Badajoz).The height vs. latitude plots show that the dust plumes reaching Granada on 20 February, 18:00 UTC and Badajoz on 21 February, 00:00 UTC, arrived at mid-levels in the lower troposphere after being uplifted in the southern slope of the Saharan Atlas from heights between 500 and 1200 m above the ground in that area (black trajectory in Fig. 5).Those trajectories end at around 3250 and 2500 m a.s.l. in Granada and Badajoz, respectively.Figure 5 also shows that the trajectories reaching Granada at about the height of 2250 m a.s.l.passed over Africa before the dust mobilization took place, while the ones below had no African origin.Finally, the trajectory reaching Granada at 4250 m a.s.l.followed the midupper-level tropospheric circulations.In Badajoz, results are analogous: at 1000 m a.s.l. and below, trajectories arrived from the north; between 1250 and 2000 m a.s.l.they arrived from Africa, but the air parcels were located north of the area where dust was observed on the morning of 20 February.The air parcels reaching Badajoz between 2250 and 3250 m a.s.l. were previously located in the area where dust was being observed, while at upper levels trajectories followed the midupper circulation pattern.
In terms of aerosol load, Fig. 6 shows the time series of AOD at 1064 nm (a) and Ångström exponent between 440 and 880 nm (b) for the entire period of the event for all sites in ICENET.We observe that the increase in AOD and decrease in Ångström exponent correspond with the arrival of the plumes at each site.The strongest part of the event occurs on 21-22 February and the second plume is observed clearly for Granada station on 22 February.For dust events following a similar pattern that the one described here during the period 2005-2010, Valenzuela et al. (2012) reported a maximum AOD of 0.98 at Granada, which is significantly lower than the maximum measured at Granada station during this event.

Synoptic scenario and context
Several atmospheric features at mid-upper levels were relevant for this episode as they promoted instability near the surface and induced dust transport in the lower free troposphere.The first feature is the amplification and break-up of a Rossby wave in the eastern Atlantic resulted in a trough that became isolated as a cut-off low over the Atlas Mountains on 19 February.A shallow cyclone then originated leeward of the Atlas Mountains.From early 21 February, the cut-off low displaced off the Moroccan coast and centered southwest of Cape St. Vincent.On 22 February it decayed bringing the Iberian Peninsula under the influence of the Azores and North African subtropical highs, with dominant zonal flows.The second is an upper-level anticyclone over a wide area centered over Niger-Chad, which during the episode intensi-fied and extended northwards to the western Mediterranean.This high pressure influenced circulation at mid-upper levels in combination with the cut-off low.The third is moisture flux at mid-tropospheric levels, which entered from the central Atlantic into the African continent below 20 • N and was transported to northern Africa (at 400-550 hPa according to the radiosoundings in the area) between the upper-level trough and the high-pressure system.The tropical air masses are well visible in the satellite imagery as an elongated cloud band moving north and eastward, and so are the convective clouds formed ahead of the band.The tropical-extratropical interaction between the advected tropical moisture and the upper-level trough located over the Atlas Mountains is linked to convective precipitation in northwestern Africa; see Knippertz (2003) and references therein.Divergence at upper levels (250 hPa) and low-level (850 hPa) convergence are found over the area where the gust front mobilized the dust on 20 February.The interaction with the Ahaggar Mountains in southern Algeria possibly enhanced convection and low-level instability.Convective precipitation was registered at several locations of eastern Spain when the cloud band passed over the area in the second half of the episode.From 22 February onwards the cloud band and local convective situations were gradually displaced to the Mediterranean, as zonal flows began to dominate.
At low levels, the low pressure that formed in the lee of the Atlas Mountains moved to the SW of Cape St. Vincent on 21 February following the upper-level instability.The low was then intensified and influenced northern Africa and most of the Iberian Peninsula.In addition, high pressures over the western Mediterranean were formed when the Rossby wave train progressed to the east and retreated poleward.Then, the North African high, which was previously located over Libya at 850 hPa, extended to Tunisia and Algeria and was gradually intensified in connection with the northward extension of the high pressures at upper levels, which arrived (along with the cloud band on its western flank) at the western Mediterranean Basin.
The advection of dust-laden air masses to the Iberian Peninsula was driven by both the low located to the SW of the Iberian Peninsula and the North African high.The presence of these two synoptic systems corresponds to one of the typical synoptic situations leading to dust transport over the Iberian Peninsula (Rodríguez et al., 2003;Escudero et al., 2005).During the episode, however, two distinct strong plumes were transported from northern Africa to the Iberian Peninsula in consecutive days and showed a different evolution.Dust mobilized by the gust front on 20 February south of the Saharan Atlas and north of the Ahaggar Mountains migrated west and northward to the Iberian Peninsula, as shown in the satellite images, forming a curve-shaped plume over Iberia due to the cyclonic shear imposed by the low.The second strong dust plume was mobilized and transported northwards on 21 February on the western side of the North African high, driven by the intensification of this high-Figure 5. Evolution of the air parcels reaching Granada on 20 February at 18:00 UTC (a) and Badajoz on 21 February at 00:00 UTC (b) at different heights.Back-trajectories were calculated from the ground level to 5000 m a.s.l., at every 250 m.Lines in grey indicate trajectories arriving at the lowest levels, with no African history; in green are trajectories that passed over the southern slope of the Saharan Atlas before the observed dust mobilization; in black are the trajectories followed by the parcels residing at the times and area where dust was observed; in blue are the trajectories residing at higher levels.One representative trajectory is shown for each evolution and the altitude interval is shown in the same color as the representative trajectory.The brown line corresponds to the ground level for the trajectories more associated to the dust advection (thick black lines).The location of the air parcels around the time of observation (12:00 UTC) of the dust plumes is shown as a red circle.pressure system, which was the dominating feature in the second half of the episode.In this second case, dust was advected mostly below the cloud band and affected the eastern part of the Iberian Peninsula as well as most of the western Mediterranean Basin.
The low-pressure system weakened on early 22 February and the region was increasingly dominated by the Azores and the North African highs.As a consequence, zonal flow swept the first dust plume along northern Spain from west to east in subsequent days.The second dust plume, which was moving northward along eastern Iberia, was also displaced to the Mediterranean.The study region was then under high pressures and the event ended.

Ceilometer data analysis
The vertical structure of the dust event described above has been monitored with ICENET: firstly, by obtaining the total attenuated backscatter profiles using the calibration factors from Table (2) and, secondly, by applying the inversion and obtaining particle backscatter profiles.In addition, the integral of the backscatter profiles multiplied by the lidar ratio is used to estimate the AOD during the event and the center of mass of the backscatter profiles is considered as an indicator of the presence of a decoupled aerosol layer (a dust plume in this case) or the entrainment of the aerosol layer into the boundary layer.All these products were calculated in near-real time and serve as an example of the promising capabilities for real-time characterization of singular events with a network of distributed ceilometers.
Figure 7 shows time series of total attenuated backscatter profiles, i.e., calibrated profiles, for the five ceilometers used in this study.From top to bottom the series correspond to Granada, Tabernas, Badajoz, Valladolid and Montsec stations, respectively.Tabernas station is covered by clouds during most of the event and Montsec station is also affected by clouds during part of the event.
Dust arrives first at the stations in Granada and Tabernas (on 20 February at 18:00 UTC).As the dust plume moves northwestward we observe the dust plume in Badajoz (on 21 February at 00:00 UTC) and Valladolid (on 21 February at 06:00 UTC).At Montsec, the dust plume is detected on 21 February at 12:00 UTC).The second plume brings the cloud band and this is visible at Tabernas station around 12:00 UTC on 21 February and a bit later at 21:00 UTC on 21 February at Montsec station.Finally, the displacement of the dust from west to east at the end of the event, when the cut-off low weakens, appears as a dust plume at Valladolid on 22 February at 15:00 UTC, at Badajoz station on 22 February at 21:00 UTC, and at Granada station on 23 February at 06:00 UTC.Tabernas and Montsec are influenced by the second dust plume and the cloud band, and this is not as clearly visible as at the other stations.Another feature that is observed in Fig. 7 is that the dust plumes, especially the first one, are entrained into the boundary layer rapidly.
After applying the inversion, a quantitative comparison of stations is possible, as shown in Fig. 8, for different stages of the dust outbreak (the specific times are shown on the vertical lines in Fig. 7).The beginning of the outbreak, when the first plume arrives at the different stations, is shown in Fig. 8a.The center of mass of the dust plumes is about 3000 m a.s.l. for all stations.Additionally, for Granada and Badajoz, we observe that the height of the peak in particle backscatter coefficient is in accordance with the backward trajectory analysis shown in Sect.4.1.The arrival of the second plume is shown in Fig. 8b for all sites on 22 February at 06:00 UTC.At this stage, we observe that Granada and Tabernas stations (which are only 100 km apart) show very different behavior in particle backscatter and also in the height of the dust plume.Finally, Fig. 8c shows the final part of the outbreak when dust is mobilized eastwards to the Mediterranean Sea.In this case, dust is below 2000 m a.s.l. for Granada and Tabernas, whereas for the rest of the stations it is still observed above 3000 m a.s.l.In general, the particle backscatter coefficient profiles indicate a stronger intensity of the event at this stage of the event, after the second plume arrives, especially for Granada and Tabernas.
For the entire dust outbreak period and all stations the integral of the backscatter profiles is shown in Fig. 9a.This parameter allows identifying the beginning of the dust event for each station.Thus, an increase is observed in the integral of the backscatter in Granada around 20 February at 19:30 UTC, at Badajoz it is detected around 21 February at 05:30 UTC, and at Valladolid and Montsec it is observed at 16:30 UTC and 17:00 UTC, respectively.Due to clouds, this increase in the integral of the backscatter is not observed in Tabernas.The influence of the dust load after the first plume masks the arrival of the second plume, but the dust mobilization towards the Mediterranean sea is observed again at Badajoz (around 22 February at 20:00 UTC) and in Montsec at 23:00 UTC.The change in the integral of backscatter to larger values is coincident with the starting time observed in the total attenuated backscatter temporal series, and it is in accordance with the satellite observations and backtrajectory analysis.Additionally, the center of mass of the particle backscatter coefficient profiles is used to monitor the evolution of the profile region with more predominance of aerosol particles.Thus, in Fig. 9b for Granada before the event, the center of mass is about 1500 m a.s.l., and when the dust arrives the center of mass is elevated to 2500 m.After 9 h the center of mass is about the same as before the event, indicating that, possibly, the dust plume is no longer decoupled, and it is entrained into the boundary layer.A similar behavior is observed for Badajoz, Valladolid and Montsec stations.Again, the second plume is not observed in changes in the center of mass, but the mobilization of dust towards the Mediterranean Sea is observed as an increase in the center of mass of the profiles for Badajoz, Valladolid and Tabernas.

Conclusions
The use of ceilometers for the characterization of optical aerosol properties is possible but, due to the weak signal, it is important to screen out profiles in order to ensure the quality of the inversion.In addition, due to the vast number of data, it is important to perform all these operations in an automated, unsupervised way and, preferably, in near-real time.The methodology proposed uses ancillary data from sun photometers in order to constraint the calibration of the ceilometers.The time series of this calibration is used to determine the quality of the inversions, selecting those that present, at the reference height, a ratio of the backscattering signal to molecular attenuated backscatter within the mean calibra-tion factor ± the standard deviation.A comparison with independent lidar measurements indicates that this method allows the automatic discrimination of the quality of the inversions with ceilometers.During this comparison a difference smaller than 15 % in backscatter coefficient is observed.Thus, it is feasible to routinely provide particle backscatter coefficient profiles with ceilometers.
The inverted profiles obtained with ceilometers could be used for elevated aerosol layer alert by setting a threshold on the particle backscatter coefficient values of the profile and are potentially useful for model assimilation and evaluation since all the processing is automated and in near-real time.
This method has been applied to a group of ceilometers (ICENET) and tested during a dust outbreak reaching Spain on 20 February 2016 and lasting until 24 February 2016.This dust event affected all ICENET stations, with two distinct plumes reaching the Iberian Peninsula following different paths and a final stage where zonal flows swept the dust towards the Mediterranean Sea.This scheme of dust mobilization is unusual for this season of the year, and the intensity, spatial coverage and duration of the event make it perfect as a test for monitoring purposes with the ceilometer network.The calibration of the ceilometers allows a qualitative monitoring of the event, while the inversions provide quantitative information.Thus, ceilometers can complement lidar stations that, in principle, would operate intermittently and with less spatial density.It is worth noting that differences have been observed on profiles 100 km apart.This reinforces the need for providing vertical profiles of aerosol optical properties with a dense spatial resolution.
Parameters extracted from the particle backscatter coefficient profiles such as the integral or the center of mass can also give a quantitative idea of the presence of an elevated aerosol layer.These parameters are expected to increase with an elevated aerosol layer, and the second one can be used as a rough indicator for the deposition velocity of an elevated aerosol layer by comparing a time series of these values.
Data availability.Sun-photometer data have been downloaded from the AERONET website (http://aeronet.gsfc.nasa.gov).All raw and processed lidar and ceilometer data are available from the corresponding author upon request.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Map of the Iberian Peninsula showing the location of the ceilometers.At Granada station (circled in red) a co-located multiwavelength Raman lidar is also available.

Figure 2 .
Figure 2. Ceilometer calibration factor (C L ) vs. normalized mean bias (NMB) (a), relative difference in center of mass (b) and coefficient of correlation (R) (c).The solid horizontal line indicates the mean C L for the dust event period and dashed lines indicate the 33 % around this mean value.

Figure 3 .
Figure 3. Lidar and ceilometer particle backscatter profiles for five cases on 23 February 2016.The first case (marked in blue) is a rejected ceilometer profile and the other four cases (marked in red) are cases with a ceilometer calibration factor within the 33 % of the median calibration factor.Errors on lidar profiles have been estimated with a Monte Carlo technique(Pappalardo et al., 2004;Mattis et al., 2016) and are on the order of 10 −9 (m −1 sr −1 ); therefore, they are not shown.

Figure 4 .
Figure 4. False-color RGB dust image from MSG-SEVIRI showing different stages of the dust outbreak (dust appears magenta).

Figure 6 .
Figure 6.Sun-photometer time series representing the AOD at 1064 nm (a) and Ångström exponent between 440 and 800 nm (b) for all sites in ICENET during the dust event.

Figure 7 .
Figure 7. Ceilometer time series of total attenuated backscatter representing the evolution of the dust outbreak between 20 and 24 February 2016 (the color scale is logarithmic).Red vertical lines indicate the time of the profiles in Fig. 8: the first line of each site indicates the times for Fig. 8a, second line those for Fig. 8b, and third line those for Fig. 8c.

Figure 8 .
Figure 8. Particle backscatter coefficient profiles for all stations at the beginning (a), middle (b) and final stage (c) of the outbreak (note that the x axis has a different scale and the profiles start at ground level).The shaded areas represent the 15 % uncertainty.

Figure 9 .
Figure 9.Time series of the integral of the particle backscatter coefficient for all stations (a) and time series of the center of mass of the backscatter profiles for all stations (b).

Table 1 .
Description of the Iberian ceilometer network sites.

Table 2 .
Mean calibration factors for ceilometers in ICENET for the period 1 May 2014 to 1 May 2016.