Assessing Desert Dust Indirect E ﬀ ects on Cloud Microphysics through a Cloud Nucleation Scheme: A Case Study over the Western Mediterranean

: In this study, the performance and characteristics of the advanced cloud nucleation scheme of Fountoukis and Nenes, embedded in the fully coupled Weather Research and Forecasting / Chemistry (WRF / Chem) model, are investigated. Furthermore, the impact of dust particles on the distribution of the cloud condensation nuclei (CCN) and the way they modify the pattern of the precipitation are also examined. For the simulation of dust particle concentration, the Georgia Tech Goddard Global Ozone Chemistry Aerosol Radiation and Transport of Air Force Weather Agency (GOCART-AFWA) is used as it includes components for the representation of dust emission and transport. The aerosol activation parameterization scheme of Fountoukis and Nenes has been implemented in the six-class WRF double-moment (WDM6) microphysics scheme, which treats the CCN distribution as a prognostic variable, but does not take into account the concentration of dust aerosols. Additionally, the presence of dust particles that may facilitate the activation of CCN into cloud or rain droplets has also been incorporated in the cumulus scheme of Grell and Freitas. The embedded scheme is assessed through a case study of signiﬁcant dust advection over the Western Mediterranean, characterized by severe rainfall. Inclusion of CCN based on prognostic dust particles leads to the suppression of precipitation over hazy areas. On the contrary, precipitation is enhanced over areas away from the dust event. The new prognostic CCN distribution improves in general the forecasting skill of the model as bias scores, the root mean square error (RMSE), false alarm ratio ( FAR ) and frequencies of missed forecasts ( FOM ) are limited when modelled data are compared against satellite, LIDAR and aircraft observations.


Introduction
Dust particles cover approximately 50% of the global load of the natural airborne particles [1,2]. Dust aerosols interact with electromagnetic radiation [3] and directly affect the radiation budget of the land-atmosphere system (aerosol direct effect) [4]. Furthermore, they absorb and scatter in various wave lengths both solar and long wave radiation [5][6][7][8]. Consequently, they modify the heating rates, the stability and the stratification of the atmosphere [9,10]. Additionally, dust interactions with radiation

The Embedded Cloud Nucleation Scheme
The WRF double-moment (WDM6) microphysics scheme [13] supports prognostic equations for water substance variables, such as cloud water, water vapor, rain, cloud ice, graupel and snow. Furthermore, it calculates the prognostic concentration number of CCN, cloud and rain water. The concentration number of ice species like graupel, snow and ice follows the diagnostic estimation of the ice phase process of Hong [70]. This simple estimation is theoretically accepted as the prediction of the ice phase number concentrations has less influence than the prediction of warm phase concentrations in deep convective cases [71]. Warm rain processes such as accretion and autoconversion are based on the methodology of Cohard [72]. The calculation of any other source/sink terms is based on the WRF single-moment (WSM6) microphysics scheme [61].
The prediction of the activated CCN concentration number uses the drop activation process based on the Twomey relationship, which links the activated CCN number and the supersaturation [20,73]. As the evaporation of cloud droplets is also considered, the corresponding CCN number is added to the total CCN concentration number. Any other sink/source mechanism apart from the CCN activation and evaporation is ignored. However, no real time CCN concentration is supported. Thus, a constant initial CCN number of 100 cm −3 [67] is used, in order to treat the CCN spatiotemporal distribution as a prognostic variable. The above value is representative of maritime aerosol particles, which can act as effective CCN [13]. The activated CCN number concentration is predicted using the drop activation by Twomey's relationship between the number of activated CCN and supersaturation: n act = (n + N C ) S S max k (1) Remote Sens. 2020, 12, 3473 4 of 25 where n act is the activated CCN, n is the initial CCN value, Nc is the cloud droplet number concentration, S is the supersaturation and k (0.6 is set in the current study) is the parameter following the observation varying from 0.3 to 1.0 [73]. The evaporation of cloud droplets is also taken into account and returns CCN particles to the overall CCN count. Any other sink/source term is ignored. However, the above approximation remains the weakest link in this scheme, as there is no real time CCN replenishment [13].
The microphysics processes of WDM6 differ significantly from those of WSM6 due to the predicted concentration number of CCN, cloud (N C ) and rain (N R ) droplets. All microphysics processes are calculated according to the following cloud-rain drop distribution: where X ∈ [C, R] indicates clouds or rain, ν x and α x are the two dispersion parameters. N X and D X are the predicted total concentration number and diameter of the X drop category, respectively. The slope parameter λ X is given by: where v R = 2 and α R = 1 for rain and ν c = 1 and α c = 3 for clouds. Moreover, q x (kg kg −1 ) is the mixing ratio for the X hydrometeor, ρ α (kg m -3 ) and ρ w = 10 3 (kg m −3 ) are the air and water density, respectively. In this study, a prognostic initial number of CCN is applied, using the dust aerosol mass concentrations of the Georgia Tech Goddard Global Ozone Chemistry Aerosol Radiation and Transport of Air Force Weather Agency (GOCART/AFWA) emission scheme of the WRF/Chem model. The number concentration of dust particles N (kg −1 ) for an n size bin is calculated as: where C is the dust mass mixing ratio (µg kg −1 ) for an m size bin. Summing up the N for all size bins (information concerning the bins is given in the following Section 3) gives the total concentration number of the dust aerosols. The calculation of the prognostic initial concentration number of CCN is based on the parameterizations of the cloud droplet activation scheme of Fountoukis and Nenes [55], which is an evolution of the study of Nenes and Seinfeld [74]. This scheme (hereafter FN) represents the aerosol number and their chemical composition with respect to the size and the calculation of the initial CCN number (the CCN spectrum) at a specific supersaturation (s) level based on Köhler theory [75]. Furthermore, the CCN spectrum exists in the dynamical conditions of an air parcel with constant updraft velocity. The CCN spectrum F s is given by: where N i is the aerosol number concentration, s is the supersaturation, s g,i the critical supersaturation, σ i is the geographic standard deviation, i is the bin number. Processed dust is assumed to have the same hygroscopic properties as sulfate. Moreover, the prognostic calculation of the cloud droplet concentration number is achieved through the computation of the maximum supersaturation (s max ) as the air parcel ascends. With s = s max , the cloud droplet concentration number is given by: Remote Sens. 2020, 12, 3473

of 25
The calculation of the s max is based on a parameterization of the water vapor balance. Full information on how the s, the s g,I and the s max are calculated can be found in the studies of Nenes and Seinfeld [74] and Fountoukis and Nenes [55]. It should be noted that CCN parameterizations of the original WDM6 scheme have not been deactivated in order to take into account the impact of sea salt particles.
As cumulus parameterization schemes are responsible for the simulation of convection, in this study, the CCN treatment by FN scheme has also been embedded in the Grell-Freitas (GF) cumulus parameterization scheme [52] in order to estimate the conversion rate and the precipitation efficiency on the CCN number concentration. The scheme includes modulation of convection by the environment (dynamic control), modulation of the environment by the convection (feedback) and a cloud model that defines the cloud properties (static control). Convection is triggered when the CAPE values are positive and the difference between cloud base (Lift Condensation Level LCL) and cloud top is 150 hPa minimum. Moreover, the LCL has to be within 50 hPa of the Level of Free Convection-LFC [76]. The original Grell and Freitas [52] cumulus scheme uses fixed values for CCN concentrations and for the conversion rate of condensation to precipitation. Thus, in this study, the prognostic CCN concentrations calculated by the implemented FN scheme are used. Furthermore, the above scheme is chosen as it is a scale aware scheme that can also be activated in the so-called "grey zones" concerning the spatial resolutions (6-10 km). Recent studies have shown that the activation of the GF scheme also in the so-called "grey zones" improves the overall precipitation pattern by reducing biases and improving root mean square error values [52,77].

Methodology
The simulations are performed using the fully coupled meteorology-chemistry Weather Research and Forecasting/Chemistry model (version 3.9) WRF/Chem [56,57]. Both components use the same horizontal and vertical grids, the same physics schemes for the sub-grid scale transport and the same transport for the preserving of the air and scalar mass. Hence, the meteorological field is simultaneously simulated with the chemical field and the interactions in between. The chemical component is based on the Georgia Tech/Goddard Global Ozon Chemistry Aerosol Radiation and Transport (GOCART) module, which simulates tropospheric aerosols, such as sulfate, mineral dust, sea salt and black and organic carbon [78,79]. As far as the dust emissions are concerned, the Air Force Weather Agency (AFWA) package is used as it supports a non-experimental and sectional dust emissions scheme using an advanced set of dust emissions parametrizations [80,81], including five dust bins of 0.7, 1.4, 2.4, 4.8 and 8.0 µm effective radii. The scheme is based on a modified version of the saltation function of Marticorena and Bergametti [82]. The saltation for large-scale particles is triggered by wind shear and leads to fine particle emissions through bombardment [83]. There is a proportionality between horizontal and vertical dust flux, which is calculated when the friction velocity exceeds a specific threshold. The horizontal dust flux is calculated according to White [84]. The above proportionality is based on clay sand and erodibility fields. However, some studies [85] estimate that the linkage of the parametrized vertical dust flux to the clay fields is the main reason for the systematic overprediction of the simulated dust concentration. Furthermore, a correction to soil moisture parameterizations is applied [86]. The distribution of the emitted dust particles is based on the brittle fragmentation theory [87].
The GOCART/AFWA scheme supports by default a dust dry deposition scheme based on [88]. The calculation of dry deposition velocity is based on the aerodynamic, the quasi-laminar layer and the surface resistances. As far as wet deposition is concerned, the scheme of Seinfeld and Pandis [75] is implemented for both in and below cloud scavenging. It is shown in literature that especially the in-cloud area plays a leading role in the scavenging of suspended particles affecting the entire spatiotemporal distribution of the dust particles [89]. However, indirect effects concerning dust aerosols with cloud interactions are not supported. Hence, as is described in Section 2.1, the nucleation scheme of Fountoukis and Nenes [55] is embedded in the WDM6 [13,67] double-moment microphysics scheme. Furthermore, the fixed CCN concentrations used in the Grell-Freitas (GF) cumulus parameterization scheme [52] are replaced by the overall CCN concentrations calculated by the Fountoukis and Nenes [55] cloud nucleation scheme.
The performance of the Fountoukis and Nenes [55] scheme is assessed through a comparison of modelled precipitation against the WORLDVIEW-NASA satellite retrievals (Atmospheric Infrared Sounder-AIRS-Level 2). AIRS is an instrument based on the second Earth Observing System polar orbiting form (EOS-Aqua) and includes a group of infrared, visible and microwave sensors. The AIRS precipitation is an estimate of daily precipitation (in mm), with a 45 km spatial resolution with coverage twice a day. Parameters of cloud top pressure, fractional cloud cover and cloud-layer relative humidity have been used. Satellite data have been used as the amount of ground stations data for the specific areas were insufficient. The nearest model grid point values to the satellite ones have been used for the evaluation. A contingency table (Table 1) according to Wilks [90] methodology is used in order to calculate bias (B), false alarm ratio (FAR) and frequencies of missed forecasts (FOM). According to Table 1 a pairs mean that both observations and modeled values exceed a specific precipitation threshold. Similarly, b pairs mean that modelled values exceed the threshold but observation do not. Furthermore, c pairs mean that the model does not exceed the precipitation threshold; however, the observations are over it. Finally, d pairs mean that both model and observation did not manage to exceed the precipitation threshold Bias greater than one corresponds to the model's overprediction, whereas bias less than one indicates the model's underprediction. For unbiased forecasts, B is equal to one. For both FAR and FOM, best scores are equal to zero, whereas the value equal to one (or 100%) corresponds to the worst performance.
The response of the atmospheric optical properties to the explicitly resolved CCN patterns is further assessed through the use of LIDAR data of the EARLINET data base. The Aerosol extinction coefficient is used (550 nm-1.5 Level) instead of atmospheric optical depth (AOD), as the relationship between AOD and CCN is more complicated. More specifically, CCN concerning aerosol cloud interactions are located at the cloud area, whereas AOD is related to the entire atmospheric vertical column. Furthermore, the air mass may interact with clouds kilometers away or hours after the AOD measurements, which may occur under clear sky conditions [91]. Modelled CCN concentrations have also been compared against aircraft measurements (French aircraft service for environmental research-SAFIRE/MNPCA-CNRM/GAME) of ChArMEx/ADRIMED campaign (Chemistry-Aerosol Mediterranean Experiment/Aerosol Direct Radiative Impact on the regional climate in the MEDiterranean Region). The biases and the root mean square errors (RMSEs) concerning the above comparisons have been calculated by: Moreover, a correlation coefficient (R) is calculated as a square root of Pearson correlation coefficient: where M i and O i are the model estimates and the observed values respectively.

Configuration of the Modeling Experiments
In order to evaluate the performance of the embedded FN scheme in the WRF/Chem model, two experiments have been conducted, during a significant dust transport episode over the Western Mediterranean, characterized by torrential rainfall over Western Europe, during 19-25 June 2013. For this period and specific area, sulfate concentrations have been observed in the dust layer, during the ChArMEx/ADRIMED (the Chemistry-Aerosol Mediterranean Experiment/Aerosol Direct Radiative Impact on the regional climate in the MEDiterranean Region) campaign [92]. Aircraft measurements of the above study have indicated that dust particles were transported inside layers located above 1 km (atmospheric boundary layer), where pollution was dominant. Furthermore, measurements have also confirmed the presence of externally sulfate-coated particles, indicating that the mixing of dust particles with pollution can be justified and thus dust could be considered as effective CCN medium. The first numerical experiment is the control run (namely, CTRL), which does not account for any aerosol indirect effect. The second one (namely, FN run) considers the dust aerosol indirect effects through the FN scheme. The value of 5 cm s −1 (0.5 Pa s −1 ) has been used as a constant vertical velocity and taken from measurements performed during that study [92]. For both experiments, a cold-started simulation was initialized, on 16 June 2013, to prepare the dust background concentration. Both experiments, namely CTRL and FN, have been hot started on 20 June 2013 at 12 UTC in 120 h of simulations up to 25 June 12 UTC.
WRF/Chem is set up in three nested domains. The coarser domain (namely D01) covers North Africa, the Mediterranean, the Saudi Arabia peninsula and Europe, with a resolution of 18 × 18 km, which consists of 540 × 300 horizontal grid points and with 38 vertical points up to 50 hPa. The remaining two nested domains cover an area of Western Europe (namely D02 and D03, respectively). In this study a two-way nesting is used. D02 covers Western Europe with a resolution of 6 × 6 km, consisting of 451 × 295 horizontal grid points, while D03 covers Corsica and Sardinia on a resolution of 2 × 2 km having 286 × 145 horizontal grid points. The third domain concerns the area where aircraft measurements of CCN concentrations took place. The cumulus parameterization scheme is deactivated and convection is resolved explicitly in D03. Simulated CCN concentrations over this domain will be compared against observational records in order to assess the performance of the embedded FN scheme. The three domains are shown in Figure 1. The European Center of Mesoscale Forecast (ECMWF) operational analyses are used for both the initial and boundary conditions on a resolution of 0.5 • × 0.5 • . The Rapid Radiative Transfer Model (RRTMG) is used for the simulation of both long-and short-wave radiations [93]. The basic parameterizations concerning the simulations of this study are shown in Table 2.

Description of the 20-25 June 2013 Synoptic Conditions
There are several variables that characterize synoptic weather conditions favorable to the dust mobilization. The most important variables are the geopotential height (usually at 500 hPa), the sealevel pressure and the surface winds. Well organized troughs (formations of low geopotential heights) lead to atmospheric instability conditions and downward momentum transport. Hence, a spatial pressure field is formed which is directly linked to the generation of significant pre and post frontal winds responsible for the dust triggering and transport [97]. The event of 20-25 June 2013 has been chosen, as it was associated with torrential rainfalls over Western Europe caused by a lowpressure system persistence over the Iberian Peninsula. On 19 June 2014 at 00 UTC ( Figure 2a), a cutoff low over the Iberian Peninsula triggered baroclinic instability conditions. At that time, an extended frontal zone was located over Western France (Figure 2b). The prevailing SW winds ahead of the cold front at the tail of the frontal zone favored moisture advection from Western Mediterranean Sea and triggered torrential rainfall to Western France occurred on 21 June. This flow pattern also favored a large-scale desert dust transport from Algeria and Morocco and led to significant increase in dust load over Western Europe ( Figure 3). Indeed, the simulated dust load presents local maxima of 3000 mg m −2 over the Western Mediterranean on 21 June at 12 UTC ( Figure 4).

Description of the 20-25 June 2013 Synoptic Conditions
There are several variables that characterize synoptic weather conditions favorable to the dust mobilization. The most important variables are the geopotential height (usually at 500 hPa), the sea-level pressure and the surface winds. Well organized troughs (formations of low geopotential heights) lead to atmospheric instability conditions and downward momentum transport. Hence, a spatial pressure field is formed which is directly linked to the generation of significant pre and post frontal winds responsible for the dust triggering and transport [97]. The event of 20-25 June 2013 has been chosen, as it was associated with torrential rainfalls over Western Europe caused by a low-pressure system persistence over the Iberian Peninsula. On 19 June 2014 at 00 UTC (Figure 2a), a cut-off low over the Iberian Peninsula triggered baroclinic instability conditions. At that time, an extended frontal zone was located over Western France (Figure 2b). The prevailing SW winds ahead of the cold front at the tail of the frontal zone favored moisture advection from Western Mediterranean Sea and triggered torrential rainfall to Western France occurred on 21 June. This flow pattern also favored a large-scale desert dust transport from Algeria and Morocco and led to significant increase in dust load over Western Europe ( Figure 3). Indeed, the simulated dust load presents local maxima of 3000 mg m −2 over the Western Mediterranean on 21 June at 12 UTC ( Figure 4). off low over the Iberian Peninsula triggered baroclinic instability conditions. At that time, an extended frontal zone was located over Western France (Figure 2b). The prevailing SW winds ahead of the cold front at the tail of the frontal zone favored moisture advection from Western Mediterranean Sea and triggered torrential rainfall to Western France occurred on 21 June. This flow pattern also favored a large-scale desert dust transport from Algeria and Morocco and led to significant increase in dust load over Western Europe (Figure 3). Indeed, the simulated dust load presents local maxima of 3000 mg m −2 over the Western Mediterranean on 21 June at 12 UTC ( Figure 4).

Comparison of NCCN, Nc and RH
The average concentrations of CCN and cloud droplets of both CTRL and FN simulations as well as their differences (FN-CTRL) for 21 June 2013, 12 UTC, up to the height of 3 km and over D01 and D02 are shown in Figure 5. As expected, the CCN concentration (NCCN) of the FN run is in general much larger than the CTRL values, which are around 100 cm −3 (Figure 5a). On the other hand, FN simulation estimates that NCCN will reach 2000 cm −3 at the pathways of the severe dust advection over the Western Mediterranean (Figure 5b). The detected differences in NCCN are up to 1000 cm −3 over France (Figure 5c), showing a dominant hazy atmosphere at that time [91,98]. The total maxima of

Comparison of NCCN, Nc and RH
The average concentrations of CCN and cloud droplets of both CTRL and FN simulations as well as their differences (FN-CTRL) for 21 June 2013, 12 UTC, up to the height of 3 km and over D01 and D02 are shown in Figure 5. As expected, the CCN concentration (NCCN) of the FN run is in general much larger than the CTRL values, which are around 100 cm −3 (Figure 5a). On the other hand, FN simulation estimates that NCCN will reach 2000 cm −3 at the pathways of the severe dust advection over the Western Mediterranean (Figure 5b). The detected differences in NCCN are up to 1000 cm −3 over

Comparison of N CCN , N c and RH
The average concentrations of CCN and cloud droplets of both CTRL and FN simulations as well as their differences (FN-CTRL) for 21 June 2013, 12 UTC, up to the height of 3 km and over D01 and D02 are shown in Figure 5. As expected, the CCN concentration (N CCN ) of the FN run is in general much larger than the CTRL values, which are around 100 cm −3 (Figure 5a). On the other hand, FN simulation estimates that N CCN will reach 2000 cm −3 at the pathways of the severe dust advection over the Western Mediterranean (Figure 5b). The detected differences in N CCN are up to 1000 cm −3 over France (Figure 5c), showing a dominant hazy atmosphere at that time [91,98]. The total maxima of the N CCN differences are profoundly located over the dust sources (e.g., over Libya) areas reaching values of 4000 cm −3 . In general, the N CCN concentrations are higher over land areas, especially over the areas that are considered as important aerosol sources [99,100].
Remote Sens. 2020, 20, x FOR PEER REVIEW 10 of 26 The above results are comparable with those reported by [101,102], reporting that NC and NCCN maxima coincide as the last ones trigger the activation of more cloud droplets [25,98,102]. In the CTRL run (Figure 5e), however, NC presents a different to FN simulation pattern with the maxima of 300 cm −3 being located away from the areas influenced by the dust advection. This is because the overall NCCN and NC patterns in the FN run are altered by the mass preservation principle. Thus, the NC dominance in FN simulation over South Western France (Figure 5f) triggers the precipitation suppression 6h later as the rain water becomes limited. Nc has also decreased over Central Europe by 400 cm −3 approximately, as this area is not characterized by significant dust advection. More specifically, NCCN of the FN run do not exceed 100 cm −3 and the prospect of a significant formation of NC is actually small. The relative humidity also plays an important role in the above mechanism. Figure 5g,h show the simulated relative humidity at FN and CTRL runs. In both runs, high values of relative humidity coincide with high or elevated NC concentrations. Figure 5i shows the FN-CTRL relative humidity differences. As expected over the area of South Western France, the relative humidity differences reach 30-50%. This increase coincides with the increase in the NC concentrations over the specific area. The abovementioned result confirms that high relative humidity creates favorable conditions for the cloud droplet activation [66,103]. Over Central Europe, the relative humidity is decreased up to 30%, contributing to the NC decrease over that specific area.

Impact on the Precipitation Pattern
Τhe increase in the CCN concentration due to the presence of dust particles affects the cloud condensation and the precipitation pattern. Figure 6a shows the D02 simulated 6h convective rainfall for the FN run, whereas Figure 6b indicates the convective rainfall differences of the FN against the CTRL run and the simulated CCN concentrations of the FN run in contours. The FN convective rainfall is enhanced by up to 6mm over the marine areas of Western France and the Biscay Bay. This means that the prognostic treatment of the CCN embedded in the cumulus parameterization scheme  (Figure 5d). The above results are comparable with those reported by [101,102], reporting that N C and N CCN maxima coincide as the last ones trigger the activation of more cloud droplets [25,98,102]. In the CTRL run (Figure 5e), however, N C presents a different to FN simulation pattern with the maxima of 300 cm −3 being located away from the areas influenced by the dust advection. This is because the overall N CCN and N C patterns in the FN run are altered by the mass preservation principle. Thus, the N C dominance in FN simulation over South Western France (Figure 5f) triggers the precipitation suppression 6 h later as the rain water becomes limited. Nc has also decreased over Central Europe by 400 cm −3 approximately, as this area is not characterized by significant dust advection. More specifically, N CCN of the FN run do not exceed 100 cm −3 and the prospect of a significant formation of N C is actually small. The relative humidity also plays an important role in the above mechanism. Figure 5g,h show the simulated relative humidity at FN and CTRL runs. In both runs, high values of relative humidity coincide with high or elevated N C concentrations. Figure 5i shows the FN-CTRL relative humidity differences. As expected over the area of South Western France, the relative humidity differences reach 30-50%. This increase coincides with the increase in the N C concentrations over the specific area. The abovementioned result confirms that high relative humidity creates favorable conditions for the cloud droplet activation [66,103]. Over Central Europe, the relative humidity is decreased up to 30%, contributing to the N C decrease over that specific area.

Impact on the Precipitation Pattern
The increase in the CCN concentration due to the presence of dust particles affects the cloud condensation and the precipitation pattern. Figure 6a shows the D02 simulated 6h convective rainfall for the FN run, whereas Figure 6b indicates the convective rainfall differences of the FN against the CTRL run and the simulated CCN concentrations of the FN run in contours. The FN convective rainfall is enhanced by up to 6 mm over the marine areas of Western France and the Biscay Bay. This means that the prognostic treatment of the CCN embedded in the cumulus parameterization scheme enhances the convective rainfall over marine areas with lower CCN concentrations, which do not exceed the 50 cm −3 over the area of the Biscay Bay. In this case, the clouds rain out fast as large rain droplets are allowed to develop and the evaporation conditions are very limited [33,40,104]. On the contrary, in the FN run, a decrease in the convective precipitation by up to 8 mm is observed over the land areas of NE France, where a hazier atmosphere dominates as the CCN concentrations vary between 450 and 650 cm −3 . Clouds in such environments are characterized by significant evaporation rates before the occurrence of precipitation [33]. Thus, the autoconversion of cloud to rain droplets and consequently to precipitation is very limited [40,105].
Remote Sens. 2020, 20, x FOR PEER REVIEW 11 of 26 enhances the convective rainfall over marine areas with lower CCN concentrations, which do not exceed the 50 cm −3 over the area of the Biscay Bay. In this case, the clouds rain out fast as large rain droplets are allowed to develop and the evaporation conditions are very limited [33,40,104]. On the contrary, in the FN run, a decrease in the convective precipitation by up to 8 mm is observed over the land areas of NE France, where a hazier atmosphere dominates as the CCN concentrations vary between 450 and 650 cm −3 . Clouds in such environments are characterized by significant evaporation rates before the occurrence of precipitation [33]. Thus, the autoconversion of cloud to rain droplets and consequently to precipitation is very limited [40,105]. To investigate the impact of the FN scheme to the non-convective rainfall, two areas in the D02 are selected, namely Areas 1,2 (Figure 7), as in these areas, significant suppression or enhancement of precipitation is noticed. Figure 8a shows the simulated 6h non-convective rainfall for the FN run whereas 8b shows the non-convective rainfall differences of the FN against the CTRL runs for the 21 st June 2013 at 18.00 UTC. The suppression -enhancement precipitation mechanism is similar to the mechanism of the convective rainfall case. The non-convective rainfall is decreased by up to 8 mm over the area of Area 1 (Western France and NE Iberian Peninsula), where higher concentrations of CCN are simulated reaching values between 450 and 1050 cm −3 .On the contrary, an enhancement in the precipitation up to 10 mm is observed over Area 2 (NE Germany and Central Poland), where CCN concentrations are significantly lower (50 cm −3 or less). The impact of the FN scheme on the precipitation some days after the dust transport event, on 24st June 2013 at 06.00 UTC, has been investigated through the FN simulated 6h non-convective rainfall and the FN-CTRL rainfall differences (Figure 8c and 8d respectively). Once more Area 2 is characterized by CCN concentrations ranging between 50 and 250 cm −3 and indicates an enhancement in the non-convective rainfall To investigate the impact of the FN scheme to the non-convective rainfall, two areas in the D02 are selected, namely Areas 1,2 (Figure 7), as in these areas, significant suppression or enhancement of precipitation is noticed. Figure 8a shows the simulated 6h non-convective rainfall for the FN run whereas 8b shows the non-convective rainfall differences of the FN against the CTRL runs for the 21st June 2013 at 18.00 UTC. The suppression -enhancement precipitation mechanism is similar to the mechanism of the convective rainfall case. The non-convective rainfall is decreased by up to 8 mm over the area of Area 1 (Western France and NE Iberian Peninsula), where higher concentrations of CCN are simulated reaching values between 450 and 1050 cm −3 . On the contrary, an enhancement in the precipitation up to 10 mm is observed over Area 2 (NE Germany and Central Poland), where CCN concentrations are significantly lower (50 cm −3 or less). The impact of the FN scheme on the precipitation some days after the dust transport event, on 24 June 2013 at 06.00 UTC, has been investigated through the FN simulated 6 h non-convective rainfall and the FN-CTRL rainfall differences (Figures 8c and 8d respectively). Once more Area 2 is characterized by CCN concentrations ranging between 50 and 250 cm −3 and indicates an enhancement in the non-convective rainfall reaching 8 mm.     The impact of the FN scheme to the non-convective precipitation pattern is furthermore investigated through vertical cross-sections of FN-CTRL differences of the cloud and rain water mixing ratios (mg kg −1 ). The cross-sections are in longitudinal orientation between points AB for Area 1 and points CD and EF for Area 2 (Figure 7). Figure 9a indicates that, for the AB area, for 21 June 2013 at 18.00 UTC, the FN enhancement to cloud water mixing ratio (shaded) reached values up to 400 mg kg −1 . The increased cloud water mixing ratio is due to the higher CCN levels in the nearby areas, as it is shown in contours in Figure 9b (at least 450-500 cm −3 ). As is mentioned above, significant CCN concentrations in the atmosphere favor evaporation processes and limit the precipitation as the autoconversion of cloud to rain water becomes relatively poor [33,104,105]. Moreover, the high relative humidity (contour lines over 90%) also creates favorable conditions for the cloud water mixing ratio increase (Figure 9a). Due to the poor cloud to rain water autoconversion and the high cloud water mixing ratio, limited amounts of rain water are shown in Figure 9b as the rain water mixing ratio is decreased up to 400 mg kg −1 , leading to the overall suppression of the precipitation for the area, as is also shown in Figure 8b. In the cross-section of Figure 9b, no decrease in the rain water in the lee side of the Pyrenees is observed, as the cross-section orientation does not include the area characterized by significant rainfall suppression (Figure 8b). The above procedure confirms the fact that cloud and rain water production are two inverse processes [106,107], highly dependent on the impact of the CCN to the atmospheric conditions. Figure 9c shows the inverse procedure over Area CD with increased precipitation. In Figure 9c, the cloud water mixing ratio is decreased up to 500 mg kg −1 for 21 June, 2013 at 18.00 UTC. The cloud water mixing ratio is decreased as the specific area is not characterized by low amounts of CCN concentrations (approx. 50 cm −3 , Figure 9d). This leads to high rates of cloud to rain water autoconversion [40,107] as fewer CCN antagonize the same amount of water. Furthermore, the lower relative humidity (50-70% approx.) also causes a weak effect of CCN. On the other hand, in Figure 9d the rain water mixing ratio is enhanced up to 500 mg kg −1 in the layer 2-3 km due to the accretion of cloud water [54].
Similar conditions occur in the EF area in Figure 9e,f for 24 June 2013 at 06.00 UTC, where the cloud water mixing ratio is suppressed, whereas the rain water mixing ratio and the overall precipitation are enhanced, under a low-relative-humidity environment.
In order to assess the impact of the FN scheme on the predictability of precipitation, an overall evaluation for both simulations is performed. The FN and CTRL bias scores, false alarm ratio (FAR) and frequencies of missed forecasts (FOM) of the simulated 24 h precipitation against the WORLDVIEW-NASA satellite retrievals (Atmospheric Infrared Sounder-AIRS-Level 2) are calculated for the Area 1 (AB-90 records) for 21 June, 2013 and Area 2 (CD-60 records, EF-52 records), for 21 and 24 June, 2013 (Tables 2-4, respectively). Nearest grid point value to the satellite ones are chosen. Table 3. FN and CTRL bias, FAR and FOM scores of the simulated 24 h precipitation against the AIRS (WORLDVIEW-NASA) satellite observations for Area AB, as is mentioned in Figure 7. The bold numbers indicate the certain thresholds in which the below scores have been estimated.    Table 4. FN and CTRL bias, FAR and FOM scores of the simulated 24 h precipitation against the AIRS (WORLDVIEW-NASA) satellite observations for Area CD, as is mentioned in Figure 7. The bold numbers indicate the certain thresholds in which the below scores have been estimated. The FN configuration improves generally the bias scores for the entire precipitation thresholds. In more details, over the area AB (Table 3) the systematic overprediction of CTRL for all precipitation thresholds is improved by the FN simulation as the bias scores are significantly reduced. FN also reduces the FAR in comparison to CTRL, which significantly overestimates the precipitation. FOM, however, does not appear any significant difference between the two schemes. As far as the areas CD (Table 4) and EF (Table 5) are concerned, the systematic underestimation in precipitation in the CTRL scheme is replaced by a slighter overestimation in the FN run. In general, the precipitation patterns, especially for heavy precipitation (e.g., at the 25 mm threshold), seem to be massively improved by the implementation of the FN scheme, which leads to an overall predictability improvement. However, for both areas, the FAR scores of the CTRL are better, especially for higher rainfall heights. On the contrary, the FN simulation significantly improves the FOM scores. Table 5. FN and CTRL bias, FAR and FOM scores of the simulated 24 h precipitation against the AIRS (WORLDVIEW-NASA) satellite observations for Area EF, as is mentioned in Figure 7. The bold numbers indicate the certain thresholds in which the below scores have been estimated.

Evaluation of the FN Scheme Through the Optical Properties.
A comparison is made between the estimated aerosol extinction coefficients of both FN and CTRL runs against the EARLINET's network observations (Figure 10a,b, respectively). The stations of Barcelona (BRC) and Potenza (POT) have been selected (Figure 7) using a sample of 96 and 82 records, respectively. Figure 10a shows the comparison for BRC on the 25 June 2013 at 17.07 UTC. Although this date is away from the dust event, it was chosen not only due to the absence of any data set available on the days before (21-23) June, but also in order to demonstrate the impact of the FN implementation to the overall aerosol mass simulation a few days after the main event. An outlier in lidar measurements at 3500m of 0.9/km has been excluded from the statistical procedure. The comparison shows a clear overestimation of the CTRL aerosol extinction coefficient values, especially at the height of 1-5 km. This overestimation is significantly reduced by the FN scheme, as is shown in Table 6. The CTRL overestimation of the extinction coefficients could be attributed to the dust concentration overprediction computed by the GOCART/AFWA scheme [85,89,108]. The CTRL bias of 0.07 reduces to 0.013 for the FN run. The root mean square error (RMSE) is also improved as the CTRL RMSE value of 0.17 is reduced to 0.14 for the FN simulation. set available on the days before (21)(22)(23) June, but also in order to demonstrate the impact of the FN implementation to the overall aerosol mass simulation a few days after the main event. An outlier in lidar measurements at 3500m of 0.9/km has been excluded from the statistical procedure. The comparison shows a clear overestimation of the CTRL aerosol extinction coefficient values, especially at the height of 1-5 km. This overestimation is significantly reduced by the FN scheme, as is shown in Table 6. The CTRL overestimation of the extinction coefficients could be attributed to the dust concentration overprediction computed by the GOCART/AFWA scheme [85,89,108]. The CTRL bias of 0.07 reduces to 0.013 for the FN run. The root mean square error (RMSE) is also improved as the CTRL RMSE value of 0.17 is reduced to 0.14 for the FN simulation.  (Figure 10b) indicates similar results. Once more the CTRL overestimation is reduced by the FN simulation. However, the differences here are greater. The outlier value of 0.14/km near the level of 4500 m has also been excluded from the comparison. Table 6 and Table 7 show the bias and the RMSE differences for BRC and POT, respectively. For BRC, the CTRL overestimation is improved by the FN simulation, as the CTRL bias of 0.070 is clearly reduced to the FN bias value of 0.13. However, there is a systematic overestimation from both simulations over both sites. This may be attributed to the fact that AFWA generally overestimates the dust aerosol concentration as it takes into account higher amounts of clay, at the dust sources parameterizations [85,108]. Furthermore, the observation and the simulation data are derived at the wavelength of 530 nm. This presents a degree of uncertainty because the correlation between the aerosol extinction coefficient and the CCN concentrations could be more satisfactory in lower wavelengths [109] . The CTRL RMSE value of 0.17 is also reduced to 0.14 in the FN scheme. Similarly, for POT, the CTRL overestimation is significantly improved by the FN simulation, as the CTRL bias value of 0.041 is reduced to 0.028. Similarly, the CTRL RMSE value of 0.066 is also reduced to 0.048. However, in the level of 4-4.5 km, an observation value presents a sudden drop reaching 0.15 km −1 . This value cannot be considered as reliable and is not taken into account.  The comparison for the Potenza station on 21 June 2013 at 22.47 UTC (Figure 10b) indicates similar results. Once more the CTRL overestimation is reduced by the FN simulation. However, the differences here are greater. The outlier value of 0.14/km near the level of 4500 m has also been excluded from the comparison. Tables 6 and 7 show the bias and the RMSE differences for BRC and POT, respectively. For BRC, the CTRL overestimation is improved by the FN simulation, as the CTRL bias of 0.070 is clearly reduced to the FN bias value of 0.13. However, there is a systematic overestimation from both simulations over both sites. This may be attributed to the fact that AFWA generally overestimates the dust aerosol concentration as it takes into account higher amounts of clay, at the dust sources parameterizations [85,108]. Furthermore, the observation and the simulation data are derived at the wavelength of 530 nm. This presents a degree of uncertainty because the correlation between the aerosol extinction coefficient and the CCN concentrations could be more satisfactory in lower wavelengths [109]. The CTRL RMSE value of 0.17 is also reduced to 0.14 in the FN scheme. Similarly, for POT, the CTRL overestimation is significantly improved by the FN simulation, as the CTRL bias value of 0.041 is reduced to 0.028. Similarly, the CTRL RMSE value of 0.066 is also reduced to 0.048. However, in the level of 4-4.5 km, an observation value presents a sudden drop reaching 0.15 km −1 . This value cannot be considered as reliable and is not taken into account.

Estimation of CCN vVertical Profiles and Comparison with Observational Data
The performance of the implemented FN scheme is also evaluated through a comparison between vertical profiles of CCN concentrations (cm −3 ) measured by aircraft (SAFIRE/MNPCA-CNRM/GAME) and those simulated by the FN (red line) and CTRL (green line) runs over eastern Sardinia (Domain03) on 22 June 2013 at 10 UTC. Figure 11 shows the locations of the aircraft measurements, while the mean geographic location of the aircraft for each sampling period is considered, in a sample of 180 records. The simulated CCN concentration for each corresponding location is calculated as the weighted average from the nearest grid point of the model. FN simulation estimates considerably higher concentrations than the CTRL (Figure 12). The FN configuration reproduces in a similar way to the observed CCN profile, with an overestimation of approximately 400-500 cm −3 mainly at heights above the atmospheric boundary layer (1000-3000 m). The excess of the estimated CCN concentration by FN simulation is directly linked to the systematic overestimation of the dust concentration by the GOCART-AFWA scheme [89,108]. Table 8 shows the bias and the RMSE differences between the FN and the CTRL runs. The underestimation of the CTRL simulation by about 280 cm −3 turns to overestimation of FN simulation by 38 cm −3 . The implementation of the FN scheme also massively improves the error of the model with 125 cm −3 instead of 358 cm −3 . Additionally, FN also improves the correlation coefficient against the aircraft observations reaching 0.93, as is shown in Figure 13 (green color-solid line).  The performance of the implemented FN scheme is also evaluated through a comparison between vertical profiles of CCN concentrations (cm −3 ) measured by aircraft (SAFIRE/MNPCA-CNRM/GAME) and those simulated by the FN (red line) and CTRL (green line) runs over eastern Sardinia (Domain 03) on 22 June 2013 at 10 UTC. Figure 11 shows the locations of the aircraft measurements, while the mean geographic location of the aircraft for each sampling period is considered, in a sample of 180 records. The simulated CCN concentration for each corresponding location is calculated as the weighted average from the nearest grid point of the model. FN simulation estimates considerably higher concentrations than the CTRL (Figure 12). The FN configuration reproduces in a similar way to the observed CCN profile, with an overestimation of approximately 400-500 cm −3 mainly at heights above the atmospheric boundary layer (1000-3000 m). The excess of the estimated CCN concentration by FN simulation is directly linked to the systematic overestimation of the dust concentration by the GOCART-AFWA scheme [89,108]. Table 8 shows the bias and the RMSE differences between the FN and the CTRL runs. The underestimation of the CTRL simulation by about 280 cm −3 turns to overestimation of FN simulation by 38 cm −3 . The implementation of the FN scheme also massively improves the error of the model with 125 cm −3 instead of 358 cm −3 . Additionally, FN also improves the correlation coefficient against the aircraft observations reaching 0.93, as is shown in Figure 13 (green color-solid line).

Discussion and Conclusions
The potential indirect effects of desert dust aerosols on cloud formation and precipitation have been investigated in this study through the implementation of the revised Fountoukis and Nenes [55] cloud nucleation scheme (FN scheme) in the fully coupled meteorology-chemistry model WRF/Chem. The impact of the above nucleation scheme on both the WDM6 [13,66] double-moment microphysics scheme and the Grell and Freitas [51] cumulus scheme has been taken into account. Dust CCN dynamical impacts were assessed in an extreme Saharan dust outbreak over the Western

Discussion and Conclusions
The potential indirect effects of desert dust aerosols on cloud formation and precipitation have been investigated in this study through the implementation of the revised Fountoukis and Nenes [55] cloud nucleation scheme (FN scheme) in the fully coupled meteorology-chemistry model WRF/Chem. The impact of the above nucleation scheme on both the WDM6 [13,66] double-moment microphysics scheme and the Grell and Freitas [51] cumulus scheme has been taken into account. Dust CCN dynamical impacts were assessed in an extreme Saharan dust outbreak over the Western Mediterranean and Europe accompanied by torrential rainfall on 20-25 June, 2013. Two experiments having as a basis the GOCART-AFWA module with the FN scheme or with the default configuration have been conducted.
As far as convective clouds are concerned, the treatment of CCN based on prognostic dust particles suppresses precipitation over hazy areas, characterized by significant evaporation and limited cloud to rain droplet conversion. More specifically, clouds in such environments are characterized by significant evaporation rates before the precipitation has occurred [33]. Furthermore, in spite of a CCN increase, the same amount of water is distributed among the particles. Thus, the concentration of cloud droplets is increased but their diameters are reduced. In this case, the transformation time of cloud to rain drops increases, especially during the early stages of the cloud formation. Hence, the autoconversion mechanism of cloud to rain droplets and also the precipitation becomes very limited [40,105]. Furthermore, cloud reflectivity is enhanced as cloud surface increases. On the other hand, convective rainfall in a clean atmosphere is enhanced as the exact inverse procedure is occurred. In such cases, the cloud rain out acts very quickly as evaporation is limited.
The non-convective precipitation mechanisms follow a similar enhancement-suppression pattern as above. Hazy environments limit precipitation as high CCN levels lead to significant amounts of cloud water and poor conversion rates of cloud to rain water mixing ratios [33,104,105]. Apart from the enhanced evaporation, high amounts of relative humidity contribute to the above mechanism. On the contrary, pristine conditions lead to the exact opposite results, as limited evaporation conditions, lower relative humidity amounts and decreased cloud water mixing ratios lead to important cloud to rain water conversion and enhancement of precipitation [106,107]. It is remarkable that the prognostic treatment of dust particles enhances significantly the precipitation over areas far away from the dust event. Enhancement of the precipitation over those areas has also been noticed some days after the event.
Comparison of FN and CTRL results of precipitation, extinction coefficients and CCN levels against observations indicated that the FN scheme generally improves the forecasting skill of the model. Bias and RMSE scores are limited in overall as the FN simulation is better correlated with the observed data. CTRL run overestimates precipitation over hazy areas and underestimates it over pristine areas showing important biases especially for thresholds of large precipitation heights. Those biases are improved by the implemented FN scheme, which suppresses the precipitation over the polluted areas and enhances it in clear atmospheric conditions. Concerning hazy areas, the FN simulation improves the FAR scores but does not significantly affect the FOM scores. In pristine areas, however, CTRL has a better performance concerning the FAR. FOM on the other hand improves significantly by the FN simulation. The above result indicates that prognostic CCN treatment is important for both hazy and pristine areas, as it alters the overall precipitation pattern. Furthermore, as far as the extinction coefficient is concerned, the CTRL overestimation is limited by the FN simulation. Moreover, the CCN underprediction on behalf of the CTRL simulation turns to a slight overprediction by the FN run.
The systematic overprediction of both the extinction coefficient and the CCN concentration may be attributed to the overestimation of the GOCART/AFWA scheme in dust aerosol levels. However, this overprediction, as well as the overestimation in precipitation, may also be attributed to the fact that the FN simulation takes into account only prognostic dust particles, since prognostic sea salt particles or other aerosols are neglected in this study. The impact of sea salt aerosols on the model performance is a plan of further investigation.
Author Contributions: K.T. contributed to the methodology and formal analysis, investigation, visualization and writing original draft; P.K. contributed to the conceptualization, methodology, supervision, investigation, writing, review and editing; A.P.; N.M. contributed to the investigation, writing, review and editing. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.