Assimilation of lidar signals : application to aerosol forecasting in the western Mediterranean basin

This paper presents a new application of assimilating lidar signals to aerosol forecasting. It aims at investigating the impact of a ground-based lidar network on the analysis and short-term forecasts of aerosols through a case study in the Mediterranean basin. To do so, we employ a data assimilation (DA) algorithm based on the optimal interpolation method developed in the P OLAIR3D chemistry transport model (CTM) of the P OLYPHEMUS air quality modelling platform. We assimilate hourly averaged normalised range-corrected lidar signals (PR 2) retrieved from 72 h period of intensive and continuous measurements performed in July 2012 by ground-based lidar systems of the European Aerosol Research Lidar Network (EARLINET) integrated into the Aerosols, Clouds, and Trace Published by Copernicus Publications on behalf of the European Geosciences Union. 12032 Y. Wang et al.: Assimilation of lidar signals gases Research InfraStructure (ACTRIS) network and an additional system in Corsica deployed in the framework of the pre-ChArMEx (Chemistry-Aerosol Mediterranean Experiment)/TRAQA (TRAnsport à longue distance et Qualité de l’Air) campaign. This lidar campaign was dedicated to demonstrating the potential operationality of a research network like EARLINET and the potential usefulness of assimilation of lidar signals to aerosol forecasts. Particles with an aerodynamic diameter lower than 2.5 μm (PM 2.5) and those with an aerodynamic diameter higher than 2.5 μm but lower than 10 μm (PM10−2.5) are analysed separately using the lidar observations at each DA step. First, we study the spatial and temporal influences of the assimilation of lidar signals on aerosol forecasting. We conduct sensitivity studies on algorithmic parameters, e.g. the horizontal correlation length (Lh) used in the background error covariance matrix (50 km, 100 km or 200 km), the altitudes at which DA is performed (0.75–3.5 km, 1.0–3.5 km or 1.5–3.5 km a.g.l.) and the assimilation period length (12 h or 24 h). We find that DA with Lh = 100 km and assimilation from 1.0 to 3.5 km a.g.l. during a 12 h assimilation period length leads to the best scores for PM10 and PM2.5 during the forecast period with reference to available measurements from surface networks. Secondly, the aerosol simulation results without and with lidar DA using the optimal parameters ( Lh = 100 km, an assimilation altitude range from 1.0 to 3.5 km a.g.l. and a 12 h DA period) are evaluated using the level 2.0 (cloud-screened and quality-assured) aerosol optical depth (AOD) data from AERONET, and mass concentration measurements (PM 10 or PM2.5) from the French air quality (BDQA) network and the EMEP-Spain/Portugal network. The results show that the simulation with DA leads to better scores than the one without DA for PM2.5, PM10 and AOD. Additionally, the comparison of model results to evaluation data indicates that the temporal impact of assimilating lidar signals is longer than 36 h after the assimilation period.


Introduction
Aerosols consist of tiny pieces of solid or liquid matter suspended in the atmosphere.They have an impact on vegetation and human health by penetrating the respiratory system, and can lead to respiratory and cardiovascular diseases (Dockery and Pope, 1996;Lauwerys et al., 2007).They also influence visibility (Wang et al., 2009) and affect the earth's environment and climate by changing the amount of incoming solar radiation and outgoing terrestrial long-wave radiation retained in the earth's system (Intergovernment Panel on Climate Control, IPCC, 2013).Furthermore, they have an indirect effect, by changing the microphysical properties of clouds (Intergovernment Panel on Climate Control, IPCC, 2013).
In order to model the transport and formation of aerosols, a variety of chemistry transport models (CTMs) have been developed (Simpson et al., 2003;Schaap et al., 2004;Hodzic et al., 2006;Sartelet et al., 2007).In air quality modelling, CTMs are often employed to forecast aerosol concentrations.For instance, the monitoring atmospheric composition and climate (MACC, http://www.gmes-atmosphere.eu/)model and the POLYPHEMUS air quality modelling system (http:// cerea.enpc.fr/en/prevision.html)perform real-time forecasts of aerosols over Europe.However, a CTM is always a simplification of the real atmosphere, and there are large uncertainties in aerosol modelling (Roustan et al., 2010).A CTM is limited in terms of spatial and temporal resolutions.It is also limited to a restricted selection of physical and chemical processes, which are often simplified or parametrised.In addition, input data are often highly uncertain.Initial and boundary conditions of pollutants are two crucial factors in forecasting.Since initial and boundary conditions are often outputs from a larger-scale simulation, or from a fixed set of climatological average values based on long-term observations, they usually lack accuracy.On the other hand, aerosol measurements provide detailed insights into the atmosphere's current state, using satellite observations on a global scale or in situ measurements from ground-based or airborne instruments.Unfortunately, although measurements can help to improve the knowledge of the atmosphere, they do not directly provide the necessary initial or boundary conditions for aerosol modelling.
The OI method was used in several studies for improving initial conditions of CTMs.For example, it was first developed to assimilate AOD (aerosol optical depth) retrieved by satellite during the Indian Ocean Experiment (INDOEX) (Collins et al., 2001).The OI method was also used in a simplified radiative transfer model by Huneeus and Boucher (2007) to assimilate synthetic observations of MODIS (Moderate Resolution Imaging Spectroradiometer) and CALIPSO (Cloud-Aerosol Lidar with Orthogonal Polarization).Adhikary et al. (2008) assimilated monthly mean AOD data from MODIS and AERONET using the OI method to produce three-dimensional distributions of AOD over Asia.Niu et al. (2008) improved dust storm forecasts (dust concentrations) over China by assimilating satellite retrieval data and surface meteorological station data.Tombette et al. (2009) used the OI method over western Europe for assimilating PM 10 (particulate matter with an aerodynamic diameter lower than 10 µm) mass concentration observations from the BDQA (Base de Données de la Qualité de l'Air) network.The OI method was also employed in a study of inverse modelling of optical observations (lidar backscatter coefficients and AOD) by chemical DA (Kahnert, 2009).Pagowski et al. (2010) used the OI over the United States of America for DA of PM 2.5 (particulate matter of an aerodynamic diameter lower than 2.5 µm) observations.Liu et al. (2011) developed a DA system using the OI method within the National Centers for Environmental Prediction (NCEP) for assimilating MODIS AOD retrieval products (at 550 nm wavelength) from both the Terra and Aqua satellites and for analysing aerosol mass concentrations.Huneeus et al. (2012) used the OI method to estimate the emission fluxes of a range of aerosol species at a global scale by assimilating daily total and fine-mode AOD at 550 nm from MODIS into a global aerosol model of intermediate complexity.The OI method was used by Schwartz et al. (2012) individually or simultaneously to assimilate AOD at 550 nm retrieved from MODIS sensors and surface PM 2.5 observations for the analysis of aerosol mass mixing ratios at each grid point.Recently, Wang et al. (2013) used the OI within an observing system simulation experiment (OSSE) to investigate the potential impact of future ground-based lidar networks on the analysis and shortterm forecasts of PM 10 over Europe.They showed a potentially powerful impact of the future lidar networks for PM 10 forecasts.Li et al. (2013) used the OI for multiple aerosol species and for prediction of PM 2.5 in the Los Angeles basin.The OI method was also employed in a mesoscale numerical weather prediction system (GRAPES/CUACE_Dust) to study dust aerosol assimilation in eastern Asian (Wang and Niu, 2013).Jiang et al. (2013) developed a DA system in the WRF-Chem model using the OI method to explore the impact of assimilating surface observations of PM 10 over China.
The EnKF method was employed to simulate severe dust storm episodes occurring in March 2002 over China by assimilating surface dust concentration observations (Lin et al., 2008).The EnKF method was used to assimilate lidar attenuated backscatter coefficients and depolarisation ratios contained in the CALIPSO Level 1B data set (Sekiyama et al., 2010).Also, a global aerosol assimilation system was developed using the EnKF method for assimilating AOD and AAE (aerosol Ångström exponent) from the AERONET network and the MODIS satellite (Schutgens et al., 2010a, b).
4D-Var was used to assimilate the lidar network Asian dust data (Sugimoto and Uno, 2009).They showed that DA is effective for both improving the model results and estimat-ing the emission in the dust source region.Benedetti et al. (2009) also used the 4D-Var method in the European Centre for Medium-Range Weather Forecasts (ECMWF) for the Global and regional Earth-system Monitoring using Satellite (GEMS) and in situ data project, in order to issue aerosol forecasts and reanalyses of aerosol fields using AOD data from satellite sensors.
In meteorology, OI was surpassed by 4D-Var or EnKF (Kalnay, 2003), but it is still a commonly used DA method in CTMs, as OI is simple to implement and is computationally cheaper than other DA methods (Wu et al., 2008).By contrast, 4D-Var assimilates observations over a time window which could yield better results (Benedetti and Fisher, 2007) when the model is reliable.However, it is more complex to implement because the adjoint of the model is required (Benedetti et al., 2009;Sugimoto and Uno, 2009).Denby et al. (2008), Pagowski and Grell (2012) and Candiani et al. (2013) compared two different DA methods, the OI and the EnKF, for aerosol forecasts.They reported that the EnKF delivers slightly better results than the OI, but the cost of implementation of the EnKF is higher than that of the OI, due to the high number of required model simulations.The OI is then employed in this paper to assimilate observations sequentially.
Several aerosol properties have been assimilated for aerosol forecasts, e.g.surface mass concentrations (Niu et al., 2008;Tombette et al., 2009;Pagowski et al., 2010;Li et al., 2013;Wang and Niu, 2013;Jiang et al., 2013), aerosol particle number size distributions (Viskari et al., 2012), AOD data from satellites or the AERONET network (Huneeus and Boucher, 2007;Adhikary et al., 2008;Benedetti et al., 2009;Schutgens et al., 2010a, b;Liu et al., 2011;Huneeus et al., 2012;Schwartz et al., 2012), lidar backscatter coefficients (Huneeus and Boucher, 2007;Kahnert, 2009;Sekiyama et al., 2010) and lidar extinction coefficients (Campbell et al., 2010;Zhang et al., 2011).Most studies showed fast-fading DA impact on aerosol forecasting, especially in the early forecast hours (Tombette et al., 2009;Jiang et al., 2013).Wang et al. (2013) found that information on the vertical profile can extend the temporal influence of DA.However, in situ surface measurements and AOD data do not provide vertically resolved information in the atmospheric column.Lidar backscatter coefficient profiles provide information on the aerosol vertical structure, but estimating the aerosol backscatter coefficient from single-wavelength elastic lidar signals only (e.g. through the Klett-Fernald method, Klett, 1985) using an a priori value of a lidar ratio (extinction-tobackscatter ratio) brings in errors of up to 30 %.No critical assumptions are needed to calculate aerosol backscatter coefficients using the multi-wavelength aerosol lidar (e.g.Raman lidars), typically under nighttime conditions (Ansmann et al., 1992), but most operational lidar stations are singlewavelength lidars.Furthermore, a multi-wavelength aerosol lidar is more costly and mainly dedicated to scientific purposes than a single-wavelength aerosol lidar, and often per-forms at one visible light wavelength (e.g.532 nm), which is not eye safe (e.g.aviation near the city).Therefore, it is more realistic to put a single-wavelength aerosol lidar system into operational service.That is why Wang et al. (2014) developed for the first time DA algorithms to assimilate normalised range corrected lidar signals (PR 2 ) directly at one wavelength (e.g.355 nm).
This paper aims at investigating the usefulness of a ground-based lidar network in analysis and short-term forecasts of aerosols based on a case study over the Mediterranean.Important DA algorithm parameters are also studied, e.g. the correlation length in the background error covariance matrix, the altitudes at which DA is performed, and the assimilation period length.
This paper is organised as follows.Section 2 describes the modelling system, i.e. the CTM POLAIR3D/POLYPHEMUS, the OI method and the experiment design.Section 3 provides a description of the observations used.DA parameter tests are conducted in Sect. 4. Results are shown and discussed in Sect. 5. Our findings are summarised in Sect.6.
The aerosol optical property module developed by Wang et al. (2014) is employed.It simulates the molecular backscatter and extinction coefficients (β m and α m ) from the Boltzmann constant, the atmospheric pressure, and temperature.The aerosol extinction and backscatter coefficients (β a and α a ) are simulated from the model aerosol concentration outputs (i.e.aerosol water content and aerosol) by estimating the particle wet diameter and the aerosol complex refractive index of a particle.Lidar signals (i.e.PR 2 normalised at a reference altitude) and AOD are simulated as functions of the molecular backscatter and extinction coefficients and the aerosol extinction and backscatter coefficients.
The modelling domain is the same as the one used for the forecasts at http://cerea.enpc.fr/en/prevision.html.It covers western Europe and parts of eastern Europe (15 • W, 35 • E × 35 • N, 70 • N, see Fig. 1), with a horizontal resolution of 0.5 • × 0.5 • .In the simulation, 14 vertical levels are considered from the ground to an altitude of 12 000 m a.g.l.
(above ground level).The heights of the cell interfaces are 0, 30,60,100,150,200,300,500,750,1000,1500,2400,3500,6000 and 12 000 m a.g.l.Meteorological inputs are interpolated from reanalysis provided every 3 h by the European Centre for Medium-Range Weather Forecasts (ECMWF).Boundary conditions are climatological conditions obtained from averaging boundary conditions from MOZART4 (Model for OZone And Related chemical Tracers version 4) (Emmons et al., 2010) over the years 2004-2008.Sea-salt emissions are assumed to be made up of 39.33 % of sodium, 55.025 % of chloride and 7.68 % of sulfate, and modelled following Monahan et al. (1986).Anthropogenic emissions of gases and aerosols are generated with the EMEP inventory for 2009.For example, the EMEP provides yearly emissions of PM 2.5 and coarse PM (PM with an aerodynamic diameter higher than 2.5 µm but lower than 10 µm).The PM 2.5 fraction is speciated into mineral dust, black carbon and primary organic aerosol.The PM coarse fraction is attributed to mineral dust.In the simulation, Saharan dust is only forced by boundary conditions.
The OI approach is used for assimilating lidar signals from the model aerosol concentration outputs (Wang et al., 2014).The analysed mass concentration x a is obtained from the equation where x b are the model mass concentrations, y is the observation vector, H is the observation operator that simulates normalised PR 2 from the mass concentrations x b , H is the tangent linear operator of H , and B and R are respectively the background and observation error covariance matrices.Wang et al. (2014) provided two algorithms based on the OI method to compute the analysed state x a .One algorithm analyses PM 10 concentrations.The other analyses PM 2.5 and PM 10−2.5 concentrations separately but simultaneously.Wang et al. (2014) reported that the latter algorithm leads to better forecasts than the former, because the model often simulates PM 2.5 better than PM 10−2.5 , and the background error variances of PM 2.5 and PM 10−2.5 are set separately in the latter algorithm.Therefore, we employ the latter algorithm in this paper.We set the background error covariance matrix as a block diagonal matrix having two main diagonal blocks.One main diagonal block is set as the background error variance matrix of PM 2.5 .Another is set as the background error variance matrix of PM 10−2.5 .We set the background error of PM 2.5 and PM 10−2.5 at 5 µg m −3 and 30 µg m −3 respectively in B, since the model simulates PM 2.5 more accurately than PM 10−2.5 (see Sect. 5).We take  R = σ 2 I, where σ is an observation standard deviation (depending on instrumental and representativeness error variances) and I is the identity matrix in the observation space.
The value of σ is different in each DA test.It is determined by a χ 2 diagnosis, in which the scalar χ 2 at each DA step is defined by On average, χ 2 should be equal to the number of observations (Ménard et al., 1999).This χ 2 diagnosis could balance observation and background errors.After DA, the analysed concentrations are redistributed over the model variables following the initial (background) chemical and size distributions (Tombette et al., 2009;Wang et al., 2013Wang et al., , 2014)).
The simulation with DA, referred to as the DA experiment, consists of two periods: an assimilation period and a forecast period.During the assimilation period, at each time step, all available lidar data retrieved in the framework of the pre-ChArMEx (Chemistry-Aerosol Mediterranean Experiment)/TRAQA (TRAnsport à longue distance et Qual-ité de l'Air) and ACTRIS/EARLINET campaigns are assimilated.During the forecast period, DA is not performed.Hence, the model mass concentrations evolve depending on initial and boundary conditions, emissions and meteorology.Concentrations can be impacted by lidar DA far from the place where lidar signals are assimilated, because analysed mass concentrations are transported by winds and diffusion (turbulence).
In regional models, uncertainties are linked to input data and parametrisations, e.g.initial and boundary conditions (Roustan et al., 2010), meteorological inputs (Dawson et al., 2007) and emissions (de Meij et al., 2006;Napelenok et al., 2006).Replacing some input data, such as boundary conditions or emissions, with another set of data which is also uncertain may either improve or deteriorate the aerosol simulations locally, depending on period/year and place.However, DA may be used to improve input data, such as initial conditions, using observations (as done in this paper).The impact of DA may vary locally with the quality of the input data used.
Table 1.Description of the lidar systems used in this study."Reso."stands for resolution."ASL" stands for a.s.l.(altitude above sea level).The letters "p" and "c" in the wavelengths stand for parallel and cross linear polarisation component respectively.Table 2. DA tests with different configurations for the evaluation of the impact of the assimilation parameters on the forecasts.L h is the horizontal correlation length used in the Balgovind approach (Balgovind et al., 1983) to define the error covariance matrix B.

Observations
In the following, we describe the observations used in this study: the lidar signals used for assimilation, and surface mass concentrations and AOD used for the DA validation.

Lidar observations
An intensive measurement effort was performed by 12 ground-based lidar stations from the ACTRIS/EARLINET network (Bösenberg et al., 2003;Pappalardo et al., 2014) in the Mediterranean basin and one station in Corsica in the framework of the pre-ChArMEx/TRAQA and AC-TRIS/EARLINET campaigns in July 2012 during 72 h.All stations were located on the northern side of the Mediterranean.One of the goals of this campaign was to locate and track aerosols in the lower and middle troposphere in the Mediterranean region and to help improve our forecast ability of CTMs using DA.The ground-based lidar stations (blue/grey and yellow star markers in Fig. 1) performed continuous measurements from 9 July at 06:00 UTC until 12 July at 06:00 UTC.The participating EARLINET stations include Athens, Barcelona, Bucharest, Evora, Granada, L'Aquila, Limassol, Madrid, Messinia, Potenza, Payerne and Clermont-Ferrand.The MISTRALS/ChArMEx station was situated at INRA (Institut National de la Recherche Agronomique), San Giuliano, at about 3 km from the eastern coastline of Corsica (see Fig. 1).Data received by the Payerne and Messinia stations are not available.Also, data received by the Limassol station are not used in this paper, because Limassol is outside of the model domain.
Table 1 shows the site coordinates and properties of the lidar systems used in this campaign.The vertical resolution of measurements ranges from 3.25 m to 30 m (depending on the lidar system).The temporal resolution of measurements ranges from 30 s to 300 s (depending on the lidar system).The raw data (except those of the Corsica station) were automatically treated by the single calculus chain (SCC) developed by the EARLINET lidar network (http: //www.earlinetasos.org)(D'Amico et al., 2012) to generate integrated profiles of range-corrected lidar signals (PR 2 ) in near real time under cloud-free conditions.The SCC is an automatic tool to get different types of aerosol products (e.g.PR 2 , aerosol extinction and backscatter coefficients) from raw lidar data.In this work, only one type of the available products, PR 2 , is used.All observations are integrated with a time resolution of 1 h, the DA time step used in this study, and normalised at an altitude in the range of the molecular zone.It is because there is almost no aerosol in the molecu-lar zone.The linear approximation of the observed lidar signal should be equal to the one of the simulated molecular signal (without aerosol contribution) in the molecular zone (Wang et al., 2014).In this paper, it is taken at 4750 m a.g.l., which corresponds to the model level of 3500-6000 km a.g.l.

Observations for validation
We employ two independent data types for DA validation: surface mass concentration measurements (i.e.PM 2.5 and PM 10 ) and AOD data.
The surface mass concentration measurements are from the BDQA (Base de Données sur la Qualité de l'Air, the French national database for air quality which covers France) network, the Barcelona network (three stations), the EMEP-Spain/Portugal network, and the EMEP-Europe database (see Fig. 1).The French and Barcelona networks (triangles in Fig. 1) provide hourly averaged mass concentration measurements of PM 2.5 and PM 10 .The EMEP-Spain/Portugal and EMEP-Europe networks (squares in Fig. 1) provide daily averaged mass concentration measurements of PM 10 .The number of used stations in the BDQA, Barcelona, EMEP-Europe and EMEP-Spain/Portugal networks, which provide PM 10 or PM 2.5 measurements in July 2012, is shown in Table 3.The BDQA network provides the most measurements, with 240 stations for PM 10 and 70 stations for PM 2.5 .
The hourly AOD data at 355 nm are derived by level 2.0 (cloud-screened and quality-assured) AOD data at 340 and 380 nm retrieved from AERONET (AErosol RObotic NETwork, http://aeronet.gsfc.nasa.gov/)following the Ångström law (Wang et al., 2014).The locations of the AERONET stations considered (e.g.stations that are close to the lidar network and that provide the level 2.0 AOD data in the pre-ChArMEx/TRAQA and ACTRIS/EARLINET campaign) are shown as orange diamonds in Fig. 1.Thirteen AERONET stations are used for validation in this paper (see Table 3).

Case study
The Mediterranean basin is the receptacle of aerosols from different origins, e.g.biogenic emissions, natural dust emissions from the Sahara (Moulin et al., 1998;Hamonou et al., 1999), anthropogenic emissions from highly populated coastal areas, marine aerosols, and wildfires.Emissions from anthropogenic and biogenic origins strongly interact to form secondary organic aerosols (Sartelet et al., 2012).The aerosol load is often high over the Mediterranean region (Putaud et al., 2010;Nabat et al., 2013).Therefore, it is a good place to test the usefulness of lidar DA to improve the forecast of CTMs.
Figure 4 shows wind fields at about 2 km a.g.l.interpolated from ECMWF data for each morning of the lidar measurement period, i.e. 9 July 2012 at 08:00 UTC, 10 July 2012 at 08:00 UTC, 11 July 2012 at 08:00 UTC, and 12 July 2012 at 08:00 UTC.Westerly or northerly winds transported pollution over the Mediterranean during the lidar campaign.Figure 5 shows the AODs at 550 nm retrieved from MSG (Meteosat Second Generation)/SEVIRI satellites (http://www.icare.univ-lille1.fr/msg/,Thieuleux et al., 2005) 15 min image averaged from all available images between 04:00 and 18:00 UTC on 9-12 July 2012, where the high AODs observed mainly in the southern part of the Mediterranean were mostly due to Saharan dust.However, penetration of the Saharan dust plume over the continent of Europe was limited, except in the south of Italy and the south and east of Spain.At the Ersa surface station in Corsica, the chemical analysis of filters from 00:00 until 12:00 UTC on 11 July 2012 did not detect Saharan dust (Nicolas, 2013), and the MIS-TRALS/ChArMEx aerosol lidar in Corsica does not show evidence of a dust layer (see Fig. 3).
To check that the penetration of the Saharan dust plume over the continent of Europe was limited, and to assess where analysed concentrations are transported to after assimilation, Fig. 6 shows 48 h backward trajectories (dashed lines) of air masses arriving at 2 km a.g.l. and 72 h forward trajectories (solid lines) of air masses departing at 2 km a.g.l. at 10 lidar stations.These data are outputs of the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT; http://ready.arl.noaa.gov/HYSPLIT.php)model (Draxler and Rolph, 2014;Rolph, 2014) using the Global Data Assimilation System (GDAS) meteorological data with a resolution of 1 • × 1 • .These backward (forward) trajectories end (start) respectively at 06:00 UTC on 9 July 2012, 10 July 2012, 11 July 2012 and 12 July 2012.They indicate that aerosols measured in Spain, Portugal and France were transported to the northeast or east.Aerosols measured by lidars at other stations (i.e.Athens, L'Aquila, Potenza, and Bucharest) were transported to the south or east.Those observations are in agreement with wind fields shown in Fig. 4. In addition, there is almost no rainfall along trajectories (not shown in this paper).

Assimilation parameter tests
In this section, we perform sensitivity tests, first on the DA period length, and then on the horizontal correlation length used in the background error covariance matrix and on the assimilation altitude range.  1 "Obs."stands for observation. 2"Sim."stands for simulation. 3"Corr."stands for correlation. 4Correlation is not presentable for three stations.

Assimilation period length
to 3 days, during which surface mass concentration observations were assimilated.They suggested that an assimilation period of 12 h is necessary to improve the aerosol forecast.In this work, two DA period lengths, 12 h and 24 h, are employed to study the impact of the assimilation period length on aerosol forecasting.The results are detailed in this section.
The 72 h period of continuous lidar measurements from 06:00 UTC 9 July to 06:00 UTC 12 July 2012 is split into three experiments of 24 h each.For the assimilation period of 12 (or 24) hours, for each of the three experiments, the lidar data are assimilated during 12 (or 24) hours, and 60 h forecasts are issued at 06:00 UTC on 10, 11 and 12 July, respectively.All DA experiments use the same parameters (i.e. the horizontal correlation length is 100 km and the assimilation altitude ranges from 1.0 to 3.5 km a.g.l.), except for the assimilation period length.
Figure 7 shows the scores, the RMSE (root mean square error) and the (Pearson) correlation calculated against the ground observations over France (the BDQA network) for PM 10 and PM 2.5 , since the BDQA network provides most measurements of PM 10 and PM 2.5 .We refer to Appendix A for the definition of statistical indicators.Overall, the simulation with lidar DA leads to better scores than the simulation without DA during the first 36 h of forecast.The improvements in DA are significant for PM 10 .The RMSE (or correlation) of PM 10 averaged over the first 36 h of forecast is 9.4 µg m −3 (or 39 %) without DA, 8.4 µg m −3 (or 49 %) with 12 h DA and 8.4 µg m −3 (or 50 %) with 24 h DA.For PM 2.5 , the improvements in DA are not significant, except for the correlation.The RMSE (or correlation) of PM 2.5 averaged over the first 36 h of forecast is 4.5 µg m −3 (or 37 %) without DA and 4.4 µg m −3 (or 43 %) with either 12 h DA or with 24 h DA.Comparing DA with 24 h of analysis (DA test "24 h DA") to 12 h of analysis (DA test "12 h DA"), the simulation with 24 h of analysis delivers slightly better scores during the forecast period (to the right of the black lines).However, the difference between DA tests "24 h DA" and "12 h DA" after 6 h forecasts is barely significant.Since the measurement period of the lidar campaign in July 2012 lasted only 72 h, and the simulation with 24 h of analysis does not lead to much better scores than the one with 12 h of analysis during the forecast period, we chose to perform DA experiments with an assimilation period of 12 h in the following to have sufficient DA experiments to evaluate the results of DA statistically.In this case, the 72 h period of continuous lidar measurements is split into six 12 h assimilation periods (should 24 h DA be chosen, the 72 h period of continuous lidar measurements would be split into only three disjoint 24 h assimilation periods).Figure 8 shows the schematic representation of these six DA experiments.Each DA experiment includes a 12 h assimilation period (grey bars) and a 60 h forecast period (white bars).All DA experiments begin either at 06:00 UTC or at 18:00 UTC on 9, 10 or 11 July 2012.Figure 9 shows the schematic representation of the lidar measurement segments assimilated in six DA experiments.At each DA step, all available lidar data retrieved from 10 lidar stations are assimilated.

Assimilation correlation length
In Table 2, the different configurations of DA are summarised, with the horizontal correlation length L h (e.g.50, 100 and 200 km) and the assimilation altitude range used.The scores (i.e.RMSE and correlation) of the different con-figurations of DA for PM 10 and PM 2.5 are shown in Fig. 10.These scores are calculated against the observations of the BDQA network.In this section, the impact of the horizontal correlation length L h of the error covariance matrix B is studied, since L h is an important parameter that determines to what horizontal extent the particle concentrations are corrected by DA.
At the beginning of the assimilation period, all simulations have the same scores, since the simulations without DA and with DA use the same initial conditions.The improvement in aerosol mass concentrations at stations over France starts 6 h after the start of the DA experiment.This delay is due to the fact that the only lidar station in France used for this study is in Corsica, away from continental France (see Figs. 4 and 6; the station in Clermont-Ferrand provided too few observations due to bad weather during the campaign, see Fig. 9).It is also because the assimilation altitude range is high: it starts higher than 1.0 km a.g.l.The analysed mass concentrations therefore take time to be transported to the ground level.We find that the correlation length L h = 200 km (yellow lines in Fig. 10) is too large, decreasing slightly the correlation coefficients for both PM 10 and PM 2.5 at French stations during the assimilation period (to the left of the black lines in Fig. 10).The scores are computed for the BDQA network (hourly data).The vertical black lines denote the separation between the assimilation period (to the left of the black lines) and the forecast (to the right of the black lines)."12 DA" ("24 DA") stands for DA with 12 (24) hours of analysis.The forecasts of "12 DA" and "24 DA" start at the same moment.The scores in the first 12 analysis hours of "24 DA" are not shown.
During the forecast period (to the right of the black lines in Fig. 10), the temporal impact of the assimilation of lidar signals lasts longer than 36 h for all DA tests.Notice that the temporal impact of assimilating surface mass concentrations was estimated to be between 6 and 12 h (Tombette et al., 2009;Jiang et al., 2013).The comparison of the DA tests with L h = 50 km (green lines in Fig. 10), L h = 100 km (red lines in Fig. 10) and L h = 200 km (yellow lines in Fig. 10) shows that using L h = 100 km leads to better forecasts than using L h = 50 km or L h = 200 km on the first forecast day.In addition, using L h = 200 km (yellow lines in Fig. 10) results in a higher RMSE than the simulation without DA for PM 2.5 on the first forecast day.This is because the analysed zone in the model is set to be isotropic (a horizontal disc, the centre of the disc being the measurement station, i.e. the lidar site), whereas the analysed zone should be horizontally anisotropic, depending on the wind direction and the aerosol spatial distribution (e.g.aerosol origins).Using L h = 200 km defines an isotropic analysed zone which is too large, leading to a decrease in correlation coefficients.On the second forecast day, using L h = 200 km leads to much better scores than using L h = 50 km or L h = 100 km for both PM 10 and PM 2.5 .Moreover, the beneficial impact of the assimilation with L h = 200 km lasts longer than 48 h.It is because using L h = 200 km leads to higher corrections around the lidar site due to the use of the Balgovind approach (Balgovind et al., 1983) (the closer to the lidar site the grid point is, the higher the correlation is).The corrections due to the higher correlation around lidar sites (far away from France) are more accurate, and impact France later.

Assimilation altitude range
The choice of the assimilation altitude range is influenced by two factors.First, as the normalisation of range-corrected lidar signals is done at high altitudes, the lower the altitude is, the higher the error in the simulated lidar signal is.It is mostly because the integration of simulated extinction coefficients from the considered altitude to the normalisation altitude leads to accumulation of errors of simulated lidar signals at high altitudes, especially in the case where highaltitude aerosol layers are not well modelled (Wang et al., 2014).Second, the numerical computations of the lidar operator H and its tangent lidar operator H (see Eq. 1) are very costly.The larger the assimilation altitude range is, the more costly the numerical computation is.Hence, in this section, we investigate the impact of the assimilation altitude range on aerosol forecasting.We perform three DA tests (0.75-3.5 km a.g.l., 1.5-3.5 km a.g.l. and a reference case 1.0-3.5 km a.g.l. in Fig. 10).As shown in Fig. 10, assimilating lidar signals from 0.75 to 3.5 km a.g.l.(magenta lines) leads to similar scores (with respect to hourly data of BDQA) as assimilating from 1.0 to 3.5 km a.g.l.(a reference case, red lines).A first explanation is that the observation variance (sum of instrumental and representativeness variances, from the χ 2 diagnosis) of the model level of 0.75-1.0km a.g.l. is set higher than those of the model levels from 1.0 to 3.5 km a.g.l.A second ex-planation is that the scores in Fig. 10 are computed using the observations from the BDQA network, where most improvements are from assimilation of lidar signals in Spain or Portugal (see Figs. 4 and 6).However, of the lidar stations in Spain, only Madrid and Granada provided data between 0.75 and 1.0 km a.g.l.(see Figs. 2 and 3).In addition, assimilating lidar signals from 1.0 to 3.5 km a.g.l.(magenta lines) leads to slightly better scores than from 1.5 to 3.5 km a.g.l.(black lines).) and the correlation of PM 10 (or PM 2.5 ) averaged for each of the six different experiments.The scores are computed for the BDQA network (hourly data).The vertical black lines denote the separation between the 12 h assimilation period (to the left of the black lines) and the 60 h forecast period (to the right of the black lines).The "DA L h = 50 km", "DA L h = 100 km" and "DA L h = 200 km" simulations correspond to an assimilation altitude range from 1.0 to 3.5 km.The "DA 0.75-3.5 km" and "DA 1.5-3.5 km" simulations correspond to L h = 100 km.

Results and discussions
The US Environmental Protection Agency (EPA) has issued minimal guidance on PM model performance evaluation measures, goals, and criteria.Boylan and Russell (2006) suggested using the mean fractional bias (MFB, %) and the mean fractional error (MFE, %), because they bound the maximum bias and error (see Appendix A).We evaluate the simulation without DA using the hourly observations from the French BDQA network with these criteria.For PM 10 (or PM 2.5 ), the MFB and MFE averaged over all six experiments are respectively −29 % and 46 % (or 6 % and 43 %).For both PM 10 ) of PM 10 (or PM 2.5 ) averaged over the different experiments without and with DA (L h = 100 km and altitude range 1.0-3.5 km).For the six successive experiments, the time origin corresponds respectively to 06:00 UTC on 9 July, 18:00 UTC on 9 July, 06:00 UTC on 10 July, 18:00 UTC on 10 July, 06:00 UTC on 11 July and 18:00 UTC on 11 July.The scores are computed for three stations around Barcelona (hourly data, see Fig. 1).The vertical black lines denote the separation between the 12 h assimilation period (to the left of the black lines) and the 60 h forecast period (to the right of the black lines).and PM 2.5 , the criteria evaluation goals are verified.However, the model simulates PM 2.5 better than PM 10 , which is slightly underestimated.This is probably because road resuspension of PM is not considered, either in the model or in the input data (e.g.boundary conditions).As a consequence, we have set a lower standard deviation for PM 2.5 (i.e. 5 µg m −3 ) than for PM 10−2.5 (i.e. 30 µg m −3 ) in the background error covariance matrix B (see Eq. 1).As discussed in Sect.4, the "DA L h = 100 km" DA test, which assimilates lidar signals retrieved from the lidar campaign from 1.0 to 3.5 km a.g.l.during 12 h with L h = 100 km, and which performs 60 h forecasts, delivers the best scores during the forecast period.Therefore, in the following, we consider the "DA L h = 100 km" DA test ("Lidar DA" hereafter).Since most improvements are in the first 36 h of forecast, we compute the scores only for this period hereafter, instead of for the whole forecast period (i.e. 60 h).

Validation with the BDQA network
For PM 10 , the averaged RMSE (correlation) over the first 36 h of forecast is 8.8 µg m −3 (40 %) without DA and 8.0 µg m −3 (49 %) with DA (see Table 4).For PM 2.5 , the averaged RMSE (correlation) over the first 36 h of forecast is 4.4 µg m −3 (39 %) without DA and 4.3 µg m −3 (44 %) with DA (see Table 4).Notice that DA improves PM 10 more efficiently than PM 2.5 .Therefore, DA would be very useful in reducing the uncertainties in the simulation due to road resuspension of coarse PM.However, these improvements are not Against the observations at BDQA stations on the southern side of 44 • N (dashed line in Fig. 1), the averaged RMSE (MFB and MFE respectively) of PM 10 over the first 36 h of forecast is 16.4 µg m −3 (−53 % and 30 % respectively) without DA and 13.7 µg m −3 (−26 % and 46 % respectively) with DA.For PM 2.5 , the averaged RMSE (MFB and MFE respectively) over the first 36 h of forecast is 7.1 µg m −3 (−20 % and 47 % respectively) without DA and 6.5 µg m −3 (−6 % and 44 % respectively) with DA.The improvements in DA are more significant by comparisons to measurements at BDQA stations south of 44 • N than at all the BDQA stations.Aerosol forecasts over these southern stations are impacted by DA of the Corsica, Spain and Portugal lidar data (see Fig. 6).
Moreover, we compare simulations with DA in the daytime (DA is performed from 06:00 to 18:00 UTC) to simulations with DA in the nighttime (DA is performed from 18:00 to 06:00 UTC).We find that they lead to similar scores (results not shown in this paper).

Validation with the Barcelona network
Figure 11 shows the time evolution of the RMSE averaged over all six experiments without and with lidar DA for PM 10 and PM 2.5 .The RMSEs are computed at three surface stations around the Barcelona lidar station (violet triangles in Fig. 1).We find that the impact of the assimilation of lidar signals is longer than 48 h in the forecast period for both PM 10 and PM 2.5 .The averaged RMSE of PM 10 over the first 36 h of forecast is 8.9 µg m −3 without DA and 7.0 µg m −3 with DA (see Table 4).The averaged RMSE of PM 2.5 over the first 36 h of forecast is 6.0 µg m −3 without DA and 4.7 µg m −3 with DA (see Table 4).We find that the aerosol error reduction around Barcelona is higher than the one over France and the south of France (estimated using the BDQA network).That is because the surface stations around Barcelona are close to the ground-based Barcelona lidar station, leading to larger benefits of DA.Furthermore, the surface stations around Barcelona are also strongly influenced by the Evora and Madrid lidar sites due to wind fields, because Barcelona is on the leeward side of these lidar sites during the lidar campaign in July 2012 (see Fig. 6).The improvements due to lidar DA associated with long-range transport of pollution from Evora and Madrid are also validated.

Validation with the EMEP-Spain/Portugal network
Figure 12 shows the scatter plots of simulated PM 10 concentrations without and with DA against PM 10 daily measurements at EMEP-Spain/Portugal stations (cyan squares in Fig. 1).The PM 10 correlation and RMSE are slightly im-proved.During the assimilation and forecast periods (72 h), the RMSE averaged over all six experiments is 6.9 µg m −3 without DA and 6.3 µg m −3 with DA (see Table 4).Compared to the simulations without DA, DA ("Lidar DA") increases the correlation from 58 % to 63 % (see Table 4).Meanwhile, the mean bias error (MBE) decreases from 3.1 to 2.3 µg m −3 (see Fig. 12).Also, we compute the statistics of the simulation results without and with DA using daily concentrations at all EMEP-Europe stations (seven stations, green squares in Fig. 1).However, since EMEP-Europe stations are far away from the lidar network, the PM 10 RMSE, correlation and bias are slightly but barely improved (not shown).

Validation with the AERONET network
Figure 13 shows the time evolution of the AOD measurements and AODs of the 36 h forecasts without DA and with DA at AERONET stations Rome (41.84 • N, 12.65 • E, 130 m a.g.l.) and Bucharest (44.35 • N, 26.03 • E, 93 m a.g.l.).The impact of assimilating lidar signals lasts about 36 h, which corresponds to the findings of Sects.5.1 and 5.2.
Figure 14 shows the scatter plots of simulated AODs without and with DA against AOD from hourly measurements of the AERONET network over the first 36 h of forecast, where only 13 AERONET stations leeward and close to the lidar network are considered (see Fig. 1).As shown by comparing the left panels of Figs. 12 and 14, the model simulates AOD better than PM 10 .This is mostly because the model simulates fine particles (PM 2.5 ) better over the modelling domain (horizontal and vertical), and PM 2.5 tends to have larger contributions to optical properties than coarse particles when no Saharan dust event occurs (Chazette et al., 2005;Randriamiarisoa et al., 2006).This is also probably because the model may simulate the integrated mass concentration better than vertically resolved mass concentrations.
As shown in Fig. 14, AODs are significantly improved in the simulation with DA for high AOD observations (few cases).When the observed AODs are larger than 0.4 (N = 262), the RMSE (MBE) is 0.23 (0.2) without DA against 0.20 (0.13) with lidar DA.It is because large AODs could be associated with transport of particles above the boundary layer, which is not well simulated by the model (probably due to large-scale model uncertainties), but which is followed by the lidars (Wang et al., 2014).It may also be that assimilation of lidar signals improves the estimation of aerosol mass concentrations more efficiently when aerosol concentrations are high, e.g. during air pollution events, that is, when the lidar signal is strong.

Conclusions
In this paper, a data assimilation (DA) algorithm based on the optimal interpolation (OI) method is used to assimilate lidar signals (normalised PR 2 ) for aerosol forecasts over Europe.The lidar data were retrieved from a 72 h period of intensive and continuous measurements in July 2012.The measurements were performed by 12 ground-based lidar stations of ACTRIS/EARLINET in the Mediterranean basin and one station in Corsica which was set up in the framework of the pre-ChArMEx/TRAQA campaign.
First, we studied the impact of the length of the lidar DA period on aerosol forecasts.We found that 24 h DA leads to slightly better forecasts than 12 h DA.However, the differences between 24 h DA and 12 h DA are small after 6 h of forecast.Furthermore, because the impact of lidar DA lasts longer than 36 h in the forecast period, we have used 12 h as the assimilation period length in this paper.Also, we conducted sensitivity studies on algorithmic parameters, e.g. the horizontal error correlation length and altitudes at which DA is performed.DA with the error correlation length L h = 100 km and assimilation from 1.0 to 3.5 km a.g.l.leads to the best scores for PM 10 and PM 2.5 during the forecast period (the evaluation was done using measurements from the BDQA network, since the BDQA network provides the most measurements for the DA validation).
The simulation results without and with lidar DA were evaluated using hourly concentration measurements from the BDQA network over France, daily concentration measure-ments from the EMEP-Spain/Portugal network and AOD measurements from the AERONET network over Europe.The results showed that the simulation with DA leads to better scores than the one without DA for aerosol forecasts (PM 2.5 , PM 10 and AOD).Moreover, the temporal impact of assimilating lidar signals is longer than 36 h, whereas this temporal impact was estimated to be shorter in previous studies that assimilated surface mass concentrations, e.g. between 6 and 12 h, by Tombette et al. (2009) and Jiang et al. (2013).When the temporal impact was estimated using only the three stations around the Barcelona lidar site, the impact lasted longer than 48 h.Additionally, since the model simulates fine particles better than coarse particles, we set a higher error in the background error covariance matrix (see Sects. 2 and 5) for coarse particles than for fine particles, leading to larger corrections by DA of coarse particle concentrations than of fine particle concentrations.
The maximum likelihood ensemble filter (MLEF) (Zupanski, 2005) or the iterative ensemble Kalman filter (IEnKF) (Bocquet andSakov, 2013a, 2014) could be used in forecasts of aerosols in place of the OI method in order to avoid the tangent linear approximation of the lidar observation operator, and would handle the nonlinearity of the lidar observation operator.They would also update and propagate the background error covariance matrix during the assimilation period.As some lidars provide measurements at several channels, we expect to have better results by assimilating a more extended lidar data set, i.e.PR 2 , at several channels.
The lidar-derived PBL height (Morille et al., 2007;Baars et al., 2008;Granados-Muñoz et al., 2012;Lesouëf et al., 2013) could also be assimilated in the model.More accurate PBL heights would improve the forecast ability of air quality models (Pielke and Uliasz, 1998), because the PBL height determines the volume in which pollutants are mixed.Finally, as Schwartz et al. (2012) have shown, simultaneous DA of different aerosol observations (PM 2.5 and AOD) produced the best overall forecasts; for future works, we may combine DA of lidar signals and mass concentration measurements for real-time forecasts of aerosols.

Figure 1 .
Figure 1.Locations of the different measurement sites used in this paper (see Tables 1 and 3 for the number of stations used in the different networks).The rectangular area delimited by the black box shows our modelling domain.The red triangles indicate the locations of the stations of the French air quality network (BDQA).The cyan squares indicate the locations of the stations of the EMEP-Spain/Portugal network.The violet triangles indicate the locations of the stations around Barcelona.The green squares indicate the locations of the EMEP-Europe stations.The orange diamonds indicate the locations of the AERONET stations.The dark blue/grey star markers indicate the locations of ACTRIS/EARLINET stations.The grey star markers indicate lidar stations without data between 9 and 12 July or outside of the forecast domain.The yellow star marker indicates the location of the Corsica lidar station.The dashed line shows the latitude of 44 • N which is used to split the French stations in Sect.5.1.

Figure 6 .Figure 6 .
Figure 6.Backward (resp.forward) trajectories of 48 h (resp.72 h) at lidar site locations (black stars) at 2 km a.g.l.ending (resp.starting) at 06:00 UTC 9 July 2012, 06:00 UTC 10 July 2012, 06:00 UTC 11 July 2012 and 06:00 UTC 12 July 2012.Data are from the HYSPLIT Model.Dashed (resp.solid) lines indicate backward (resp.forward) trajectories, where the 12 h spacing is given by the discs.The backward trajectories pertain to the source attribution problem of the lidar measurements whereas the forward trajectories show the propagation of the DA updates around lidar locations.

Figure 7 .
Figure 7.The top (or bottom) panel shows the time evolution of the RMSE (µg m −3 ) and the correlation of PM 10 (or PM 2.5 ) averaged over the different DA experiments for three experiment types: one without DA, one with 12 h of DA and one with 24 h of DA.The scores are computed for the BDQA network (hourly data).The vertical black lines denote the separation between the assimilation period (to the left of the black lines) and the forecast (to the right of the black lines)."12 DA" ("24 DA") stands for DA with 12 (24) hours of analysis.The forecasts of "12 DA" and "24 DA" start at the same moment.The scores in the first 12 analysis hours of "24 DA" are not shown.

Figure 9 .
Figure 9. Schematic representation of the lidar measurement segments assimilated (black segments) during the assimilation period for six DA experiments."Cler.-Ferr."stands for Clermont-Ferrand.

Figure 10 .
Figure10.The top (or bottom) panel shows the time evolution of the RMSE (µg m −3 ) and the correlation of PM 10 (or PM 2.5 ) averaged for each of the six different experiments.The scores are computed for the BDQA network (hourly data).The vertical black lines denote the separation between the 12 h assimilation period (to the left of the black lines) and the 60 h forecast period (to the right of the black lines).The "DA L h = 50 km", "DA L h = 100 km" and "DA L h = 200 km" simulations correspond to an assimilation altitude range from 1.0 to 3.5 km.The "DA 0.75-3.5 km" and "DA 1.5-3.5 km" simulations correspond to L h = 100 km.

Figure 11 .
Figure11.The top (or bottom) panel shows the time of the RMSE (µg m −3 ) of PM 10 (or PM 2.5 ) averaged over the different experiments without and with DA (L h = 100 km and altitude range 1.0-3.5 km).For the six successive experiments, the time origin corresponds respectively to 06:00 UTC on 9 July, 18:00 UTC on 9 July, 06:00 UTC on 10 July, 18:00 UTC on 10 July, 06:00 UTC on 11 July and 18:00 UTC on 11 July.The scores are computed for three stations around Barcelona (hourly data, see Fig.1).The vertical black lines denote the separation between the 12 h assimilation period (to the left of the black lines) and the 60 h forecast period (to the right of the black lines).

Figure 12 .Figure 12 .Figure 13 .Figure 13 .
Figure 12.Scatter plots of simulated PM 10 mass concentrations without DA (left panel) and with DA (right panel) against daily PM 10 measurements at several EMEP-Spain/Portugal stations.

Figure 14 .Figure 14 .
Figure 14.Scatter plots of simulated AODs without DA (left panel) and with DA (right panel) against AOD hourly measurements at different AERONET stations over the first 36 h of forecast.

Table 3 .
Number of stations used for PM 2.5 , PM 10 or AOD in the different networks.

Table 4 .
Statistics (see Appendix A) of the simulation results for the different networks.

www.atmos-chem-phys.net/14/12031/2014/ Atmos. Chem. Phys., 14, 12031-12053, 2014 12046 Y. Wang et al.: Assimilation of lidar signals very
significant, especially for PM 2.5 , because most BDQA stations are far away from the lidar network.Improvements are more significant when using only stations in the south of France.