Aerosol forecast over the Mediterranean area during July 2013 ( ADRIMED / CHARMEX )

The ADRIMED (Aerosol Direct Radiative Impact on the regional climate in the MEDiterranean region) project was dedicated to study the atmospheric composition during the summer 2013 in the European Mediterranean region. During its campaign experiment part, the WRF (Weather Research and Forecast Model) and CHIMERE models were used in the forecast mode in order to decide whether intensive observation periods should be triggered. Each day, a simulation of 4 days was performed, corresponding to (D−1) to (D+2) forecast leads. The goal of this study was to determine whether the model forecast spread is lower or greater than the model biases compared to observations. It is shown that the differences between observations and the model are always higher than those between the forecasts. Among all forcing types used in the chemistry-transport model, it is shown that the strong bias and other related low forecast scores are mainly due to the forecast accuracy of the wind speed, which is used both for the mineral dust emissions (a threshold process) and for the long-range transport of aerosol: the surface wind speed forecast spread can reach 50 %, leading to mineral dust emission forecast spread of up to 30 %. These variations are responsible for a moderate forecast spread of the surface PM10 (a few percentage points) and for a large spread (more than 50 %) in the mineral dust concentration at higher altitudes, leading to a mean AOD (aerosol optical depth) forecast spread of ±10 %.

1 Introduction The regional air quality originally focused on photochemical pollution such as ozone and nitrogen dioxides, (Fenger, 2009).This interest was partly motivated by the European "Air Quality Directives" of 1996 that specified policies to reduce air pollution, by gaseous species only at that time (Monks et al., 2009).More recently, the need for a better understanding of aerosols was taken into account in this reg-Correspondence to: Laurent Menut, menut@lmd.polytechnique.frulation framework.While the particulate matter with a diameter less than 10µm (called PM 10 ) has been controlled 10 for many years, the last ten years showed intensification of aerosol monitoring, in particular through the added routine measurements of PM 2.5 (European Union, 2008).In this context, the Mediterranean is well known as a hot spot for its high aerosol concentrations and high spatial and temporal 15 variability (Millan et al., 2005).
Aerosol sources and sink studies remain difficult, since the particulate matter includes lots of different components: several chemical species or materials (organic matter, sulfates, nitrates, ammonia, mineral dust, sea salt etc.), several sizes 20 and shapes, several origins in space, lifetimes, potential direct and indirect effects on radiation, cloud formation, etc.
In order to reduce a potential damage due to the too high aerosol concentrations, it is thus necessary to improve our knowledge of all these aspects (Carslaw et al., 2010).

25
A way to reduce atmospheric pollution is to accurately forecast atmospheric concentrations in order to be able to act at the right time and place to reduce the anthropogenic part of the emissions.This remains a challenge today, and forecast systems often miss large pollution events.Currently, 30 the main effort in Europe to forecast air quality is conducted with the MACC-II system (Marécal et al., 2015), a continuation of the first European multi-model forecast platform (Hollingsworth et al., 2008).This platform itself builds on the first air quality forecast system in Europe, based on the 35 CHIMERE model (Rouïl et al., 2009).
Some previous studies tried to identify and reduce the forecast error.The study of Pérez et al. (2006) is one of the first to explore the interest to couple mineral dust concentrations and radiation to improve aerosol forecast.Manders et al. (2009) 40 quantified the capability of the LOTOS-EUROS system to forecast PM 10 .By reducing some systematic identified biases, Borrego et al. (2011) showed the forecast could be improved over Portugal.Another way to improve forecast is to reduce biases by increasing realism in the aerosol representa-tion, as presented by Mulcahy et al. (2014) for the Met Office global numerical weather prediction model.More recently, several studies showed that data assimilation can reduce the forecast error by constraining the forecast initial conditions, as in (Niu et al., 2008) and (Curier et al., 2012).In all these studies, the bias and variability were considered together.Other frameworks provide daily experimental forecast such as DREAM (Pérez et al., 2007) and SKYRON (Spyrou et al., 2013), mainly focusing on mineral dust.
The goal of this study is to estimate the relative contributions of two modelling aspects, the bias and the variability, by comparing several forecasts to observations.The main question is to what extent the differences between observed and modelled concentrations are caused by modelling errors and by the nonlinear variability of the atmospheric system?
To answer this question, we use the same measurements and model configurations as the ones presented in (Menut et al., 2015).The added value of the present study is the use of this modelling platform in a forecast mode.Section 2 presents the ADRIMED project and the observation sites used.Section 3 presents the modelling system and the forecast setup.Sections 4 and 5 present the forecasted meteorological fields and emissions, respectively.Section 6 presents aerosol optical depth and concentration forecast results.Conclusions and perspectives are presented in Section 7.

The ADRIMED project and the observations used
In this study, we take advantage from the CHARMEX program (Dulac et al., 2013), and more precisely from the ADRIMED project studying the atmospheric composition during June and July 2013 over the Mediterranean area (Mallet, 2014).During this period the ADRIMED project experimental part was measuring the atmospheric composition in the Western-Mediterranean region.At the same time regional models were running real-time forecasts to help optimize the time and location of these measurements.The modeling goal in this case was not to analyze meteorology and chemical composition during a long time period, but to quickly provide an insight in the current state of the atmosphere and its probable evolution over the next few days.Table 1 summarizes the measurement site locations used in this study, giving the longitude, latitude, and altitude above the sea level (ASL) for each site.
To compare the meteorological variables, the European Climate Gridded dataset (E-OBS) daily averaged data were used (Haylock et al., 2008).The dataset contains data for 2m temperature and precipitation rate collected from several thousands of meteorological stations throughout Europe and the Mediterranean area.The data are processed through a series of quality tests to remove errors and unrealistic values.The Aerosol Optical Depth (AOD) was compared using the hourly measurements of the AERONET (AErosol RObotic NETwork) photometers (Dubovik and King, 2000).
The AOD data are recorded by numerous stations deployed around the world.Several quality levels are available in the AERONET database (http://aeronet.gsfc.nasa.gov/); in this 100 study the level 2.0 is used.The PM 10 surface concentrations are measured by the EEA (European Environmental Agency (Guerreiro et al., 2013)

The modeling system
The modeling system is composed of several models: the WRF regional meteorological model, the CHIMERE chemistry-transport model, and additional individual models 110 for emission flux estimations.All these models are integrated in a modelling platform usable both in analysis and forecast mode.This section first describes WRF and CHIMERE models and then the forecast modeling platform.Note that the model configuration (domains, simulated period, model 115 setup) used here in a forecast context is the same as the one used in an analysis context presented in (Menut et al., 2015).

The meteorological model WRF
The meteorological variables are simulated with the WRF regional model, version 3.5.1.The model is used in its non-120 hydrostatic configuration, with a constant horizontal resolu-tion of 60 km × 60 km and 28 vertical levels from the surface to 50 hPa.The Single Moment-5 class microphysics scheme is used allowing for mixed phase processes and super cooled water (Hong et al., 2004).The radiation scheme is the 125 RRTMG scheme with the MCICA method of random cloud overlap (Mlawer et al., 1997).The surface layer scheme is based on the Monin-Obukhov theory with the Carslon-Boland viscous sub-layer.The land surface physics is calculated using the Noah Land Surface Model scheme with 130 four soil temperature and moisture layers (Chen and Dudhia, 2001).The planetary boundary layer physics is treated using the Yonsei University scheme (Hong et al., 2006), and the cumulus parameterization is based on the ensemble scheme of Grell and Devenyi (2002). 135 The global boundary condition fields used are those of the National Centers for Environmental Prediction (NCEP), Global Forecast System (GFS), (Sun et al., 2010).In order to preserve both large-scale circulation features and smallscale gradients and variability, the 'spectral nudging' was 140 used.This nudging method has been already evaluated in regional models Von Storch et al. (2000).In this study, the spectral nudging was applied to all wavelengths greater than ≈2000 km (wavenumbers less than 3 in latitude and longitude) for wind, temperature, and humidity above 850 hPa.

145
This configuration allows the regional model to create its own structures within the boundary layer and yet to follow the large-scale meteorological fields.

The chemistry-transport model CHIMERE
CHIMERE is a chemistry-transport model able to simulate 150 concentration fields of gaseous and aerosols species at a regional scale.The model is off-line, which means that it requires pre-calculated meteorological fields.In this study, we used the version fully described by Menut et al. (2013a).The horizontal domain is the same as the one of WRF.For the vertical grid, the 28 vertical levels are projected onto the 20 levels of the CHIMERE mesh.
The gaseous species are calculated using the MEL-CHIOR 2 scheme, and the aerosols are parameterized according to Bessagnet et al. (2004).This module takes into account species such as sulfate, nitrate, ammonia, primary organic (OC) and black carbon (BC), secondary organic aerosols (SOA), sea salt, dust, and water.These aerosols are represented using nine bins with diameters ranging from 40 nm to 20 µm.The life cycle of these aerosols is completely represented, with nucleation of sulfuric acid, coagulation, adsorption/desorption, wet and dry deposition and scavenging.The scavenging is represented both by coagulation with cloud droplets and precipitation.The formation of SOA is also taken into account.
The anthropogenic emissions are estimated using the same methodology as the one described by Menut et al. (2012) but with the HTAP masses as input data.These masses were prepared by the EDGAR Team, using inventories based on MICS-Asia, EPA-US/Canada, and TNO databases (http://edgar.jrc.ec.europa.eu/htapv2).Biogenic emissions are calculated using the MEGAN scheme (Guenther et al., 2006) that provides fluxes of isoprene, terpene, and pinenes.In addition to this version, several processes were improved and added in the framework of this study.First, the mineral 180 dust emissions are now calculated using new soil and surface databases, as described by Menut et al. (2013b).Second, emission fluxes produced by vegetation fires are estimated using the new high resolution fire model presented by Turquety et al. (2014).And finally, the photolysis rates are 185 explicitly calculated using the FastJ radiation module (Wild et al., 2000) fully described by Mailler et al. (2015).

The forecast configuration
Even though the WRF and CHIMERE models are regularly updated, the forecast configuration of these models remains 190 the same and was previously used in many studies, as listed in (Menut and Bessagnet, 2010).More precisely, this forecast configuration was used during the ESCOMPTE project in the south of France (Menut et al., 2005) and during the AMMA experimental campaign for mineral dust aerosols 195 in western Africa (Menut et al., 2009).CHIMERE is also used in an operational context since 2003 for the PREVAIR French air quality forecast (Honoré et al., 2008;Rouïl et al., 2009) and in the MACC European project (Inness et al., 2013).

200
This forecast system is presented in Figure 1 .The first step is to calculate forecasted regional meteorology.The global GFS/NCEP forecast fields are used to force the regional WRF3.5.1 model from (D-1) (i.e the day before) to (D+2) (two days in advance).The WRF results are then used 205 for several calculations: (i) the surface emission fluxes, (ii) the transport and mixing of gaseous and aerosol species with CHIMERE.For the specific case of the vegetation fire emissions, satellite observations of fire activity (MODIS near-real time detection) during the previous day are analyzed to de-210 rive the corresponding burned area.These are then used as input to the high resolution fire emissions model (Turquety et al., 2014), assuming fires will continue burning during the first 72h of the forecast period.The biogenic and mineral dust emission fluxes depend on the meteorology, while 215 the anthropogenic emissions are only dependent on the week day.The initial conditions for gas and aerosol concentrations are taken from the forecast of the day before.In practice, this means that the system was launched several days before the first day for the first forecast of the period in order to have a 220 correct spin up.
In this study, the simulation was performed from 10 June to 5 July 2013.Each day, a simulation of four days is performed, from (D-1) to (D+2).For each modelled period, meteorological parameters, gas and aerosols species are cal-225 culated hourly on the domain grid.Thus, for each of these parameters, each grid cell and each hour of the period, this Fig. 1.The forecast modelling system.This system includes the download of global meteorological fields, the simulations of the regional models WRF and CHIMERE, and the calculation of numerous emissions fluxes for gas and aerosols species and corresponding to anthropogenic, biogenic, vegetation fires, sea salt and mineral dust emissions.Each day, four days are modelled and the current day (D+0) is used as initialization for the next day forecast (D-1).
allows to have four different values.By comparing these four values, we can quantify the forecast variability.Our analysis focuses on the period from the 14th to the 26th June 2013, 230 identified as the period with the most interesting pollution events during the ADRIMED project.

Calculation of the statistical scores
To compare the forecast results with observations, the following statistical scores are used.The variables O t and M t 235 stand for the observed and modelled values, respectively, at time t.The mean value X N is with N the total number of data used for the calculation.To quantify the temporal variability of the model compared to 240 the observations, the Pearson's product moment correlation coefficient R is calculated as: The Pearson correlation coefficient is the ratio of the covariance between two data sets O t and M t and the product of 245 their two standard deviations.A value of 1 is a complete positive correlation.Similarly, a value of -1 represents a complete negative correlation.
To quantify the mean differences between observations and model results, the bias and the root mean square error 250 (RMSE) are estimated as: For the precipitation amount, it is more convenient to use statistical scores based on the Hit Rate.In terms of 255 its relevance in chemistry-transport modelling, the key factors are space and time variability.In the presence of precipitation, the whole aerosol column is scavenged, and even if the precipitation rate is under-or overestimated, aerosols are quickly deposited.The Hit Rate score is de-260 fined as the following: for a threshold arbitrarily chosen as P r T =0.1mm/day (i.e., there is precipitation this day at this site), the event is considered as true if P r > P r T .Every time this condition is true for both observations and the model, an increase of an "a" value occurs.Every time the condition is 265 true for the observations and false for the model, an increase of a "c" value occurs.The Hit Rate HR is then defined as: The target value for the Hit Rate is 1, meaning that the model was able to capture all the observed events.

Predictability of meteorological parameters
Due to many processes, atmospheric concentrations of trace gases and aerosols are very sensitive to the meteorological fields.First, some of the sources are directly dependent on the near-surface meteorology: (i) mineral dust emissions depend on the surface wind speed, (ii) biogenic emissions depend on temperature and radiation, and (iii) fire emissions depend on the soil moisture (for fire efficiency) and the boundary layer dynamics (for the pyroconvection).Second, during the transport atmospheric species will be influenced by: (i) wind, pressure, humidity, and temperature for the boundary layer dynamics and tropospheric long-range transport, and (ii) clouds and radiation attenuation for the photochemistry.Finally, the sinks of atmospheric species are mainly (i) surface layer turbulence, acting on gas and aerosols dry deposition, and (ii) precipitation via aerosol scavenging.In order to understand the different impacts of meteorological variability on the aerosol concentrations, we focus on temperature, wind speed, and precipitation.The forecast bias and spread for the 2m temperature (T 2m ) are examined at the locations where E-OBS data are available, and the results are presented in Table 2 .In general, the correlations between measurements and modelled values are hith, with values between 0.74 (Zorita, D+2) and 0.98 295 (Aranjuez, Cordoba, Gap).Only one location, Bastia, shows a positive bias (with values from 0.7 to 0.79).All other locations show negative biases ranging between -1.58 (Cordoba, D+2) and 4.02 (Baceno, D-1).This shows in general that the model underestimates the mean daily 2m temperature over 300 the whole simulation domain.This result is consistent with the previous study of Wyszogrodzki et al. (2013) reporting a negative bias in WRF simulations over the United States.
The biases and correlations are found to fluctuate depending on the forecast range.However, these fluctuations are 305 fairly low, with no significant trends in terms of the impact on atmospheric pollutant concentrations.In addition, it appears that the differences between observations and model are always higher than those between several forecast leads; this implies that the model is generally biased, and that the 310 chaotic character of the forecast is low compared to this bias.This was recently discussed by Zhang et al. (2013) who showed that the meteorological forecast accuracy with the WRF model strongly depends on the predictability of the lower-atmospheric boundary layer, especially when synoptic 315 forcing is weak.These conditions are the most common in case of air pollution peaks.
Figure 2 shows the time series of the modelled 2m temperature differences between the forecasts.Results are presented for three sites, Banizoumbou, Bastia, and Lampedusa, 320 which are of interest in terms of other variables, such as the mineral dust, sea salt, and biogenic emissions.Note that these percentages are calculated using temperature values in Kelvin.The maximum differences are calculated for Banizoumbou: over the whole period, values range from ≈ -2 to 325 +2 % (for a mean value of 300K, a variability of ± 6 K).In Bastia, the maximum differences are lower: ≈ -0.5 to +0.5 %.Finally, in Lampedusa the differences may be considered as negligible with values less than 0.2% (less than 0.6 K).

Wind speed and direction 330
The wind speed is a key variable in meteorology and chemistry-transport modelling.Close to the surface (represented by the 10m wind speed), it drives mineral dust and sea salt emissions, the diurnal cycle of the boundary layer convection, and the dry deposition.At higher altitudes, it 335 determines the horizontal transport.In order to quantify the wind speed spread between the forecasts, times series and vertical profiles are presented and discussed in the next section.

Time series of differences 340
The forecast spread of |U | 10m is quantified at the same sites as those used for the 2m temperature.The results are shown in Figure 3 as the percentage of differences between the (D-1) forecast and the other forecasts (D+0, D+1, and Compared to the 2m temperature, the 10m wind speed variability between forecasts is higher.There is no systematic bias; the differences range from 0 to 250%, at wind speeds between 0.1 and 10 ms −1 .For the site of Banizoumbou, mineral dust emissions are sensitive to the wind speed.It is known that saltation occurs for wind speed values up to ≈ 7 ms −1 (even though this absolute value can depend on the soil texture and the landuse).A variability of ± 1 ms −1 (low in absolute value) can have a large impact on mineral dust emission fluxes.For the sites of Bastia and Lampedusa, the forecast differences are lower: being situated on Islands, these sites have a more stable 10m wind speed than over hot and dry land (Banizoumbou).Even though the values are lower, they remain high in terms of differences: up to 150% in Bastia and 130% in Lampedusa.

Vertical profiles
Figure 4 presents vertical profiles of mean wind speed and direction for two locations, Bastia and Lampedusa.The pro-365 files are shown for the whole atmospheric column modelled by CHIMERE, from the surface to 8000 m AGL for the 21 June 2013 at 12:00 UTC.The first result is that the spread between forecasts is higher for the wind speed than for the wind direction.This spread is observed at all altitudes and 370 thus would have an impact both for surface emissions and long-range transport.For example, the wind speed at Lampedusa at around 3000 m AGL ranges from 2 (D-1) to 10 m s −1 (D+2).The forecasted aerosol plumes could be advected too fast in the middle of the Mediterranean in the (D+2) forecast.

Precipitation rates
Table 3 presents the Hit Rates and biases between the E-OBS observations and modelled values.For the sites where precipitation amount was observed and / or modelled, the results show that this variable is correctly modelled in terms of time frequency, but less well in terms of magnitude.In general, when a precipitation event is observed, it is often reproduced by the model.The precipitation intensity appears to be more difficult to simulate, often with a factor of 2 differences (under-or overestimated).For a chemistry-transport model, independent of the meteorology, the time occurrence is more important than the magnitude: the scavenging schemes lead to the total cleaning of the atmospheric column when a precipitation event is diagnosed.This meteorological parameter remains difficult to model but is known to have a large impact on forecast accuracy Eder et al. (2006).The Hit Rate and the spread between forecasts show this parameter explains some discrepancies between forecasts and observations.

Predictability of emissions
The predictability of emissions is quantified for the mineral 395 dust and biogenic emissions.Anthropogenic emissions are not hourly or daily meteorology-dependent, and their time variability is therefore not considered here.For the fire emissions, the model is not able to forecast the burned areas in advance.Each day, the burned areas of the day before are used 400 for the whole period to forecast: the main varying parameter is kept constant.In addition, no significant fire events occurred in June 2013.The fires emission variability is thus not considered neither.

Mineral dust emissions 405
Mineral dust emissions depend on the soil texture, the surface with the landuse and the surface layer wind speed.At the regional scale and over a few days, there is no variability of the soil and surfaces characteristics.On the other hand, the surface layer wind speed can vary a lot.Mineral dust emissions 410 are strongly dependent on the wind speed and thus the corresponding friction velocity u * (Menut et al., 2013b).These dynamical variables act in a non-linear way: the mineral dust emission occurs only if the friction velocity is greater than a threshold value u T * , itself depending on the surface charac-415 teristics.This means that for a small change , in the friction velocity (parameterized using the 10m wind speed), the mineral dust emission could be either zero (if u * = u T * -) or nonzero (if u * = u T * + ). Figure 5 presents two maps for the mineral dust fluxes.

420
The map for the 20 June 2013 is shown as an example, after daily cumulating the hourly fluxes calculated by the model.For this day, the emissions mainly occur over western Africa and Saudi Arabia.Depending on the location, these fluxes range from 0.1 to more than 20 g m −2 day −1 .The second map shows the difference between the fluxes calculated for the (D-1) and (D+2) forecasts.For the region of the highest fluxes, the absolute differences are quite large, i.e., of the same order of magnitude as the flux itself.In order to quantify the forecast spread in a synthetic way, 430 the mineral dust emission fluxes are daily cumulated over the whole simulation domain.The values are presented in Figure 6 (top) and expressed in Tg day −1 .These results show that the fluxes are close between the forecasts: the two main peaks are modelled for the 25 and 28 June with the same order of magnitude.Since the fluxes depend mainly on the wind speed, the latter means that the model is stable at the synoptic scale, and the mean large-scale wind patterns are reproduced regardless of the forecast lead.The same results are seen in terms of the relative differences in Figure 6 (bottom).The sign of these differences varies in time showing large day-to-day variability.The maximum values of differences are ± 30% of the maximum daily flux.Logically, the longer the forecast lead (i.e., for (D+1) and (D+2)), the higher the differences.The largest differences do not occur 445 for the highest absolute values: dust emission being a threshold process, when a high wind speed is forecasted, this is generally true for all forecasts, and the emission fluxes are simulated in a similar way.But, when the wind speed is close to the threshold, a large spread between the forecast leads can 450 occur, as for the 22 June for example.The fluxes correspond to the minimum over the whole period, but the differences are the largest.

Biogenic emissions
The biogenic emissions are sensitive to the temperature and 455 the photosynthetically active radiation (PAR).Over vegetative areas, some changes in these meteorological values could impact the isoprene and terpene emission fluxes.As for the dust emissions, the biogenic emissions are cumulated over the whole simulation domain.The time series are pre-460 sented in Figure 7 (top).A moderate day to day variability is modelled over the whole period: starting with a low value of 2.2 10 9 molecules day −1 , a maximum of 2.4 10 9 molecules day −1 is reached on 17 June, followed by a monotonic decrease to 1.8 10 9 molecules day −1 .The relative differences 465 (%) are shown in Figure 7 (bottom).For all forecasts, the same tendency is observed: the longer the forecast lead, the larger the spread in the flux differences.The differences are moderate, between -2 and +6%.This is consistent with the low differences between the forecasts of the 2m temperature, 470 the latter being the main driver for biogenic emissions.

Predictability of aerosol
Using the meteorological variables and emission fluxes analyzed above, the hourly concentrations of gaseous and aerosol species are simulated with the CHIMERE model.
Here we focus on the aerosol forecast.First, surface concentrations of PM 10 are compared to observations using statistical scores.The sea-salt and mineral dust vertical profiles are discussed.Finally, the Aerosol Optical Depth (AOD) forecast analysis is presented.

PM 10 statistical scores
Results for the PM 10 are presented in Table 4 .The correlations are lower and show higher variability compared to the meteorological variables.The low values are mainly due to the short study period: the contribution of the long-range transport to the aerosol variability is not represented.The mean bias varies from 0 to 14 µg m −3 , being within the range of regional chemistry-transport models (REF).These scores are very different from site to site, but for each site they remain close for different forecast leads.This indicates that the main errors in the forecast system are still caused by the aerosol representation and its forcings, rather than their chaotic character during a forecast.An improvement in the aerosol representation was also shown to improve the forecast score in the Met-Office weather prediction model (Mulcahy et al., 2014).Table 4. Scores for the comparison between observed and modelled hourly PM10 surface concentrations (µg m −3 ).For each site, n obs is the number of hourly valid PM10 observed surface concentrations.P M10 is the mean PM10 value averaged over the whole period for observations and model.The statistical scores presented are the correlation R, the RMSE, and the bias.For each site, the highest correlation is shown in bold.

Vertical profiles of sea-salt and mineral dust
The PM 10 represent the agregation of numerous aerosol species.In order to better understand the forecast variability, vertical profiles are presented for the two dominant species  highest concentrations peak at an altitude of 3000m AGL.For this date, the concentration maximum is ≈ 17 µg m −3 at Bastia and ≈ 350 µg m −3 at Lampedusa.In the two cases, the peaks correspond to the long-range transport of African dust emissions.While the forecast spread is moderate at Bastia, it is large at Lampedusa, with values ranging from 200 to 350 µg m −3 .These differences result from the spread previously discussed for the wind: directly involved in both emission and transport, the wind speed forecast spread impacts the concentration spread at some altitude.A correct representation of the altitude of these dense layers is a crucial point, as previously shown by Wang et al. (2014), using lidar data assimilation to improve aerosol forecast.
The sea salt vertical profiles show that the highest concentrations are close to the surface.This makes sense, since these two sites are on islands in the Mediterranean sea and thus close to the emission sources.Compared to mineral dust, the absolute values of the concentrations are low.But depending on the forecast lead from (D-1) to (D+2), the vari-ability can be high and of the same order of magnitude as the concentrations.In this case, the forecast spread can be directly related to the 10m wind speed used in the model for the emission flux calculations.

Aerosol Optical Depth 530
The Aerosol Optical Depth in another way to represent the aerosol concentration evolution over a large domain.By vertically integrating the aerosol concentrations that are optically active in a specific wavelength (500 nm in this study), the AOD can be an indicator of the daily evolution of 535 aerosols related to the long-range transport.In addition, the dense network of AERONET enables to quantify the realism of the aerosol transport modelling for numerous locations.As for the previous parameters, the correlation values vary a lot between the studied locations.This represents the model's ability to reproduce dense plumes at the right time and place in the domain.But for one location, the values remain close between the forecasts.For example at Forth Crete, the mean AOD is between 0.115 and 0.120, with an observed mean value of 0.099.For all locations, the bias is mainly positive except the Banizoumbou site.The two highest bias values are at Banzoumbou and Dakar in Africa, close to the mineral dust sources.These mineral dust con-550 centrations are the main contributors to the AOD over this region during this summertime period.For the other sites, the bias remains lower than 0.1, and the correlations are ≈ 0.9 at Izana, 0.6 at Forth Crete and 0.8 at Lampedusa, showing the model's ability to capture aerosol plumes far from the 555 main African sources.
The spread between the forecasts is also seen on the AOD maps, Figure 9 .The daily averaged AOD is shown for the 20 June as an example.This day was identified as the one with a dense plume of mineral dust spreading from Africa 560 to the South of Europe.The highest AOD peaks are located in Western Africa and Saudi Aarabia, with maximum values of ≈ 1.8.The plume over Europe shows values between 0.1 and 1.The three other maps represent the absolute difference between the daily averaged map of 20 June (D-1) and the forecasts for the same day: (D+0), (D+1), and (D+2).Logically, the longer the forecast lead, the greater the differences between them.
First, these maps show that the largest differences are located in Africa, where mineral dusts are emitted and where 570 the highest AOD are calculated, such as the hot spots located in Senegal and Yemen.The differences appear as plumes, reflecting the fact that they are caused by both emissions and transport.Another interesting point is that these differences are not spatially homogeneous.The differences represent 575 "dipoles" of the opposite sign, and they increase with the forecast lead.The largest gradients are located where the highest AOD are simulated.These locations correspond to the largest emissions and transport of mineral dust.The latter are very sensitive to the wind speed and direction, and the 580 gradients reflect the impact of the wind direction variability between the forecast leads.The long-range transport can also lead to the differences of the opposite sign: the longer the transport of dense plumes, the more pronounced are the differences.Finally, with AOD values ranging between 0 to 585 2, the absolute differences between all forecasts can reach ± 0.1 (≈ 10%).

Conclusions
This study was dedicated to the quantification of the spread between several aerosol forecasts over the Euro-590 Mediterranena area.This was done in the framework of the ADRIMED campaign (Mallet, 2014), as part of the CHARMEX project (Dulac et al., 2013).The studied period, the domain, and the model set-up are the same as those presented in (Menut et al., 2015).In the present study, the 595 model was run every day for four-day-long simulations.By comparing several forecasts between them and with observations, we quantified the relative impacts of the model biases and the chaotic character of a forecast on the forecast accuracy.

600
In order to quantify the forecast accuracy of aerosols, several forcing parameters are studied.For the meteorological parameters, it was shown that the 2m temperature is mainly biased, but well correlated to the measurements and with a low spread between the forecasts.The precipitation is simu-605 lated moderately well: on average two events of three are reproduced, and precipitation rate is biased.But since its main effect is fast scavenging of the atmospheric column, this parameter is modelled sufficiently well to ensure a low impact of the forecast lead on the aerosol content.On the other hand, 610 the wind was found to have a high variability between the forecasts.The 10m wind speed can have a day-to-day variability of ±150%.The mineral dust and biogenic emissions were also studied, depending on the forecast range.The biogenic emissions show a low variability between the forecasts, 615 due to the forecast stability of the temperature.But the mineral dust emission forecast is highly variable, with values of ±40% between the forecast leads.This is a direct effect of the wind forecast variability, acting on both emission fluxes (a threshold process) and the long-range transport of aerosol.

620
The forecast spread of aerosol concentrations was presented in terms of surface PM 10 , vertical profiles of sea-salt and mineral dust, and Aerosol Optical Depth.The surface PM 10 are compared to Airbase measurements in Europe: the correlation is moderate (from 0 to 0.49), and the bias varies 625 from 0 to -14 µg m −3 .However these scores are weakly variable, the bias and the correlation remaining stable with increasing the forecast lead.The forecasts are more variable in terms of vertical profiles: for sea-salt, a wide spread of 100% is found close to the surface (and thus to their maritime 630 emissions), whereas mineral dust concentrations are strongly variable (± 50%) at certain altitude, in the center of dense plumes.The AOD were compared using the Aeronet measurements.The correlations are higher than for PM 10 , and the bias is weakly variable between the forecasts.

635
Finally, there are two main conclusions for this study: (i) the differences between observations and the model remain higher than between the forecasts.When high differences between the model results and the observations occur, they are mainly due to the model biases rather than forecast lead.

640
(ii) among all studied variables, the highest variability of the forecast is due to the wind speed and direction.The wind is at the origin of mineral dust and sea salt emissions, as well as the long-range transport of these long-lived species; therefore the differences in the forecasted wind speed and direc-645 tion are at the origin of the spread between the aerosol concentration forecasts.

Fig. 2 .
Fig.2.Time series of hourly modelled 2m temperature differences for several sites.Each line corresponds to a difference between the forecast (D+0), (D+1) or (D+2) and the simulation for the day before (D-1).Results are expressed as percentage of differences.

Fig. 3 .
Fig. 3. Time series of relative differences (%) in hourly 10m wind speed (m/s) for sites Banizoumbou, Bastia, and Lampedusa for several forecast leads.The data are shown only for the values where |U | (D−1) > 2ms −1 , to avoid too large and unrealistic values at low wind speed.

Fig. 4 .
Fig. 4. Modelled vertical profiles of the wind speed (m s −1 ) for the cells corresponding to the locations of Bastia and Lampedusa.Profiles are presented for the 21 June 2013 at 12:00 UTC and for the four forecasts, from (D-1) to (D+2).

Fig. 6 .
Fig. 6.Top: time series of daily mineral dust fluxes, spatially cumulated over the modelled domain.Units are in Tg day −1 .Bottom: relative differences between the fluxes (%).

Fig. 7 .
Fig. 7. (top) Time series of the daily isoprene fluxes, spatially cumulated over the modelled domain.Units are in 10 9 molecules day −1 .(bottom) Differences between the flux expressed percentages.

500
in the budget: sea salt and mineral dust.The profiles are extracted from the model outputs for the Bastia and Lampedusa locations for the 21 June 2013, 12:00 UTC.The results are presented in Figure8.The mineral dust vertical profile shows low concentrations 505 close to the surface, with values lower than 5 µg m −3 .The

Fig. 8 .
Fig. 8. Vertical profiles of mineral dust and sea salt concentrations (µg m −3 ).In each figure, the four forecasts are presented from (D-1) to (D+2).Results are presented for Bastia and Lampedusa for the 21 June 2013.

Fig. 9 .
Fig. 9. Map of modelled (D-1) AOD for the 20 June 2013 and maps of AOD differences between the several forecasts.

Table 1 .
) running the AirBase database.It contains hourly surface concentration measurements and information submitted by the participating countries through-Characteristics of the AirBase and AERONET stations used in this study.Note that the AirBase Italian stations of Chitignano, Baceno, Schivenoglia and Vercelli provide daily averaged values, while all other stations provide hourly (but not regular) measurements.The altitude is in meters Above Sea Level (ASL).

Table 2 .
Scores for the modelled hourly 2m temperature compared to the measurements.For each forecast, the correlation R (0 to 1) and the bias (in o C) are presented.

Table 3 .
Scores for the modelled daily cumulated total precipitation (mm/day) compared to the measurements.For each forecast, Hit Rate HR and bias are presented.

Table 5 .
Table5presents the statistical scores comparing the observed and modelled hourly AOD.Scores for the comparison between observed and modelled hourly AOD.For each site, n obs is the number of hourly valid AERONET AOD observations.AOD is the AOD averaged over the whole period for observations and the model.The statistical scores presented are the correlation R, the root mean square error RMSE, and the bias.