Aerosol optical properties over Europe: an evaluation of the AQMEII Phase 3 simulations against satellite observations

The main uncertainties in estimates of changes in the Earth’s energy budget are related to the role of atmospheric aerosols. These changes are caused mainly by aerosol-radiation (ARI) and aerosol-cloud interactions (ACI), which heavily depend on aerosol properties. From the 1980s, many international modelling initiatives have studied atmospheric aerosols and their climate effects. Phase 3 of the Air Quality Model Evaluation International Initiative (AQMEII) focuses on evaluating and intercomparing regional and linked global/regional modelling systems by collaborating with the Task Force on the Hemispheric 5 Transport of Air Pollution Phase 2 (HTAP2) initiative. Within this framework, the main aim of this work was to evaluate the representation of aerosol optical depth (AOD) and the Ångström exponent (AE) by the AQMEII Phase 3 simulations over Europe. The evaluation was made using satellite data from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on board the Terra and Aqua platforms. The results indicated that the skills of AQMEII simulations in the AOD representation produced fewer errors than in the AE. Regardless of the models and emissions used, models were skilful at 10 representing the low and medium AOD values observed (below 0.5). However, high values (close to 1.0) were underestimated for biomass burning episodes, and were overestimated for desert dust contributions, related mainly to emission and boundary conditions. Despite this behaviour, the spatial and temporal variability of this variable was well-represented by all the models. Generally, the AE evaluation showed more serious errors than the AOD evaluation. Moreover, the observed variability of this parameter was strongly underestimated in all the simulations. 15 1 Atmos. Chem. Phys. Discuss., https://doi.org/10.5194/acp-2017-1119 Manuscript under review for journal Atmos. Chem. Phys. Discussion started: 5 March 2018 c © Author(s) 2018. CC BY 4.0 License.

On the other hand, several previous studies had evaluated modelled aerosol optical properties against satellite data from a global view. In Ghan et al. (2001), simulated AOD were within a factor of 2 of AVHRR (Advanced Very High Resolution Radiometer) products and the behaviour of the Ångström Exponent (AE), estimated from POLDER (POLarization and Directionality of the Earth's Reflectances) and SeaWiFS (Sea-Viewing Wide Field-of-View Sensor), was similar to that simulated.
Otherwise, both the simulated AOD in Chin et al. (2002) and Reddy et al. (2005) were reproduced with most of the notable 5 features in TOMS (Total Ozone Mapping Spectrometer), AVHRR and MODIS. Moreover, Ginoux et al. (2006) revealed sensitivity to humidity when evaluated against satellite data. Kinne et al. (2003) compared aerosol modules from seven models with MODIS and TOMS, and found large discrepancies over tropical and Southern Hemisphere oceans due to sea salt treatment. Kinne et al. (2006) also discovered a lower simulated AOD among 20 different modules from the AEROCOM Project (0.11 to 0.14) than the satellite AOD composite of MODIS, MISR (Operational Microwave Integrated Retrieval System), AVHRR, 10 TOMS, and POLDER retrievals (0.15).
More recent studies are Colarco et al. (2010), who assessed simulated AOD versus MODIS and MIRS, and found similar seasonal and regional variability and magnitude over downwinds of the Saharan dust plume, a high bias in sulphate-dominated regions of North America and Europe, and a better agreement over ocean when the sea salt burden was reduced by a factor of 2. Furthermore, Zhang et al. (2012) reported a relative difference in AE of 13.8% with a negative (positive) bias over high 15 latitude regions (oceans), but correlated well for AOD in comparison with MODIS. Finally, Liu et al. (2012) evaluated longterm simulations compared with the satellite composite derived by Kinne et al. (2006) and identified a low bias for AOD, but a good representation of the observed geographical and temporal variations of aerosol optical properties.
Similar studies to ours, which made a seasonal comparison over Europe, are: Jeuken et al. (2001), in which the comparison of AOD calculations with ARST-2 (Along Track Scanning Radiometer 2) on board the European ERS-2 satellite showed 20 an average difference of 0.17-0.19, but a good representation of the observed patterns; Solmon et al. (2006), in which the simulated AOD presented a general underestimation (more pronounced over the Mediterranean Basin), but within the range of AERONET and MIRS over northern Europe, and followed spatial patterns of MODIS and TOMS over both Europe and Africa.
As the above-mentioned aerosol optical properties influence ARI and ACI, a good representation of them is, thus, a key issue to reduce the uncertainty of aerosol effects on the Earth's climate system. Curci et al. (2017) used AQMEII Phase 25 3 simulations to evaluate the black carbon absorption against AERONET but no works have evaluated against satellite the seasonal representation of optical properties over Europe by the regional models involved in AQMEII Phase 3. For this reason, our main study aim was to evaluate the representation of two main aerosol optical properties, AOD and AE, using the models of the AQMEII Phase 3 initiative over Europe. The evaluation was made by using the remote-sensing observations from the MODIS sensor. Section 2 provides a brief description of the observational and models data, and the evaluation methodology. 30 Section 3 presents the evaluation results. Finally, Section 4 summarises the main conclusions reached.
particles, taking into account their hygroscopic growth. The optical properties used in the Mie computations originate from the OPAC dataset (Hess et al., 1998).
The ES1 simulation was run by the Regional Atmospheric Modelling Group at the University of Murcia (UMU, Spain). They used the Weather Research Forecasting model online coupled with Chemistry (Grell et al., 2005, WRF-Chem), version 3.6.1.
Meteorological inputs were driven by ECMWF analysis fields. The aerosol module based on the Modal Aerosol Dynamics 5 Model for Europe (Ackermann et al., 1998, MADE), in which secondary organic aerosols (SOA) were incorporated by the Secondary Organic Aerosol Model (Schell et al., 2001, SORGAM), was used. The employed gas phase chemistry model was the Regional Acid Deposition Model, version 2 (Stockwell et al., 1990, RADM2), with 57 chemical species and 158 gas phase reactions, of which 21 are photolytic. Anthropogenic emissions were MACC emissions. Biogenic VOC emissions were computed by applying the Model of Emissions of Gases and Aerosols from the Nature (MEGAN) emissions model (Guenther,10 2006), version 2.04. The MADE/SORGAM sea salt (Gong, 2003) and dust (Shaw et al., 2008) emissions were used.
The IT1 simulation was made at the Ricerca sul Sistema Energetico (RSE, Italy) using the Weather Research Forecasting (WRF) model coupled with the Comprehensive Air Quality Model with Extensions (CAMx), version 6.10. Meteorological inputs were generated using WRF-Chem version 3.4.1. Anthropogenic emission were MACC and biogenic were computed by MEGAN. WRF-Chem was adopted to predict GOCART (Goddard Chemistry Aerosol Radiation and Transport) dust emissions 15 (Ginoux et al., 2001) along with meteorology. Sea salt emissions were computed using published algorithms (de Leeuw et al., 2000;Gong, 2003). The WRF-CAMx pre-processor (ENVIRON, 2014, version 4.2) was used to create the CAMx ready input files by collapsing the 33 vertical layers used by WRF to 14 layers in CAMx, but by maintaining the layers up to 230 m above ground level identical. Aerosol optical properties were estimated by means of the Aerosol Optical DEpth Module (Landi, 2013, AODEM) post-processing tool that was coupled to CAMx regional model. AODEM calculated the optical properties 20 (e.g. AOD, extinction and scattering coefficients, and particle number concentrations) at different wavelengths and size bins starting from the aerosol mass concentration predicted by CAMx. In this work, the Mie theory was applied by dividing the size range (40 nm to 10 µm) into 10 bins, and calculating the hygroscopic growth of each aerosol species in each bin with the Hanel formula. Moreover, particles were assumed to be internally mixed.
The IT2 simulations were run at the University of L'Aquila (Italy) using WRF-Chem (Grell et al., 2005), version 3.6. The 25 modified MADE/VBS aerosol scheme  was included in this version. This scheme is based on MADE to treat inorganic aerosols, along with the VBS approach (Ahmadov et al., 2012). MADE/VBS allows a better representation of the SOA mass. The Regional Atmospheric Chemistry Mechanism -Earth System Research Laboratory (RACM-ESRL) gas phase chemical mechanism (Kim et al., 2009) was used. Anthropogenic emissions were MACC emissions which had been adapted to the chemical mechanism used following the method of Tuccella et al. (2012). As in the IT1 and ES1 simulations, 30 biogenic emissions were calculated online by the MEGAN model (Guenther, 2006). Finally, the meteorological analyses used to initialise WRF were provided by the ECMWF with a horizontal resolution of 0.5 • every 6 h. IT2_M-ARI was run with ARI, while large-scale clouds were solved by a simple module. IT2_M-ARI+ACI took into account ARI and ACI, while aqueous chemistry was solved in convective clouds. As well as ES1, IT2 simulations used the MADE/SORGAM sea salt and dust emissions.  Barnard et al. (2010) andChapman et al. (2009). The composite aerosol optical properties were determined by the Mie theory, summation over all the size bins and wet particles diameters. An overall refractive index for a given size bin, as determined by an volume averaging of complex indexes of refraction associated with each chemical constituent of the aerosol, was used. The inclusion of ACI and ARI in WRF-Chem is described in Chapman et al. (2009) A multimodel ensemble (henceforth referred to as ENSEMBLE) of the available simulations was also evaluated. The results presented herein did not intend to represent an ensemble of opportunity, but were merely calculated as the mean of all the participating simulations. As part of the AQMEII Phase 3 initiative, the available variables of aerosol optical properties were AOD at 470, 550 and 670 nm. (2) AE between 550 and 860 nm over the ocean estimated by the DT algorithm. Although this "combined" AOD product has not yet been validated, it offers a wider coverage. The preliminary estimated error (EE) for the used AE products was 0.45 in the pixels with an AOD > 0.2 (Levy et al., 2013). The selection of this observational data was due to results found by Palacios-Peña et al. (2017a).

5
As Terra and Aqua are in Sun-synchronous orbits around the Earth, MODIS does not provide data over the entire studied domain for each time step. According to Levy et al. (2013), who have established that there is no significant difference between MODIS/AERONET comparability for Terra and Aqua data, we combined the hourly data from both satellites in order to obtain a whole year of data with a wider coverage for each time step than by using the Terra and Aqua data separately.

10
Simulations (Table 1) and satellite data had a different spatial resolution. Henceforth and beforehand, all the data were preprocessed and bilinearly interpolated to a common working grid with a resolution 0.25 • .
As mentioned above, our objective was to evaluate the representation of the main aerosol optical properties: AOD and AE.
MODIS provides AOD at 550 nm and AE between 550 and 860 nm, but from simulations, the available variables were AODs at different wavelengths. Thus in order to evaluate AE from simulations, this variable had to be estimated through the Ångström 15 empirical expression (Ångström, 1929, eq. 1), where λ is the wavelength and β is Ångström's turbidity coefficient.
By rationing equation 1 at two different wavelengths and taking algorithms, AE can be computed from the spectral AOD values (Eck et al., 1999, eq. 2). Hence it is possible to estimate AE between two known wavelengths, and to also use this AE to estimate AOD at other different wavelengths. However, as established Ignatov et al. (1998), retrievals of AE under AOD 20 conditions lower than 0.1 are highly uncertain. For this reason, we chose the criteria to estimate AE over areas with AOD > 0.1. Moreover and according to the EE for the AE products of MODIS, we set the AE values range between -0.5 and 4.0.
Once all the data had the same resolution, and following Equation 2, the simulated AOD and AE were at the MODIS wavelengths. Then the hourly data were evaluated using classical statistics such as: the mean of the individual model-prediction 25 error or bias (e i ); i.e., the mean bias error (MBE); the mean of the absolute error (MAE); and the coefficient of determination (r), according to Willmott et al. (1985), Weil et al. (1992) and Willmott and Matsuura (2005). Before the statistical analysis, a mask that showed the areas where the satellite observations were higher than the 10% of their maximum, was implemented. This mask was implemented in order to carry out a statistical analysis with a number of observations which show more confident results. Figures S1 and S2, in the appendix, show the number of observations used when the mask was implemented.
The MAE was calculated as in Equation 4 and provides an estimation of the magnitude of the error independently of over-5 or underestimation.
The temporal determination coefficient was estimated as in Equation 5 and was used as a measure of the strength of the linear relationship between two variables, in our case, the satellite and simulations values.
Finally, the Kernel probability density functions (PDF) with a broadband of 0.05 were used to evaluate the skills of the simulations to reproduce the variability (temporal and spatial) of the studied variables (AOD and AE).

Results
The next section evaluates the skills of the different AQMEII Phase 3 simulations with respect to the representation AOD and AE representation. The first section shows the model evaluation for the AOD representation and the second for AE. The 15 numerical result of each case is indicated by the numbers represented in each figure. Finally, the skills of the simulations to reproduce the variability of the studied variables (AOD and AE) are analysed using the PDF of each variable.

Model evaluation of the AOD representation
AOD is defined as the integrated extinction coefficient over a vertical atmospheric column and indicates to what degree aerosols avoid light transmission. AOD is not a direct function of the atmospheric load of particles, but can provide us an approximate 20 idea of both atmospheric load of particles and the interaction of these particles with radiation.  When seasonal figures were analysed (Figure 1), the highest values (around 1) were found over the southern part of the domain for all seasons due to frequent Saharan desert dust outbreaks, which affects the Mediterranean Region. Moreover, these desert dust outbreaks were more frequent and strong for spring and summer because the mean AOD values were higher.

5
In spring, the highest temporal mean AOD values, between 0.5 and 1, reached the southern part of the domain. These mean values can be higher than 1 over specific areas.
In summer, the temporal mean AOD values over the southern part of the domain fell between 0.5 and 0.8, but the higher mean AOD values (above 1.5) were found over a large area in Russia and its surrounding areas. During this period, a heat wave and wildfires occurred over this area, which explained this fact. For these two reasons, this season presented the highest mean 10 AOD values when both space and time were considered.
However for autumn and winter, high mean AOD values were also found over the southern part of the domain, but were lower than 0.5. The lower mean AOD value when considering space and time was found in autumn. It is noteworthy that ice, snow and clouds were avoided for the MODIS sensor, and aerosol properties were not retrieved over the areas where they were present (http://darktarget.gsfc.nasa.gov/), which explains the gap over the northern part of the terrestrial domain in winter and  Figure S3 at SM). Over north-western areas and Russia (due to the wildfires that occurred in summer), these simulations underestimated AOD (negative MBE values around -0.05). One main issue is that no clear differences were found when HTAP and MACC emissions were used.

30
When a different chemistry model was employed, minor differences in the error of simulations were found between those made using the CAMx chemistry model (IT1_MACC) and the WRF-Chem (IT2 simulations). These differences were of a similar order of magnitude to that of the differences between the IT2 simulations by including ARI and ACI. However, the One notable point in the underestimation was presented by all the simulations over the south-eastern part of the domain (represented as a blue spot), centred over Iraq. The ES1_MACC simulation did not show this underestimation because of its 15 high AOD values, but presented lower overestimation values (close to 0) over this area than over its surroundings. This small spot was also notable when MAE (Figure 2) was analysed.
In winter (the first column in Figures 1 and 2), all the simulations presented a weak underestimation over the Atlantic Ocean, except for IT1_MACC, which presented a weak overestimation in the northern part (MBE mean of 0.03). The above-mentioned blue spot was clearly defined over a small south-easterly area and was stronger during this season, even for the ES1_MACC In spring (the second column in Figures 1 and 2), the underestimation of AOD was similar to that in winter, but with steeper values. All the simulations presented an overestimation (with different degrees) over the southern part of the domain, the Balkan Peninsula and southern Russia. This overestimation was larger and stronger for the ES1_MACC simulation (MBE of 0.21) and once again presented higher MAE values (0.29). All the simulations, except for IT1_MACC, presented a weak underestimation 30 over the Atlantic Ocean. The IT simulations gave fewer errors than the rest. As in winter, a small south-easterly area (the blue spot) appeared, but was consistent with the maximum MAE values for the IT simulations.
The underestimation of AOD due to the wildfire emissions over Russia and the surrounding areas was one of the most important issues in summer (the third column in Figures 1 and 2). This underestimation was larger and stronger for the IT2 simulations, and was smaller and weaker for the FI1 simulations. Moreover, the aforementioned small area in the south-eastern  In winter the highest determination values (close to 1.0) were found over the north-eastern part of the African continent. In 20 spring, these high values were found over central and eastern parts of Europe and North Africa. In summer, the highest values were mainly over Russian and the surrounding areas and a part of the Atlantic Ocean in the south-western part of the domain.
Finally in autumn, the coefficient of determination values were lower than for the other seasons, and were mainly over the Mediterranean sea and the Atlantic Ocean.

25
AE is a parameter that indicates the relationship between the size of the particles suspended in the atmosphere and the wavelength of the incident light, and, although there is not a direct correspondence between aerosol size and AE, provides an idea of the size of particles. Low AE values are related to coarse particles, such as desert dust or sea salt, and high values are associated with fine particles, such as anthropogenic source particles or biomass burning. The AE values are between 0 (or even slightly negative in coarse mode aerosol episodes) and 4 (Boucher, 2015).  were uniformly distributed in spring over the southern part of the domain, while in summer (JAS, the third column in the first row in the AE figures) the lowest AE values (lower than 0.5) were found mainly over the southern Atlantic Ocean. Values between 2.0 and 2.5 were estimated over north-eastern coasts and the north of the Black and Caspian sea, which were lower than in summer.

5
In this section, the simulations run with the available data were less than they were for AOD. The errors of the simulations made for the whole year are shown in Figure S4 at SM.
Throughout  mation of the IT1_MACC simulations was weaker, but the underestimation was stronger and over a large area over the North and Baltic Seas. The IT2_M-ARI+ACI simulation also produced an overestimation over most of the domain, but it was weaker than that presented in spring. Notwithstanding, this simulation presented a small area of underestimation over the Baltic sea.
However, ENSEMBLE displayed a general underestimation that lowered from the coast to offshore.
In autumn (the fourth column in Figures 4 and 5), the behaviour of simulations was similar to that shown in winter.

10
FI1_HTAP produced a general, but weaker underestimation than in spring and summer.  Figure 6 shows the results of the determination coefficient, which were significant at 90%. These results were worse throughout the year (column third in Figure S4 at SM) than when the different seasons were evaluated. FI1_HTAP and IT1_MACC showed relatively high values (around 0.5) over the Mediterranean Sea, but over this area, all the other simulations presented values above 0.25. Even though the determination values were higher during the seasonal study, it was very difficult to find a clear coefficient of determination pattern. During each season FI1_HTAP was used to present the highest determination values 20 and ENSEMBLE the lowest ones.

Variability
A good approach to evaluate the spatial and temporal variability of a variable is the Probability Density Function (PDF).
This represents the density of counts for each value of the variable. In order to study how the AQMEII Phase 3 simulations represented the variability of AOD and AE, the PDF of both variables for each studied season are shown in Figure 7. In that  Simulations FI1 and IT2 displayed analogous PDFs with the highest probability for the lower AOD simulated values than those observed. In winter (JFM) and autumn (OND), the IT2 simulations presented higher probabilities for the lower AOD values than the FI1 simulations. However in spring (AMJ), these four simulations gave almost equal PDFs.
The ES1 simulation showed a remarkable representation of AOD in all seasons. The PDFs for this simulation estimated higher probabilities for the high AOD values than the other simulations and the observed values. For this reason, the probability 5 of low AOD values was lower than for the rest. ENSEMBLE displayed in JFM, AMJ and OND a high probability for the lower AOD simulated values than those observed, but with ENSEMBLE, the probability for the higher AOD values was higher than for those observed.
The PDFs of the AOD representation were different in summer (JAS). For this season, both IT2 were the simulations that displayed the nearest behaviour to the observed PDF values. As seen in the SM, these simulations displayed the lowest 10 MAE compared with the observed standard deviation. All the simulations and ENSEMBLE in this season presented a higher probability for the high AOD values than those observed.
The second column of Figure 7 represents the PDFs for the AE values. As for AOD, winter and autumn presented similar PDFs. The observed AE values showed a high probability for the low AE values, around 0.5, and a low probability for the high AE values. For spring and summer, the PDFs for the observed values were tray-shaped, with a high probability for the AE 15 values between 0.5 and 1.5.
The behaviour of the other simulations and ENSEMBLE was similar in all seasons. The FI1_HTAP simulation displayed a high probability for very low AE values, between 0 and 0.6. IT1_MACC gave similar PDFs for all the seasons, and a high probability was found for the AE values between 1 and 1.5. FI1_HTAP and IT1_MACC were the simulations with the narrowest PDFs, which indicates that these simulations produced a strong underestimation of the observed variability of AE 20 (as also indicated in the evaluation of the temporal standard deviation shown in the SM).
The PDFs for ES1_MACC, IT2_M-ARI+ACI and ENSEMBLE were wider than for the other two simulations. ES1_MACC showed a high probability for the AE values, between 0 and 1.5. IT2_M-ARI+ACI for the AE values was between 0.5 and 1.5.
ENSEMBLE showed a high probability that ranged from 0 to 1.5. Notwithstanding, these PDFs were narrower than the PDF for the observed values, thus all the simulations underestimated the representation of the AE values. This is observed in the 25 SM, where the estimation of the MBE of the standard deviation gave negative results for all the seasons and simulations.

Summary and Conclusions
Although AQMEII Phase 3 focuses on evaluating and intercomparing regional and linked global/regional modelling systems, an evaluation of the simulations of the front observations was needed. Solazzo et al. (2017) analysed the performance of models for different meteorological variables and chemical species. In order to perform a more detailed analysis of the models' 30 performance, this work focused on evaluating the aerosol optical properties representation by means of AQMEII Phase 3 simulations using satellite sensors. The evaluation of these variables is important because they strongly influence ARI and ACI and, thus, influence the atmospheric aerosol effect on the climate system. injected closer to the ground, making the average smoke transport distance even smaller than for the fixed emission height.
Also Soares et al. (2015) point out, referring to Wooster et al. (2005), that MODIS is not sensitive enough to register the fire radiative power of small or smoldering fires, and thus large fraction of those is missed in the emission data, including also strongly emitting peat fires. The 2010 Russian fires included some huge fires, but also numerous small ones over large areas, and a large fraction of those was probably missed by MODIS. 15 Moreover a high underestimation was produced for all the simulations, irrespectively of the used meteorological and chemical model, over a small area in the south-eastern part of the domain (the above-called "blue spot"). This can be explained by the fact that the emissions inventories used herein only covered European areas (see the Emission Map in the SM), thus the emissions over that area were not considered.
The AOD over the southern part of the domain was overestimated and related mainly to the high dust concentrations ac-20 cording to the boundary conditions. In line with this, Solazzo et al. (2017) found that the error in primary species as dust was strongly affected by the emission and boundary conditions in the AQMEII Phase 3 simulations.
On the whole, the FI1 simulations, which used the ECMWF files as SILAM model input, gave high AOD values because SILAM is known to have slower dry particle deposition than other models. This could explain that, although the band quiet crudely represents size distribution, AOD is also very sensitive. IT1, which used the WRF meteorological conditions as CAMx The AE satellite values were obtained only over sea. High AE values, which indicate fine particles, were found near central European coasts which were higher in summer, probably due to the wildfires that occurred then. Low AE values, which indicate coarse particles, were observed over the southern part of the domain, close to the Saharan desert and over the Atlantic Ocean.

5
It was also noteworthy that the AE values over the Atlantic Ocean were generally much higher in spring and summer than in autumn and winter. This means that the aerosol particles over ocean areas and near the coast in warm months were apparently finer than in colder months. This might be related to two different hypotheses: weaker winds in warm months or hygroscopic growth, which could be greater in cold months generally because of higher relative humidity (RH).
AE modelling skills were lower than for AOD (larger errors  low AE values and thus compensating the tendency for producing high PM 2.5 /PM 10 ratios.
It was not possible to find any clear spatial pattern for the coefficient of determination of the AE representation. One striking fact in this case was that using the mean of all the simulations as ENSEMBLE did not improve this statistical figure. In fact the worse determination results were found for ENSEMBLE. The FI1_HTAP simulation showed the highest determination values.

25
Thus despite the high underestimation of the AE values, it displayed a good skill in the temporal AE representation.
As mentioned above, a good approach to evaluate the spatial and temporal variability of a variable is PDF. A wide PDF indicates wide variability for the studied variable, and a narrow and low variability. For the AOD representation, all the simulations presented similar PDF to the observed values. The behaviour of all the simulations was similar in winter, spring and autumn; simulations FI1 and IT2 presented higher probabilities for lower AOD values than those observed; ES1 presented 30 higher probabilities for high AOD values than those observed due to the above-explained lack of dust gravitational settling.
Finally, IT1 presented the most skilful PDF. Given the probability of obtaining AOD values around 0.5 being higher, the IT2 simulations presented the best skills in summer. One general conclusion was reached from the PDF of the AE values. For this variable, all the simulations in all the studied seasons underestimated temporal and spatial variability. In conclusion, the skills of all the simulations in the AOD representation produced lower errors than in the AE representation.
For AOD, low and medium values were well-represented, but high values presented larger errors. High values due to dust were overestimated because of an overestimation in the boundary conditions. The high AOD values due to biomass burning were underestimated, which should be ascribed to an understated injection height of the total biomass burning emissions or directly to underestimated emissions. Other high AOD values were underestimated because the emissions which produced these high 5 values were not considered. The errors in the AOD representation evidenced the strong influence of emissions and boundary conditions in the estimation of aerosol optical properties. Generally speaking, the models' skills to represent the variability of AOD were acceptable. For AE, the SILAM simulation underestimated the observed values and the WRF coupled CAMx simulation and the WRF-Chem simulations were those with the best skills in the representation of this variable. But for all the simulations, the variability of this variable was underestimated.

10
Following these results, further studies are needed to improve the representation of aerosol optical properties, along with other properties such as atmospheric distribution, hygroscopicity, or the ability to act as cloud condensation nuclei (CCN) and ice nuclei (IN). The matter noted in the representation of aerosol properties can help to gain a better representation of ARI and ACI and aerosol effects on meteorology and climate, and could reduce the grave uncertainty in the estimations of changes in the Earth's radiation budget due to aerosols and clouds.

Data availability
The outputs from the simulations can be obtained by emailing to rbianconi@enviroware.com. MODIS data are publicly available on the MODIS Atmosphere website (https://modis-atmos.gsfc.nasa.gov/MOD04_L2/acquiring.html).
Acknowledgements. This work was conducted under the support of the AQMEII/HTAP Phase III initiative. The authors acknowledge Project REPAIR-CGL2014-59677-R of the Spanish Ministry of Economy and Competitiveness and the FEDER European programme for support 20 to conduct this research. This work has been also possible thanks to fellowship 19677/EE/14, funded by the "Fundación Séneca-Agencia de Ciencia y Tecnología de la Región de Murcia", and the Programme "Jiménez de la