Added Value of Aerosol-Cloud Interactions for Representing Aerosol Optical Depth in an Online Coupled Climate-Chemistry Model over Europe

: Aerosol-cloud interactions (ACI) represent one of the most important sources of uncertainties in climate modelling. In this sense, realistic simulations of ACI are needed for a better understanding of the complex interactions between air pollution and the climate system. This work quantiﬁes the added value of including ACI in an online coupled climate/chemistry model (WRF-Chem, 0.44 ◦ horizontal resolution, years 2003 to 2010) in order to assess whether there is an improvement in the representation of aerosol optical depth (AOD). Modelling results for each species have been evaluated against the Copernicus Atmosphere Monitoring Service (CAMS) reanalysis, and AOD at 675 nm has been compared to AERONET data. Results indicate that the improvements of the monthly biases are around 8% for total AOD550 when including ACI, reaching 20% for the monthly bias in AOD550 coming from dust. Moreover, the temporal representation of AOD550 largely improves (increase in the Pearson time correlation coefﬁcients), ranging from 6% to 20% depending on the chemical species considered. The beneﬁts from this improvement overcome the problems derived from the high computational time required in ACI simulations (eight times higher with respect to simulations not including aerosol-cloud interactions). Author Contributions: Conceptualization, P.J.-G., J.P.M. and L.P.-P.; methodology, L.P.-P., J.P.M., J.M.L.-R., S.J., J.J.G.-N., R.L.-P., J.R. and P.J.-G.; software, J.M.L.-R. and S.J.; validation, P.J.-G. and L.P.-P.; formal analysis, L.P.-P., J.P.M., J.M.L.-R., S.J., J.J.G.-N., R.L.-P., J.R. and P.J.-G.; investigation, L.P.-P., J.P.M., J.M.L.-R., S.J., J.J.G.-N., R.L.-P., J.R. and P.J.-G.; writing–original draft preparation, L.P.-P. and P.J.-G.; writing–review and editing, L.P.-P., J.P.M., J.M.L.-R., S.J., J.J.G.-N., R.L.-P., J.R. and P.J.-G..; visualization, L.P.-P., J.J.G.-N.; supervision, P.J.-G.; project administration, P.J.-G.; funding acquisition, P.J.-G., J.P.M. All authors have read and agreed to the published version of the


Introduction
The presence of pollutants in the atmosphere, especially aerosols, impacts on the radiative balance, and modifies atmospheric conditions (temperature, precipitation, clouds, atmospheric dynamics) [1]. The modification by atmospheric aerosols of radiation and clouds (among other meteorological variables) represents a key topic when characterizing the climate of the Earth. The potential impacts of aerosol-radiation-cloud interactions may include (e.g., [2][3][4][5]): (1) exceedances of the aerosol critical levels causing changes in the patterns of extreme events; (2) high variation of levels of air pollutants; (3) changes in atmospheric dynamics, etc. Therefore, air quality problems interact (presenting both positive and negative feedbacks) with regional and global climate systems, so they must be evaluated within an integrated framework [6]. However, given the complexity and the lack of appropriate computer resources, traditional climate models strongly simplify the representation of aerosol-radiation-cloud interactions. This has historically led to two independent disciplines: the climate modelling and the air quality communities, creating separate modelling systems that are only loosely coupled (the so-called offline models) [7].
It is nowadays widely known that realistic simulations of the combined aerosol-radiation (ARI) and aerosol-cloud interactions (ACI) require models where aerosols, meteorology/climatology, radiation, clouds, and chemistry are coupled in a fully interactive manner. Otherwise, offline models could lead to large inconsistencies. For instance, the constant droplet number assumed by offline models in order to close the equations describing cloud droplet radiation interactions forces the use of empirical relationships between aerosol number and the cloud droplet nucleation term [8]. The usage of such relations may induce inconsistencies between cloud mass and droplet number in a model [9]. Hence, fully coupled (online) climate/chemistry models must be used nowadays in order to reduce (or, at least, to characterize) the uncertainties in the representation of aerosol-cloud processes.
Thus, motivated by a better understanding of the impact of ARI and/or ACI on a regional scale, online chemistry/climate models with full couplings between the air quality and meteorological/climatological components have been developed and applied recently (e.g., [10][11][12]). A growing number of studies have recently started to address the coupling of aerosols and clouds in regional chemistry/climate models, coordinated through collaborative European projects from FP7 and Horizon 2020 (e.g., MACC), International programmes such as AQMEII [13][14][15][16][17][18]; or COST Actions like ES1004-EuMetChem [12,19,20]. However, as observed in the scientific literature, the application of regional chemistry/climate on-line coupled models is currently biased towards studies of isolated episodes (e.g., [21][22][23][24]) because of the demanding computational resources. Hence, previous works do not allow a representation of aerosol-cloud interactions for climate-representative time periods.
Therefore, the main objective of this work is to evaluate the added value provided by the inclusion of ACI with respect to a base case (where only ARI effects are taken into account) during a period covering the years 2003-2010. For that, a fully coupled (online) climate/chemistry model (WRF-Chem) allowing ARI and ACI effects has been used.

Model Configuration
The regional climate simulations were carried out using the Weather Research Forecast (WRF) version 3.6.1 [25] and WRF model coupled with Chemistry (WRF-Chem) version 3.6.1 [26]. This model is widely used by the scientific community. For the historical period (2003-2010), the ERA20C reanalysis with a horizontal resolution of approximately 125 km [27] was used. The spatial Euro-CORDEX domain [28] covers Europe with a spatial resolution of 0.44 • both in latitude and longitude. The meteorological boundary conditions were updated every 6 h [27], while aerosol/gas boundary conditions were taken from MOZART-4/GEOS-5 [29]. Outputs were recorded every hour. Twenty-nine inhomogeneous sigma levels were used vertically. The top level was defined at 50 hPa. The evolution of CO 2 , CH 4 , and N 2 O greenhouse gases were considered according to the recommendation of Jerez et al. [30]. The physics configuration consisted of: the Lin microphysics scheme [31], the Noah land surface model [32], the RRTM radiative scheme [33], the Grell 3D ensemble cumulus scheme [34,35] and the University of Yonsei boundary layer scheme [36].
The boundary conditions from the global model do not provide information on the Saharan dust. For this reason, we implemented two one-way nested domains. An outer 150 km spatial resolution covered the most important areas of Saharan dust emission [41][42][43]. As the task of this domain is to provide the boundary conditions, spectral nudging was used in the parent domain so that the introduction of dust coincides with the synoptic situation. The nudging is applied between the top of the model and the planetary boundary layer to wind (vorticity and divergence), water vapour mixing ratio and temperature. The internal domain covers all of Europe and is established by Euro-CORDEX, as described before.

Aerosol Optical Depth (AOD) Simulations in WRF-Chem
WRF-Chem model run with the GOCART aerosol module [44]. The simulated aerosols included five species, namely sulphate, mineral dust, sea salt aerosol, organic matter, and black carbon. Aerosol Optical Depth (AOD) has been externally prognosticated at 550 nm (AOD550) for different chemical species, including sulphate, mineral dust, sea salt aerosol, organic matter, and black carbon. AOD550 obtained from these species is denoted hereafter as AODSU, AODDU, AODSS, AODOM, and AODBC, respectively. Total AOD is abbreviated as AODTO. The GOCART scheme was coupled with RACM-KPP [45] as the gas-phase chemistry option; the Fast-J photolysis module was used as photolysis option [46]. Biogenic emissions were calculated using the Guenther scheme [47]. AOD550 was estimated using the reconstructed mass-extinction method [48,49], which has also been applied in other works covering Europe (e.g., [50,51]). The sum of the absorption (β a ) and the attenuation by scattering (β s ) gives the extinction coefficient (β ext ), which is a function of the wavelength (λ). The total AOD550 is the vertical integration (in our case, the addition) of the product of the extinction and the layer thickness (∆Z i ) for each i layer, that is: where β ext (m −1 ) is estimated using the following semiempirical approach depending on the aerosol mass and the relative humidity (RH): Two sets of climatic simulations were run: (1) including only aerosol-radiation interactions (ARI), which was considered to be the base case in our study; and (2) simulations including aerosol-cloud interactions (ARI+ACI, denoted afterwards as ACI for the sake of brevity). WRF-Chem aerosol-cloud-radiation interactions are explained in detail in Palacios-Peña et al. [12]. For ARI [46,52], each chemical constituent of the aerosol was associated with a complex index of refraction. The overall refractive index for a given size bin was determined by volume averaging. The Mie theory and the summation over all size bins were used to determine the composite aerosol optical properties. Wet particle diameters were used in the calculations. Finally, aerosol optical properties are transferred to the shortwave radiation scheme. ACI in WRF-Chem [52] were implemented by linking the simulated cloud droplet number with the microphysics schemes. Therefore, the droplet number affected both the calculated droplet mean radius and the cloud optical depth.
It should be noted that the WRF-Chem version used (3.6.1) does not allow a full coupling of aerosol-cloud interactions. For instance, the model distinguishes wet scavenging for large-scale and sub-grid stratiform and sub-grid convective clouds. Convective wet scavenging (conv_tr_wetscav) is not available in WRF-Chem when using GOCART. Moreover, in-cloud wet scavenging or cloud chemistry is also not available. So, the full description of ACI is not implemented in WRF-Chem. However, the microphysics implemented in the simulations rely on the Lin scheme [31,53], a single moment scheme that turns into a two moments scheme in the simulations denoted as ACI. Further detail on ACI implemented in the simulations are detailed below.
As aforementioned, the Lin scheme is a single moment scheme based on Lin et al. [31], including some modifications, such as saturation adjustment [54] and ice sedimentation, which is related to the sedimentation of small ice crystals [55]. It includes six classes of hydrometeors: water vapour, cloud water, rain, cloud ice, snow, and graupel. This scheme was one of the first to parameterize snow, graupel, and mixed-phase processes (such as the Bergeron process and hail growth by riming), and it has been widely used in numerical weather studies.
The one-moment microphysical scheme is unsuitable for assessing the aerosol-clouds interactions as it only predicts the mass of cloud droplets and does not represent the number concentration of cloud droplets [56]. The the prediction of two moments provides a more robust treatment of the particle size distributions, which is key for computing the microphysical process rates and cloud/precipitation evolution. Therefore, prediction of additional moments allows greater flexibility in representing size distributions and hence microphysical process rates.
In this sense, although the Lin microphysics is presented as a single moment scheme, the WRF-Chem model allows to transform the single into a double moment scheme. A prognostic treatment of cloud droplet number was added [57], which treats water vapour and cloud water, rain, cloud ice, snow, and graupel. The autoconversion of cloud droplets to rain droplets depends on droplet number [58]. Droplet-number nucleation and (complete) evaporation rates correspond to the aerosol activation and resuspension rates. Ice nuclei based on predicted particulates are not treated. However, ice clouds are included via the prescribed ice nuclei distribution following the Lin scheme. Finally, the interactions of clouds and incoming solar radiation have been implemented by linking simulated cloud droplet number with the Goddard shortwave radiation scheme, representing the first indirect effect, and with Lin microphysics, which represents the second indirect effect [25]. Therefore, droplet number will affect both the calculated droplet mean radius and cloud optical depth.

Copernicus Atmosphere Monitoring Service (CAMS) Reanalysis
The estimations of AODSU, AODDU, AODSS, AODOM, AODBC, and AODTO (at 550 nm) were compared to the monthly information provided by Copernicus Atmosphere Monitoring Service (CAMS) reanalysis [59]. This reanalysis is the latest global reanalysis dataset of atmospheric composition produced by the European Centre for Medium-Range Weather Forecasts (ECMWF), consisting of three-dimensional time-consistent atmospheric composition fields, including aerosols and chemical species with a horizontal resolution of 80 km. This reanalysis has been extensively evaluated for aerosols and gas-phase pollutants against observations during the years 2003-2016 [59,60], finding an error generally ranging from 10% to 20%. Focusing on AOD, CAMS reanalysis improves the skill of other existing chemistry reanalyses when compared against AERONET observations, indicating a very slight negative monthly bias of −0.007 over Europe [59].

AErosol RObotic NETwork (AERONET)
In addition to CAMS, daily data from the AErosol RObotic NETwork (AERONET) stations available during the target period (2003-2010) have been used to evaluate the AOD representation by CAMS and regional simulations (ARI-basecase and ACI). AOD at 675 nm (AOD675) is the evaluated variable because this wavelength is provided by a large number of stations (data is scarce at 550 nm). Model data (AOD550) has been estimated at 675 nm by using the Ångström exponent. Figure 1 shows AERONET stations represented by a number. The name, as well as the coordinates of the stations used, are listed in Table 1. AERONET observations come from Level 2.0 and their error is (under cloud-free conditions) < ±0.01 for λ > 440 nm and < ±0.02 for shorter wavelengths [61].  Table 1, where stations are listed by code longitude, latitude, name, and corresponding number in this figure.

Previous Validation Results
The results of the modelling system presented here have been extensively evaluated in previous works. In this sense, the simulations have been compared to the remote-sensing data from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors aboard the Terra and Aqua platforms (MOD04_L2 and MYD04_L2) from collection six (C6) [18], finding that large underestimations in the simulations can be related to the misrepresentation of black carbon coming from biomass burning episodes. Also, the slight underestimations found for AODTO can be ascribed to a misrepresentation of the aerosol vertical profile modelled in these simulations [12].
The evaluation of meteorological fields indicate that maximum and mean temperatures were underestimated in the simulations [19] when compared to the E-OBS database [62]. In addition, these authors found no straightforward conclusion with respect to the improvement (or not) of the bias when introducing ACI interactions with respect to ARI results (bias of −0.696 K in mean temperature in ARI simulations and −0.642). However, very important improvements were found for the spatiotemporal representation of the temperature when ACI simulations were compared to ARI.
Regarding cloud properties, Baró et al. [20] found that an underestimation(overestimation) of cloud fraction (CF) is observed over land(sea) areas. The model simulations present a positive bias close to 40% over the sea, while a negative bias over −35% is over the land. While the negative bias may be due to the general underestimation in the cloud condensation nuclei representation by the models (that can be slightly corrected when introducing ACI in the simulations) [63], the overestimations found offshore might be related to satellite retrievals missing thin clouds [20,64]. Also, these authors found an underestimation of the cloud liquid water path (CLWP) related to the lower number of cloud droplets provided by the Lin microphysics schemes when compared to other microphysics schemes (e.g., [15,65]).  (29,Switzerland). This is an expected fact since these are background stations, located principally in rural areas where lower aerosol loads from anthropogenic origin are expected. A particular station is Bordeaux (27, France) with mean AOD675 values of 0.04 as a consequence of the high precipitation frequency fluxes over this area [66]. On the other hand, the highest AOD values (above 0.20) were found in Tarbes (30, France); Gloria (85, Romania); and Rossfeld (45, Germany) because these areas present an important anthropogenic pollution, especially associated to black carbon [39,67]. Medium AOD values are found in the many European cities as Belsk (40, Poland); Kyiv (19, Ukraine); Mainz, Leipzig, and Karlrsuhe (7, 26, 52, Germany), Paris metropolitan area (10, 12, 13, France); southwest of Great Britain and Benelux cities and all the cities bordering the Mediterranean. These levels could also be attributed to high aerosol loads from anthropogenic origin, since these areas are far from natural influences. However, in more southern locations like the Mediterranean, AOD675 levels could also come from recurrent desert dust outbreaks [68,69].

AOD AERONET Evaluation
When mean bias error (MBE) and mean average error (MAE) are analyzed (Figure 2b,c) there is a general slight overestimation. Only 13% of the stations used underestimate AOD675 for the CAMS reanalysis. The skill of the model improves in regional simulations (ARI and ACI) with respect to CAMS reanalysis for both statistical figures. Temporal and spatial MBE mean is 0.023 for CAMS, while the bias error reduces to 0.015 for both regional simulations. On the other hand, MAE mean is 0.059 for CAMS, 0.048 for ARI (base case) and 0.049 for ACI. According to these results, the ARI and ACI simulation present a significant improvement in the absolute error with respect to the global CAMS reanalysis, while no significant differences are found between ARI and ACI.
The aforementioned analysis is representative of the average behavior of most of the stations, but some of them should be highlighted. In Bordeaux (27,  The fact that bias and error are identical in absolute value means that the bias is always produced in the same direction (always under-or overestimation) for each station. A small number of stations exhibit an MBE for regional coupled simulations higher than CAMS (3% and 31% of the stations), but these differences are usually negligible (±0.01).
With respect to the time correlation coefficient (r), the time series have been normalized before estimating the correlation (that is, the annual cycle has been removed, so that the correlation is not influenced by the seasonality of the signal). Results indicate that there is not a clear behaviour that can lead to the definition of some added value of regional simulations over CAMS or ARI vs. ACI (Figure 2d). More than 50% of AERONET stations display low correlation values with modelled data. However, correlations are higher for regional climate simulations with WRF-Chem (ARI and ACI). Around 10% of stations show values higher than 0.5 when CAMS is evaluated while 15% of the stations do for ARI and 13% for ACI. Figure 3 depicts the bias error of the different modelled AOD550 climatologies when assessed against the CAMS reanalysis, the bias of the ARI-base case simulations and the improvement of the error with respect to the base case simulation when including ACI. The improvement of the error is defined as:

Biases and Improvement of the Absolute Errors in the Simulations Including Aerosol-Cloud Interactions
where s denotes the AOD at 550 nm for each chemical species s (total AOD, sulphate, dust, sea salt, organic matter, or black carbon). The AOD provided by CAMS reanalysis is taken in this section as the benchmark to compare our simulations. Henceforth, Equation (4) reflects that a positive number is indicative of a higher error in the ARI-base case simulation; that is, the error against CAMS is reduced when introducing ACI in our simulations. Conversely, negative values of IE indicate worse skills when aerosol-cloud interactions are included in regional simulations. The evaluation results compare the climatologies taken from CAMS reanalysis (Figure 3a) with the base case simulation. In general, the bias error of the simulations (Figure 3b) show an excellent agreement between all the species and CAMS-derived AOD variables, but with a general underestimation of AODTO (−7.62%) over the entire domain. This result is in agreement with those of previous works [20,70]. The bias for AODTO is higher in northern Africa, caused by the pervasive underestimation of the aerosol load coming from the Sahara desert (AODDU is underestimated by 0.02-0.04 in this area, representing an underestimation of over 20% of AOD). However, AODDU in northern Europe is slightly overestimated, especially over the northeastern part of the domain, caused by the unrealistic high values of dust coming from the Gobbi desert through the boundary conditions [12,18]. AODSU presents a slight underestimation (Table 2) both in ARI (−0.64 × 10 −3 ) and ACI simulations (−0.55 × 10 −3 ), especially over eastern Europe, which can be related to a strong underestimation of winter sulphate aerosol [71,72]. Similarly, AODSS and AODBC depict a tendency towards underestimation in ARI-base case simulation (−4.79 and −1.89 × 10 −3 , respectively), while AODOM is overestimated in ARI simulation (3.21 × 10 −3 averaged over the domain).
With respect to the improvement of the error (IE), our results in Figure 3 and Table 2 indicate a positive value, meaning that a reduction in the error is found for ACI simulations for all chemical species. ACI have a strong effect on particulate matter and AOD coming from different species, also conditioning the emissions of natural aerosols (sea salt, mineral dust, which strongly depend on wind speed) [73,74], atmospheric transport, and chemistry of large emitting sources, such as plumes from forest fires and cities [75]. It is noticeable, the improvement of AODDU, the error of which improves by 0.05 × 10 −3 , that is, an improvement of 20.22% in the error of ARI-base case simulations. This improvement, as discussed later, is especially caused by the reduction of the error over northern Africa. Also, AODSU reduces the error by 14.06%, as a consequence of an enhanced oxidation of sulphate (higher concentrations of hydroxyl radical) when including ACI. For the rest of the chemical components, inclusion of ACI reduces the error by lower than 5% (2.09% for AODSS, 3.42% for AODOM, and 2.11% for AODBC).  Figure 4a shows the time correlation coefficients (r) between the base case simulation (ARI) and CAMS reanalysis. It should be pointed out that, as done in Section 3.1, the annual cycle of the series have been removed before estimating the r. In general, the best correlation coefficients are estimated for the southern and south-eastern part of the target domain, with r > 0.80 for all chemical species, while the poorest correlations are calculated over the Scandinavian peninsula and the northeastern part of the domain, where r < 0.1 for several species as AODOM and AODDU. The r coefficients, averaged for the entire domain (Table 3)  Moreover, the lowest correlation coefficient in the case of AODOM is caused by the limited skill of the model to reproduce the optical depth associated to organic aerosols in the northeastern part of the domain (r < 0.1); this is also where the model presents the highest bias errors for this species. Overseas, the AODOM ranges between 0.50 and 0.70, pointing out to a good representation of the time variation of AODOM over the Mediterranean and the Atlantic Ocean. Regarding the improvement of the correlation coefficient when including ACI in the simulations (Figure 4b), the behaviour is similar to the case for the IE. The total correlation coefficient improves by 8%, averaged for all the domain, with slight decreases in the values of r over the Netherlands and nortwestern Germany. The largest improvement in the simulation is found for AODDU, where r increases by nearly 20% when including aerosol-cloud interactions (Table 3), mainly because of the improvements of the correlation coefficients over western and southern Europe (slight decreases of the r are found over the central part of the simulation domain). There is a better correlation in the rest of the chemical components when ACI are introduced in the models, ranging from an improvement of 6% for AODBC and AODSU to 9% for AODSS and AODOM. The better correlation coefficients found do not present a homogeneous behaviour over the entire domain. For instance, for AODOM strongly increases over the British Isles, while decreases over the Iberian Peninsula. For AODBC, it is strongly noticeable that the largest improvements are found over those areas with the largest emitting sources (central and eastern Europe), with increases of 0.10 over those areas where the correlation coefficient is lower (0.20-0.30), thus showing improvements of 50% in the correlation coefficient.

Discussion
The base case simulations (ARI) for the period 2003-2010 underestimate AODTO by 4% when compared with CAMS reanalysis (AODTO of 0.207 in CAMS, bias error of 0.008), which is in the same range as the bias found in other works. For instance, Baró et al. [20] found that, over Europe, moderately to largely underestimations are expected for aerosols loadings and AOD for most simulations from regional chemistry transport models. In this same sense, Lapina et al. [76] found out that GEOS-Chem underestimates the mean AOD value for the marine environment by 21%.
Other works (such as Balzarini et al. [77]) indicate that, albeit different models trust different chemical mechanisms for reproducing aerosol optical properties, the AODTO is generally underestimated. These discrepancies could be studied in further works through sensitivity analysis, but the general trend for underestimations of AODTO in WRF-Chem simulations can be assigned to underpredictions in the aerosol dry mass by the models, an underestimation of the fraction of particles for a given mass, or to a misleading representation of the water associated with aerosols [69].
Focusing on the comparison with remote sensing observations, Pozzer et al. [70] compared AOD climate simulations to AERONET observations, finding that biases vary between 0.02 and 0.03. Hence, regional online coupled models improve those values, with AOD biases of 0.015 for both ARI and ACI simulations, a bias 50% lower than that found for global models as a consequence of the higher resolution used in regional models. In addition, Palacios et al. [18] used an ensemble of models from AQMEII Phase 3 [17] to assess AODTO against in-situ and remote sensing observations, finding out that the biases of the ensemble of AQMEII simulations when representing AOD are in the order of 0.05. In this sense, the bias observed in our evaluation is again lower than those found in the literature, pointing out to an accurate skill of the model when reproducing aerosol climatologies over the European area. However, this improvement is not that evident for the correlation coefficient, which does not show a clear spatial pattern. This fact could be ascribed to the use of a monthly temporal resolution for the estimation of r, reducing the number of occurrences and hence hampering the estimation of an accurate correlation.
Apart from AODTO, scientific literature is much more limited when assessing chemical components of the AOD. In the case of AODSU, the underestimation found in ARI-base case could be related to two causes: (1) the pervasive tendency of regional models to be geographically biased by strong winter underestimations of sulphate concentrations at eastern European locations (it is widely known that the winter underestimation of sulphate is a common issue in most models that operate over Europe, which represent a direct couplet of sulfur chemistry with photochemistry) [71]; and (2) the fact that sulphate formation in the modelling systems is often limited by oxidant availability [72]. With respect to the AODOM, the inclusion of a complex organic aerosols mechanism in models (such as the mechanism used here) leads to a better agreement over certain areas over the sea (for instance, over the Atlantic ocean, as observed here in our simulations) but leads to overestimates of OM in other regions [76].
With respect to the improvement of the simulations when including ACI, there is a general improvement in the error of AODTO and correlation coefficients (around 8% for AODTO, over 20% in the case of AODDU) when compared with CAMS. Since the model generally underpredicts the CAMS concentrations, the inclusion of aerosol-cloud interactions can have an important effect increasing relative humidity [24], modifying the size of the particles due to hygroscopic growth, and hence, the AOD predictions, especially for hydrophilic species.
The improvement of the errors found here for the improvement of the AOD during the 2003-2010 period agrees with the results of other works for episodic studies (e.g., [18,19,23]). These works indicate that minor, but statistically significant, improvements are observed in AODTO and other climatic variables such as 2-m temperature when ACI processes and feedbacks are taken into account. For instance, improvements in AODTO around 12% during an episode of the Saharan dust outbreak over Europe in October 2010 when including ACI in an ensemble of models run under the EuMetChem initiative [12,19]. Moreover, these latter authors found an improvement of the error of 11% when ACI were included in WRF-Chem during the Russian 2010 heatwave episode (July 2010), even over areas far from wildfires emissions [12].
Hence, the IE is especially noticeable in the case of natural aerosols (dust and sea salt, of which emissions strongly depend on wind speed), as indicated in previous scientific literature [78,79]. As seen in Figure 3, there is a strong underestimation of AODDU over the southern part of the domain. ACI simulations provide a slight increase of wind speed over the Atlantic and southern part of the domain. This is a striking feature, found not only in our simulations but in other WRF-Chem simulations, including ACI (e.g., [19,24]), as a consequence of the strong increase in solar radiation over those target areas. Simulated cloud droplet numbers are between 5 and 100 cm 3 when ACI are taken into account over remote areas as the Atlantic ocean. If simulated aerosol concentrations are not taken into account for cloud droplet formation (as this is the case of our ARI-base case simulations), a standard value of 250 cloud droplets per cm 3 is assumed for the entire modelling domain, coming from the standard configuration of WRF and WRF-Chem [16]. Hence, the reduced number of cloud droplets over the Atlantic reduced cloudiness, increases solar radiation and the higher solar radiation enhances temperature and therefore winds, which promote the emission of dust and sea salt and therefore increase the AODDU and AODSS in our ACI simulations, reducing the negative bias found for dust over northern Africa and for sea salt over the entire domain.
Finally, it should be highlighted that the results found here may point to some benefits of including ACI when representing aerosol climatologies over Europe. The improvement of the skills in the simulations (around 8% for AODTO, reaching 20% for dust bias) could make climate scientists consider the inclusion of ACI in climate simulations, despite the larger computation required (higher by a factor of eight with respect to the base case simulations) and considering the feedbacks of air pollution in the calculation of meteorological fields in coupled chemistry/climate regional models.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

ACI
Aerosol-cloud interactions AOD Aerosol optical depth AOD550 Aerosol optical depth at 550 nm AOD675 Aerosol