Saharan dust intrusions over the northern Mediterranean region in the frame of EARLINET (2014-2017): Properties and impact in radiative forcing

Remote sensing measurements of aerosols using depolarization Raman Lidar systems from 4 EARLINET (European Aerosol research Lidar Network) stations are used for a comprehensive analysis of Saharan dust events 25 over the Mediterranean basin in the period 2014-2017. In this period, we selected to study 51 dust events regarding the geometrical, optical and microphysical properties of dust particles, classifying them and assessing their radiative forcing effect on the atmosphere. From West to East, the stations of Granada, Potenza, Athens and Limassol were selected as representative Mediterranean cities regularly affected by Saharan dust intrusions. Emphasis was given on lidar measurements in the visible (532 nm) and specifically on the consistency of the particle linear 30 depolarization ratio (δp532), the extinction-to-backscatter lidar ratio (LR532) and the Aerosol Optical Thickness (AOT532) within the observed dust layers. We found mean δp532 values of 0.24±0.05, 0.26±0.06, 0.28±0.05 and 0.28±0.04, mean LR532 values of 52±8 sr, 51±9 sr, 52±9 sr and 49±6 sr, and mean AOT532 values of 0.40±0.31, 0.11±0.07, 0.12±0.10 and 0.32±0.17, for Granada, Potenza, Athens and Limassol, respectively. The mean layer thickness values were found to range from ~1700 to ~3400 m. Additionally, based also on a previous aerosol type 35 classification scheme provided by airborne High Spectral Resolution Lidar (HSRL) observations and on air mass backward trajectory analysis, a clustering analysis was performed in order to identify the major mixing aerosol types over the studied area. Furthermore, a synergy of lidar measurements and modeling was used to deeply analyze the solar and thermal radiative forcing of airborne dust. In total, a cooling behavior in the solar range and a significantly lower heating behavior in the thermal range was estimated. Depending on the dust optical and 40 geometrical properties, the load intensity and the solar zenith angle (SZA), the estimated solar radiative forcing values range from -59 to -22 W m at the surface and from -24 to -1 W m at the top of the atmosphere (TOA). Similarly, in the thermal spectral range these values range from +2 to +4 W m for the surface and from +1 to +3 W m for the TOA. Finally, the radiative forcing seems to be inversely proportional to the dust mixing ratio, since higher absolute values are estimated for less mixed dust layers. 45 https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c © Author(s) 2020. CC BY 4.0 License.


Abstract
Remote sensing measurements of aerosols using depolarization Raman Lidar systems from 4 EARLINET (European Aerosol research Lidar Network) stations are used for a comprehensive analysis of Saharan dust events 25 over the Mediterranean basin in the period 2014-2017. In this period, we selected to study 51 dust events regarding the geometrical, optical and microphysical properties of dust particles, classifying them and assessing their radiative forcing effect on the atmosphere. From West to East, the stations of Granada, Potenza, Athens and Limassol were selected as representative Mediterranean cities regularly affected by Saharan dust intrusions. Emphasis was given on lidar measurements in the visible (532 nm) and specifically on the consistency of the particle linear 30 depolarization ratio (δ p532 ), the extinction-to-backscatter lidar ratio (LR 532 ) and the Aerosol Optical Thickness (AOT 532 ) within the observed dust layers. We found mean 532 values of 0.24±0.05, 0.26±0.06, 0.28±0.05 and 0.28±0.04, mean 532 values of 52±8 sr, 51±9 sr, 52±9 sr and 49±6 sr, and mean AOT 532 values of 0.40±0.31, 0.11±0.07, 0.12±0.10 and 0.32±0.17, for Granada, Potenza, Athens and Limassol, respectively. The mean layer thickness values were found to range from ~1700 to ~3400 m. Additionally, based also on a previous aerosol type 35 classification scheme provided by airborne High Spectral Resolution Lidar (HSRL) observations and on air mass backward trajectory analysis, a clustering analysis was performed in order to identify the major mixing aerosol types over the studied area. Furthermore, a synergy of lidar measurements and modeling was used to deeply analyze the solar and thermal radiative forcing of airborne dust. In total, a cooling behavior in the solar range and a significantly lower heating behavior in the thermal range was estimated. Depending on the dust optical and 40 geometrical properties, the load intensity and the solar zenith angle (SZA), the estimated solar radiative forcing values range from -59 to -22 W m -2 at the surface and from -24 to -1 W m -2 at the top of the atmosphere (TOA). Similarly, in the thermal spectral range these values range from +2 to +4 W m -2 for the surface and from +1 to +3 W m -2 for the TOA. Finally, the radiative forcing seems to be inversely proportional to the dust mixing ratio, since higher absolute values are estimated for less mixed dust layers. 45

Introduction
The Saharan desert is one of the major dust sources globally, with dust advections to the Mediterranean countries to be modulated by meteorology along rather regular seasonal patterns (Mona et al., 2012). For instance, in the Western Mediterranean region the African dust occurrence is higher in summer (Salvador et al., 2014), while in the central Mediterranean region, spring and summer are, usually, associated with dust aerosol loads extending up to 50 altitudes of 3-4 km ( Barnaba and Gobbi, 2004). In the Eastern Mediterranean, the main dust transport occurs from spring to autumn (Papayannis et al., 2009;Nisantzi et al., 2015;Soupiona et al., 2018) as a result of the high cyclonic activity over the northern Africa during these periods (Flaounas et al., 2015). Considering also that the Mediterranean basin is a region of high evaporation, low precipitation and remarkable solar activity, the transportation of aerosols accompanied with aging and mixing processes make this area a study of interest for 55 present and future climate change effects (Michaelides et al., 2018).
It is well documented that mineral dust highly influences the atmospheric radiative balance through scattering and absorption processes (direct effect) as well as cloud nucleation, formation and lifetime (indirect effects) as shown in IPCC (2014). Considerable uncertainties in quantifying the global direct radiative effects of aerosols arise from the variability of their spatio-temporal distribution and the aging/mixing processes that can affect their optical, 60 chemical and microphysical properties. Therefore, the magnitude and even the sign of the dust aerosol solar radiative forcing are highly uncertain as they strongly depend on their optical properties, their size distribution and their complex refractive index (CRI) values. Papadimas et al. (2012) reported that the aerosol optical depth seems to be the main parameter for modifying the regional aerosol radiative effects (under cloud-free conditions) and, that on an annual basis, aerosols can induce a significant "planetary" cooling over the broader Mediterranean basin. 65 Other studies (Quijano et al., 2000;Tegen et al., 2010) have shown that the presence of clouds and the surface albedo are also unquestionable parameters affecting the net solar radiative transfer at the top of the atmosphere. However, a comprehensive analysis from ground-based aerosol optical properties to vertical profiles of short-and long-wave radiation estimations in the Mediterranean region has been reported so far only in a few papers (Sicard et al., 2014;Meloni et al., 2003;Gkikas et al. 2018). 70 Although there have been a lot of studies about Saharan dust optical properties based on the lidar technique (Papayannis et al., 2009;Mona et al., 2012;Navas-Guzmán et al., 2013;Soupiona et al., 2018), systematic longterm statistical studies are quite scarce since few aerosol depolarization data are available covering long periods. In Pérez et al. (2006) it is proposed that a regional atmospheric dust model, with integrated dust and atmospheric radiation modules, represents a promising approach for further improvements in numerical weather prediction 75 practice and radiative impact assessment over dust-affected areas, especially in the Mediterranean. Hence, an indepth study of the role of aerosols on the radiative forcing into different regions in the Mediterranean basin is still needed. While a synergy of ground-based lidar measurements and modelling seems very promising for obtaining radiative forcing estimations of dust aerosols, the use of inputs from regional models could also contribute for such estimations in areas where measurements are unavailable. 80 This paper aims to fill some of the aforementioned gaps by combining statistical lidar data of aerosol optical and microphysical properties with radiative transfer estimations and is organized as follows. A brief summary of the four selected EARLINET Mediterranean stations is given in Sect. 2 along with the data selection of the dust cases. Section 3 includes the description of the methodologies applied and the tools and models used for retrieving the aerosol optical and microphysical properties and their radiative forcing. The evaluation of the retrieved aerosol mass 85 concentration profiles and of the ground level radiation are also presented. The results of the aerosol optical, geometrical and microphysical properties of the individual dust layers and the clusters, as well as the relevant radiative forcing calculations over the studied areas are discussed in Sect.4. Finally, concluding remarks are given in Sect. 5.
The European Aerosol Research Lidar Network (EARLINET, https://www.earlinet.org/) established in 2000, provides an excellent opportunity to offer a large collection of quality assured ground-based data of the vertical distribution of the aerosol optical properties over Europe. These measurements meet absolute accuracy standards to achieve the desired confidence in aerosol radiative forcing calculations. Currently, the network includes 31 active lidar stations distributed over Europe, providing information of aerosol vertical distributions on a continental scale. 95 In this paper, level 2 data of four stations from the EARLINET database including aerosol backscatter (βaer) extinction (αaer) coefficient and depolarization ratio (δaer) profiles as a function of height above mean sea level (a.s.l.) were collected and further analyzed as described below to finally estimate their role in radiative transfer calculations in the Mediterranean region.

EARLINET stations 100
Four EARLINET stations affected by typical Saharan dust intrusions in the Mediterranean were selected (listed from West to East): Granada (Spain), Potenza (Italy), Athens (Greece) and Limassol (Cyprus) during a four year (2014-2017) common period of aerosol depolarization Raman lidar data obtained at 532 nm. Table 1 summarizes the basic information about these lidar systems for each location. Except the Limassol station that provides data only at 532 nm, the other three stations are equipped with a depolarization Raman-lidar system able 105 to provide extensive aerosol properties, namely three βaer (355, 532, 1064 nm) and two αaer (355, 532 nm) as well as aerosol intensive properties namely the backscatter and extinction related Ångström exponents (AE α355/532 , AE β355/532 , AE β532/1064 nm), the lidar ratio (LR), and additionally the linear volume ( 532 ) and particle depolarization ratio ( 532 ) at 532 nm. By using the Raman technique, as proposed by Ansmann et al., (1992), we can retrieve the αaer and βaer vertical profiles, with systematic uncertainties of ∼5-15% and ∼10-25%, respectively 110 (Ansmann et al., 1992;Mattis et al., 2002). Therefore, the corresponding systematic uncertainty of the retrieved lidar ratio values is of the order of 11-30%, while the mean uncertainty for AEα and AE β is of order 7-21% and 14-35% respectively, estimated by propagation error calculations.

Selection of dust events
The selection of the lidar dust cases was based on the values of the aerosol optical properties 532 and 532 (Groß et al., 2013). Since pure dust layers are rare over the Mediterranean cities due to continental contamination by urban, polluted or even biomass burning (BB) aerosols, a sufficiently lower 532 value with respect to the pure dust values (e.g. Freudenthaler et al., 2009) should be considered to characterize an aerosol layer as a dusty one. 120 Based on previous studies, the respective 532 values for long-range transported mixtures over the Mediterranean area are expected to range between 35-75 sr (Mona et al., 2006;Papayannis et al., 2008;Tesche et al., 2009;Groß et al., 2011;Ansmann et al., 2012;Nisantzi et al., 2015;Soupiona et al., 2018). Consequently, from the total set of https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
Saharan dust events per station listed in the EARLINET database for the period 2014-2017, we considered for further analysis only the data meeting three basic criteria: a) δ p532 ≥ 0.16 in the free troposphere, b) 35 sr ≤ 125 LR 532 ≤ 75 sr in the free troposphere and c) the thickness of the detected layer to be 500 m, at least. The critical height (in meters a.s.l.) in which the first criterion was met, was considered to be the base of the dust layer. This assumption was deemed necessary to be made since usually, the lofted dust layers cannot be distinguished from the top of the Planetary Boundary layer (PBL), while the presence of urban haze and pollution decreases drastically the values down to 0.03-0.10 (Gobbi et al., 2000;Groß et al., 2013). The top of the dust layer was estimated as the 130 height where the signals were similar to the molecular scattering in the free troposphere. For some cases of the Athens station, where depolarization measurements were unavailable, the values of the base and top were calculated from the Raman lidar signals, following the procedure proposed by Mona et al. (2006).
Moreover, we performed a careful investigation of the air mass origin and dust transport path by means of backward trajectory analysis. This analysis was carried out using the HYbrid Single-Particle Lagrangian Integrated 135 Trajectory (HYSPLIT) model (https://ready.arl.noaa.gov/HYSPLIT_traj.php, (Stein et al., 2015) together with the GDAS (Global Data Analysis System) meteorological files (spatial resolution of 1° × 1°, every 3 hours) as data input. The kinematic back-trajectories were calculated using the vertical velocity component given by the meteorological model with a 96-120 hours pathway (4-5 days back). Modis/Terra information (https://firms.modaps.eosdis.nasa.gov/map) was also taken into account for the corresponding hot spots of possible 140 fires and thermal anomalies along the trajectories (https://firms.modaps.eosdis.nasa.gov/map, not shown here).
Thus, we ended up with 51 individual cases in total (15 for Granada, 18 for Potenza, 12 for Athens and 6 for Limassol). For the region of Cyprus, the situation is more complex since Middle East dust outbreaks occur also frequently in addition to the Saharan dust events (Nisantzi et al., 2015;Kokkalis et al., 2018;Solomos et al., 2019). On top of that, dust particles originating from Middle East proved to have different lidar ratio values than the 145 corresponding observations over Saharan desert (Mamouri et al., 2013;Kim et al., 2020). Taking this into account, dust cases over the Limassol station originating from Middle East regions were excluded from our study.
The air mass trajectory analysis based on HYSPLIT for each station reveals the origin of each observed layer ( Fig.1). In the majority of cases, air masses originate from the west and northwest Africa (Morocco, Mauritania, Algeria and Tunisia). At a first glance two occurrences seem to dominate: i) trajectories that travel directly from 150 the source to the observation stations and ii) trajectories that circulate over the Mediterranean or the Atlantic Ocean (for the Granada cases), Europe and Balkans or even Turkey (for the Limassol cases) before reaching the observation stations.

Methodologies, tools and data evaluation
In order to perform simulations for further investigating the behavior of the transported dust aerosols and their impacts we used different methods with a variety of tools and models. In this section we present our efforts 160 for retrieving vertical dust mass concentration profiles, aerosol microphysical properties and radiative forcing results. Evaluation of the models used was also implemented.

Dust mass concentration lidar retrievals
To retrieve the aerosol dust mass concentration profiles, we used the β 532 and the δ p532 coefficients. Furthermore, by assuming that we have two aerosol types (dust and non-dust) inside the calculated 532 values, 165 we separated the backscatter profiles in two components: the first arising from the contribution of the weakly depolarizing particles ( = 0.05 for non-dust particles) and the second from the contribution of strongly depolarizing particles ( = 0.31 for dust particles). Then, the dust-related backscatter coefficient at 532 nm was obtained, following the procedure described by Tesche et al. (2009). The estimation of the height-resolved mass concentration (in kg m −3 ) of dust particles was based on the procedure described by Ansmann et al. (2012), 170 using the following equation, where in our study the coarse-particle mass density (ρd) was assumed equal to 2.6 gm −3 and a mean volume-to-AOT ratio for coarse mode particles, ⁄ was calculated from AERONET measurements (https://aeronet.gsfc.nasa.gov/) for each station during the period 2014-2017. We report these values as follows: 175 0.71 ± 0.37 μm for Potenza, 0.80 ± 0.29 μm for Granada, 0.94 ± 0.50 μm for Athens, and 0.87 ± 0.27 μm for Limassol. The mean values of the whole studied period were calculated, since only few of the studied cases were common in EARLINET and AERONET database. For LRd the mean values of 52 ± 8 sr, 51 ± 9, 52 ± 9 sr and 49 ± 6 sr were used per site, respectively as calculated from our findings.

The Spheroidal Inversion eXperiments (SphInX) software tool 180
The SphInX software provides an automated process to carry out microphysical retrievals from synthetic and real lidar data inputs and further to evaluate statistically the inversion outcomes. It has been developed at the University of Potsdam (Samaras, 2016) within the Initial Training for atmospheric Remote Sensing (ITaRS) project (2012)(2013)(2014)(2015)(2016). The methodology applied here for spheroid-particle approximation is the same as presented in Soupiona et al. (2019). Raman lidar observations were used as inputs for specific heights within the layers and 185 averaged to produce the 6-point dataset of the so-called 3β par + 2α par + 1δ setup. All cases fulfilling this setup were treated in parallel for retrieving their microphysical properties.

Atmospheric dust cycle model (BCS-DREAM8b)
The BSC-DREAM8b model (Basart et al., 2012), operated by the Barcelona Supercomputer Center (BSC-CNS: www.bsc.es) provides operational forecasts since May 2009, is also participating in the northern Africa-Middle East-Europe (NA-ME-E) node of the SDS-WAS program. The BSC-DREAM8b is a regional model designed to simulate and predict the atmospheric cycle of mineral dust aerosols. It is one of the most widely used 195 and evaluated models for dust studies over northern Africa and Europe (cf. Jiménez-Guerrero et al., 2008;Papayannis et al., 2009;Basart et al., 2012;Amiridis et al., 2013;Tsekeri et al., 2017). The present analysis includes vertical profiles of dust mass concentration simulations (0.3º x 0.3º horizontal resolution, 24 vertical levels, corresponding to the studied cases and for time periods close to the measurement times, usually at 18:00 and 00:00 https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.

Radiative forcing simulations
The aerosol effects on solar and terrestrial radiation are usually quantified through the so-called aerosol radiative forcing (ARF). The dust radiative forcing (DRF) defined here as the perturbation in flux in the atmosphere caused by the presence of the dust aerosols in relation to that calculated under clear sky conditions, can be expressed 205 as (Quijano et al., 2000;Sicard et al., 2014;Mishra et al., 2014): where, the net flux, at a level z, is the difference between the downwelling, and upwelling flux, ↓ and ↑, respectively: These fluxes (in W m -2 ) are calculated separately for SW and LW radiation sources and assuming that the amount of the incoming solar radiation at the TOA is equal for both cases with and without the presence of dust aerosols. Therefore, the net DRF, ( ), is expressed as: Based on this definition, the DRF at a given altitude will be positive when the aerosols cause a heating effect and 215 negative when they cause a cooling effect. Finally, the DRF within the atmosphere ( ) can be defined as the net difference between DRF at the top of the atmosphere (TOA) and the bottom of the atmosphere (BOA), denoted here as and , respectively: Our analysis was based on these equations for estimating the radiative forcing by means of direct and diffuse 220 irradiances of an accurate radiative transfer model combining lidar measurements and dust concentration simulations.

The radiative transfer model (LibRadtran)
In this study, the downwelling and upwelling shortwave (280-2500 nm) and longwave (2.5-40 µm) irradiances at TOA and BOA levels have been simulated with the LibRadtran radiative transfer model version 2.0.2. 225 (Emde et al., 2016). This software package contains numerous tools to perform various aspects of atmospheric radiative transfer calculations. In our study, the uvspec program that calculates the radiation field in the Earth's atmosphere was implemented for the disort radiative transfer equation (1-D geometry). Mid-latitude conditions (AFGL Atmospheric Constituent Profiles, 0-120 km) and a typical surface albedo value (0.16) for urban areas (Dhakal, 2002) in the SW range were taken into account based also on visual observations. The OPAC library 4.0 230 (Koepke et al., 2015) was used for desert spheroids (T-matrix calculations) to determine aerosols' radiative properties in the aforementioned wavelength ranges. The non-spherical approximation is given by typical particle size dependent aspect ratio distributions of spheroids, derived from measurements at observation campaigns. In our study, the mineral particles of each case are given to the model as spheroids for the mineral accumulation mode (MIAM) with RMIAM ∈ [0.005, 20] in μm. 235 A set of four simulations was carried out per case of the studied dust events. The first two simulations refer to clear-sky atmospheres with background/baseline aerosol conditions (default properties: rural type aerosol in the boundary layer, background aerosol above 2 km, spring-summer conditions and a visibility of 50 km), for SW and for LW range, treated separately. The other two simulations correspond to dust loaded atmosphere, also in SW and LW ranges, respectively, for which vertical profiles of different components of dust layers were used as additional 240 https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
inputs. These inputs have been obtained with three different Schemes: A) vertical mass concentration profiles simulated by the BCS-DREAM8b model, B) vertical dust mass concentration profiles retrieved based on the 532 and 532 profiles, and, C) vertical profiles of α532 along with the respective mean AOT532 value. The flow chart in Fig. 2 depicts these three Schemes applied to create the input files for the dust-loaded atmospheric conditions used in LibRadtran software package (Emde et al., 2016). For Scheme C conversion factors from OPAC were used in 245 order to convert the αpar and the corresponding AOT from 532 nm to 10 μm (peak, within the atmospheric window). The conversion was based on an adaptive inversion algorithm of Shang et al. (2018) who presented a way to convert extinction coefficients at different wavelengths by using Angstrom exponent values derived from AOTs.
For all these Schemes in this study, 30 vertical levels have been used between ground and 120 km height with a spatial vertical resolution of 0.5 km starting from ground level (BOA) to 2 km and from 5 to 10 km, a 250 resolution of 0.2 km from 2 to 5 km, due to the presence of the dust layers within this height range and additionally at the heights of 20 and 120 km (TOA). All simulations were performed for three different Solar Zenith Angles (SZA), 25 o 45 o and 65 o covering a typical diurnal cycle for radiative forcing estimates at mid-latitudes. All cases were treated for cloud-free conditions. Except the altitude in km (z out ), the additional outputs that have been implemented in our Schemes are: the direct horizontal irradiance (e dir ), the global irradiance (e glo ) the diffuse 255 downward irradiance (e dn ), the diffuse upward irradiance (e up ), and the heating rates (heat) in K day -1 , as described by Mayer et al. (2017).

Radiation data set 260
LibRadtran irradiance outputs have been validated against reference solar irradiance pyranometer measurements at the Earth's surface (Kosmopoulos et al., 2018). For this study, solar radiation data measured by pyranometers were available only for the Granada and Athens stations. The reference solar radiation data set consists of one-minute simultaneous measurements of horizontal global and diffuse irradiance measured with two CMP11 pyranometers, at Granada, and two CMP21, at Athens (NOA's actinometric station in Penteli area). For those diffuse irradiance measurements taken using a shadow ring, the model proposed by Drummond (1956) has been apply in order to correct for the diffuse radiation intercepted by the ring, as suggested the manufacturer (Kipp & Zonen, 2004). 275

Evaluation of aerosol mass concentration vertical profiles
Before using the aerosol mass concentrations vertical profiles from i) BSC-DREAM8b model simulations (Scheme A) and ii) hourly averaged lidar measurements (Scheme B) as inputs to the LibRadtran model, we performed a day-by-day comparison between them. Due to the different spatial vertical resolution between the modeled and the lidar profiles, both profiles were degraded to the fixed height levels of the OPAC dataset (0, 1, 2, 280 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 35 km). Figure 3 shows the Taylor's diagram of the mass concentration simulated by the BSC-DREAM8b model against the lidar-retrieved ones. The azimuthal angle presents the correlation coefficient; the radial distance presents the normalized SD of the simulated patterns, the distance from the point on the x-axis presents the root mean square error (RMSE), while the black point at the (1,1) cross section, indicates the lidar retrieved aerosol mass values 285 representing the reference point. The normalization of the SD is performed with respect to the calculated values. It is evident that the shape of the vertical distribution was better reproduced than the intensity of each event.
Specifically, in the 66 % of the cases there is a good correlation (r > 0.6), in the shape of the vertical distribution, while in 96 % of the cases the model gives lower concentration values (normalized < 1). Therefore, we report a mean underestimation of the mean mass concentration values of the BSC-DREAM8b of the order of 31%. Also, 290 the mean center of mass (in km) estimated from BSC-DREAM8b profiles is 0.6 km lower than the one calculated from the lidar measurements ( 2.6 ± 1.0 km and 3.2 ± 1.1 km respectively). However, we should take into consideration: i) the spatial resolution, where the lidar observations are considered as point measurements while the simulations represent uniform pixels of 0.3 o resolution and ii) the temporal resolution, where the lidar retrieved profiles are hourly averaged, while the model derived profiles are instantaneous results, saved every 6 hours. 295 Therefore, our observations are in line with previous studies of Mona et al. (2014) and Binietoglou et al. (2015) where have reported discrepancies concerning the base, the top layer height and extinction profiles and good agreement in terms of profile shape, between the BSC-DREAM8b and observations.

Evaluation of ground level LibRadtran outputs
The evaluation of the performance of the model was undertaken by statistical means. The relative root mean square error (rRMSE), the relative mean bias error (rMBE) the correlation coefficient (r) and the normalized SD https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
were calculated in order to numerically quantify the performance of the global irradiance calculated from 305 pyranometers and simulated from the three Schemes. Table 2 shows the statistical results for the modeled global irradiance values versus the reference pyranometer measurements for both locations and the threes schemes. All scheme simulations perform remarkably well, with rRMSE values ranging from 8.3 to 16.2% and rMBE values between 0 and 15.2%. In general, the rRMSE is slightly higher at Granada, mainly for the Scheme A. According to this statistic, the LibRadtran output with the best performance are those obtained with the Scheme C as input 310 followed by Scheme B and A, respectively. This order is the same attending to the rMBE values with the exception of the Scheme A at Athens. The correlation coefficient r depicts the good performance of the radiative transfer model for the three schemes and the two locations models. All simulations present a r > 0.95 with minor differences (below a 10%) in the normalized SD values respect to the pyranometer global irradiance values. A slight overestimation is observed for all scheme outputs at Granada (norm SD > 1). Conversely, this overestimation is no 315 longer evident in the modeled global irradiance for Athens. It is important to notice the good performance of the Scheme B despite the important calculus involved in it.

Aerosol geometrical and optical properties per site
For each case studied, the mean δ p532 , mean LR 532 and mean AOT 532 values were calculated inside the dust layers as shown in Fig.4 (a-d). The corresponding standard deviation (SD) values give an indication of the variability of these parameters from base to the top of the dust layer. Figure 4a shows the aerosol geometrical properties for 325 the detected layers one by one, per station and per year. The mean values of the base and top height of the dust layers per station, along with their SD are marked with horizontal bounded lines. At the four sites (Granada, Potenza, Athens and Limassol) mean layer thicknesses of 3392 ± 1458 m, 2150 ± 1082 m, 1872 ± 816 m and 1716 ± 567 m were calculated respectively. We also found mean δ p532 values of 0.24 ± 0.05, 0.26 ± 0.06, 0.28 ± 0.05 and 0.28 ± 0.04 (Fig. 4b), and mean 532 values of 0.40 ± 0.31, 0.11 ± 0.07, 0.12 ± 0.10 and 0.32 ± 0.17 330 (Fig. 4d), respectively. Similar mean 532 values of around 51 sr (Fig. 4c) were found for all stations. The Granada station holds the minimum mean value for layers' base height (1567 ± 788 m) and the maximum for top height (4960 ± 975 m) and layer's thickness. Concerning the LR values no remarkable deviations were observed among the four stations having mean values around 50 sr, which are in very good agreement with literature findings (Tesche et al., 2009;Ansmann et al., 2012;Groß et al., 2011;. The largest mean AOT value, equal to 0.40 ± 0.31, 335 observed over Granada station is in accordance to the geometrical properties of Fig.4Figure a that depicts thick dust layers in the majority of the cases. https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License. Considering Granada's station as representative of the west Mediterranean region, Potenza of the central Mediterranean region, Athens and Limassol stations of the eastern Mediterranean region, a dust aerosol mode classification per region can be made. For this purpose, the mean AOT532 versus the AE β532/1064 giving an indication of the aerosol particle size in the atmospheric column relevant to the aerosol load for each region are 345 shown in Fig. 5. A wide spread of the AOT values at moderate to low AE β532/1064 values (between 0 and 0.6) observed in the Western Mediterranean region, demonstrates that the dust size distribution in this area is dominated by coarse mode particles during events of different intensities. On the other hand, the presence of dust layers, in the central and eastern Mediterranean regions can be associated with higher AE β532/1064 values (even up to 1.5) and consequently, with the presence of fine mode particles and lower dust loads. Our findings verify that the longer the 350 dust transport is, more likely it is for the dust aerosol properties to be modified, since dust can frequently be mixed with background aerosols in the eastern Mediterranean (Groß et al., 2019).
In terms of the aerosol size distribution, the scatter plot of Fig. 5 allowed us to perform a k-means clustering (Arthur and Vassilvitskii, 2007) in order to define three physically interpretable aerosol size distributions: a) fine mode, with AE β532/1064 > 0.6, b) coarse mode, with AE β532/1064 ≤ 0.6 and AOT532 between 0 and 0.2, c) whilst 355 AE β532/1064 values smaller than 0.6 attributed to large AOTs (between 0.2 and 0.8) are representative of extreme dust events. It seems that the majority of these extreme dust outbreaks occur over the Western Mediterranean region, more likely due to its location close to the African continent. with AOTs > 0.2 (675 nm) and values around zero. More studies referring to the occurrence of extreme dust 360 events over the aforementioned area can be found in literature (Cachorro et al., 2008;Guerrero-Rascado et al., 2009).

Clustering per mixing state
Based on the High Spectral Resolution Lidar (HSRL) classification presented by Groß et al. (2013), we plotted the intensive aerosol quantities 532 versus 532 and we identified three of the six existing clusters in our data (Fig 6). The first cluster (green marks and error bars) represents a mixing state of Saharan dust and BB aerosols having a large spreading in mean LR values and low mean 532 values (40 sr ≤ LR 532 ≤ 75 sr, 0.15 ≤ δ p532 ≤ 370 0.19). The second one, (red marks and error bars) is attributed to mixed Saharan dust, where dust aerosols are dominant, but urban/continental, marine or even pollen aerosols are also possibly present (40 sr ≤ LR 532 ≤65 sr, 0.20 ≤ δ p532 ≤ 0.29). The third cluster (orange marks and error bars) is attributed to pure Saharan dust aerosols (45 sr ≤ LR 532 ≤ 60 sr, 0.30 ≤ δ p532 ≤ 0.36). The most populated and consequently, the most common, among those three clusters is the red one, as expected, due to the frequent mixing of dust aerosols with continental ones 375 (Papayannis et al, 2008). The range of our measured δ p532 values as indicated with the horizontal error bars in Fig.  6, overlap between the three identified aerosol clusters, showing a more realistic transition from one cluster to the other, bridging the gap specifically between green and red clusters from the HSRL classification.
https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.  Table 3 summarizes the mean values of the aerosol geometrical, optical, and microphysical properties of the three identified clusters. A synergistic approach of HYSPLIT (trajectories of 120 hours backward for each case) and Google Earth (distance calculator) allowed us to estimate the distance travelled (in km) to the respective sites and the mixing hours per cluster. Specifically, the term of mixing refers to the hours the air masses travelled after 385 leaving the African continent. We can see that the Saharan dust cluster of the air masses that present the lowest mixing with other air masses (26 ± 13 hours), compared to the other clusters (41 ± 26 hours for the BB & mixtures and 31 ± 13 hours for the mixed Saharan dust cluster). Moreover, the air masses of the Saharan dust cluster seem to travel faster than those of the other clusters, although covering a greater distance (4845 ± 2825 km) at the same time (within 120 hours). Now, the main difference between the two remaining clusters (BB & mixtures/Mixed 390 Saharan Dust) is attributed to the mixing hours. The air masses of the latter cluster remain 15 hours longer and circulate over the Mediterranean and Europe, so they are probably enriched with other types of aerosols.
Concerning the aerosol optical properties, the β532 and α532 show lower values for BB & dust and for mixed Saharan dust cases (1.10 ± 0.15 10 −3 km -1 sr -1 , 0.47 ± 0.28 km -1 and 1.24 ± 0.80 10 −3 km -1 sr -1 , 0.74 ± 0.48 km -1 respectively) and higher values (1.54 ± 1.05 10 −3 km -1 sr -1 , 0.80 ± 0.27 km -1 ) for the Saharan dust cluster. 395 Therefore, higher AOT values (0.32 ± 0.25) were found for the latter cluster compared to the others, due to the higher dust burden of these events over the affected sites. The highest 532 values (0.32 ± 0.02), indicate the arid origin and the coarse mode of pure Saharan dust layers (Freudenthaler et al., 2009) of the corresponding cluster. No direct information can be extracted from the similar 532 values about the mixing state of the aerosol layer, except that the range of the SD narrows as the mixing decreases. However, for the cases that observations at 355 nm were 400 available, it seems that the lidar ratio color ratio (namely the LR 355 / LR 532 ) converges to unit for the Saharan Dust cluster, indicating the absence of spectral dependence for the case of pure dust (Müller et al., 2007;Veselovskii et al., 2020). For these cases also, the AE β355/532 becomes closer to zero taking mean value of 0.35 ± 0.45.
We also summarise the changes in mean microphysical properties estimated with SphInX for the three identified clusters. The BB & Saharan dust cluster takes the lower mean Reff value (0.293 ± 0.074 μm) due to the 405 fine structure of BB aerosols included in the layer, while a mean Reff of 0.360 ± 0.081 μm corresponds to the cluster of mixed Saharan dust and a slightly larger value (0.387 ± 0.070 μm) corresponds to the Saharan dust cluster. The values for RRI, IRI and SSA at 532 nm were similar for the two clusters that not include BB aerosols, whilst the presence of BB aerosols of the first cluster leads to higher RRI and IRI values and lower SSA, results https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
that are in good agreement with the ones reported in Petzold et al. (2011) over Dakar, for mineral dust and dust 410 mixed with anthropogenic pollution.

Regional dust solar radiative forcing (DRF) 415
As mentioned previously, there is shortage of papers in the literature about the role of dust on the Earth's radiation budget. Since very few in situ direct measurements of DRF effects and heat fluxes are available especially in the Mediterranean (Bauer et al., 2011;Meloni et al., 2018) we are restricted to perform simulations to quantify the role of dust aerosols on the radiative forcing in the studied regions. The mean DRF is calculated during this simulation, calling twice the LibRadtran radiation code: with and without dust aerosols. For all cases the vertical 420 profile of DRF starting from ground level/bottom of atmosphere (BOA) up to the top of atmosphere (TOA) in the SW and LW ranges were simulated using the three aforementioned Schemes.
A negative forcing of aerosols both at the BOA and TOA is noted in the SW range, as presented in Fig.  7Figure a, which depicts the mean DRF of all cases per scheme. All results within the SW range represent less absorbing aerosols with a cooling behavior. Depending on the dust optical properties and load intensity, DRF values 425 at the BOA range from -40 to -13 W m -2 at SZA 25 o , from -43 to -14 W m -2 at SZA 45 o and from -44 to -15 W m -2 at 65 o . At the TOA, the corresponding ranges per SZA are -9.5 to -1.4 W m -2 (25 o ), -16 to -3.3 W m -2 (45 o ) and -24.3 to -6.9 W m -2 (65 o ). Similarly, in the SZA independent LW range (thermal spectral range), the DRF values range from +1.6 to +4.6 W m -2 for the BOA and from +0.8 to +3.6 W m -2 for the TOA. Our estimations are consistent with results obtained by other literature findings for Saharan dust aerosols over the Mediterranean region. 430 Specifically, Sicard et al. (2014) found that the SW RF at the BOA has always a cooling effect varying from -93.1 to -0.5 W m −2 while the corresponding LW RF has always a heating effect varying from +2.8 to +10.2 W m −2 . They also concluded that dust aerosols at the TOA have a cooling effect in the SW spectral range with a RF ranging from -24.6 to -1.3 W m −2 while at the TOA the LW RF varies between +0.6 and +5.8 W m −2 . Meloni et al. (2003) found at the island of Lampedusa instantaneous RF of -70.8 W m −2 at the BOA and -7.5 W m −2 at the TOA within the 435 range 300-800 nm for an event with AOT of 0.51 at 415 nm. For the same location and for another strong Saharan dust outbreak (AOT500=0.59), Meloni et al. (2015) reported a total (SW + LW) radiative forcing of -48.9 W m -2 at the BOA, -40.5 W m -2 at TOA, and +8.4 W m -2 in the atmosphere for SZA=55.1 o . A negative radiative effect https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.
reaching down to -34.8 W m -2 at the surface in the Mediterranean area was also recently reported by Gkikas et al. (2018) during the forecast period March 2000-February 2013. 440 Variations among these values are expected since they strongly depend on the different aerosol AOTs, mass estimations and extinction values. Estimations retrieved from Scheme B are expected to give higher values compared to those given from Scheme A as revealed also by Fig.3. The DRF at the LW spectral region is opposite in sign and significantly lower in absolute values than in the SW region. The difference between the TOA and BOA DRF, with the former to be only weakly perturbed and the latter to be quite stronger, can be attributed to the heating 445 within the troposphere, since the presence of the dust mainly leads to a displace of surface's radiative heating into the dust layer. Low reflected solar flux is partially offset by the absorption of upwelling LW radiation. Finally, in the LW spectral region, the mean DRF values at the BOA (Scheme A: +1.6±1.6 W m -2 , Scheme B: +4.6±4.7 W m -2 , and Scheme C: +2.9±9.4 W m -2 ) are higher than those at the TOA (Scheme A: +0.8±0.9 W m -2 , Scheme B: +3.6±4.4 W m -2 , and Scheme C: +1.2±6.2 W m -2 ) due to the fact that the main source of LW radiation (Earth's 450 surface) is close to the aerosol layers, mainly observed between 2 and 4 km a.s.l.. As a result, the DRF Atm is positive during the diurnal circle, yielding net radiative heating of the dust layer.
The mean net heating rate within the atmosphere, calculated by adding algebraically both rates in the SW and LW spectral ranges is presented in Fig. 7b. Here, the net heating rate is clearly dependent on the available solar radiation, that increases with SZA due to the low incoming solar radiation reaching the BOA at afternoon hours 455 (SZA 65 o ). Our estimations are in accordance with the fact that as the SZA increases, the optical path of the SW radiation grows significantly, increasing the attenuation of the direct radiation while generating a higher fraction of the diffuse radiation. This effect is more pronounced at the BOA, in which, the intensity of the heating rate is reduced with increasing SZA, since fewer photons are available to heat the dust layers.  ). For Scheme C, we report higher values of this parameter during the diurnal circle. More precisely, the net heating rate is almost 1.5 times higher at 25 o , 2 times higher at 45 o and around 0.8 times higher at 65 o , compared to the aforementioned Schemes. Greater sensitivity in the SZA appears in Scheme B as it results from the line slope. 465 Figure 7: Mean values of a) SW and LW DRF at BOA and TOA and b) the net heating rate within the atmosphere, along with their SD for the three Schemes applied in the total set of the studied cases. The inserted box depicts the line slope.
In order to further explain the difference of sign in the net heating rate of Scheme C, compared to the two others presented in Fig. 7b, we plotted the aforementioned parameter along with the base layer height, the AOT532 470 and the layer thickness of each case as presented in Fig. 8 (Table 3). Focusing on the SW range, the cooling effect for Scheme A of the Saharan dust cluster is up to 3 times 500 higher compared to the BB and Saharan dust one, whilst the cooling effect for Scheme C of the former cluster is up to 10 times higher compared to the latter. The cooling effect of Scheme B becomes also stronger with the decreasing mixing state but in a lower magnitude (the former cluster is almost 2 times higher compared to the latter). Hence, even though the cases included in the Saharan dust cluster usually take higher mass concentration values than the other cases, as predicted by BSC-DREAM8b (Scheme A), still it seemingly underestimates the intensity of strong 505 https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License. transported dust episodes over the observation stations. On the contrary, Scheme C is the most sensitive to the mixing state. To explain this result one should consider that on the one hand, spheroidal particles such as dust have larger surface area than spherical ones such as BB aerosols leading to larger AOTs (Haapanala et al., 2012) and consequently to increased negative DRF and on the other hand, the Schemes A and B involve greater assumptions concerning dust particles than Scheme C. 510 Finally, our interest is focused on the vertical DRF profiles from the surface (a.s.l.) up to 10 km height in the 515 free troposphere, where airborne dust is usually found, as estimated by Scheme C at 45 o SZA per station. The DRF profiles, in the SW region, presented in Figs. 10 (a-d), follow the aerosol extinction vertical structure. The DRF values at the BOA are high in absolute values with a cooling behavior and decreases with increasing height, while the magnitude is proportional to the aerosol load in the whole atmospheric column. Specifically, the DRF ranges from -150.0 to -1.9 W m -2 for Granada, from -38.1 to -3.7 W m -2 for Potenza, from -64.8 to -13.2 W m -2 for Athens 520 and from -90.3 to -28.4 W m -2 for Limassol. The corresponding ranges of α532 are 0.286-0.029 km -1 , 0.268-0.088 km -1 , 0.135-0.078 km -1 and 0.547-0.214 km -1 , respectively. Peaks in α532 are observed usually between 2 and 6 km a.s.l. indicating the intrusion of dust that corresponds to a decrease in the solar radiation reaching the surface.

B O A 2 5 B O A 4 5 B O A 6 5 T O A 2 5 T O A 4 5 T O A 6 5 B O A T O
https://doi.org/10.5194/acp-2020-611 Preprint. Discussion started: 1 July 2020 c Author(s) 2020. CC BY 4.0 License.

Conclusions
The characteristics of the dust aerosol optical, geometrical, and radiative properties over the Mediterranean region, were presented at this study. A total of 51 independent aerosol lidar measurements of Saharan dust events studied over 4 southern European cities were carefully selected and analyzed. The dust layers were usually observed 530 between ~1.6 and ~5 km with 532 and 532 values ranging from 0.16 to 0.35 and from 35 to 73 sr respectively, depending on the air mass mixing state. Significantly high AOT 532 values (0.40 ± 0.31) were found for Granada indicating that the dust outbreaks occurring over this area are more intense. Results of 532 versus 532 are presented in order to elucidate the difference of pure dust and dust mixtures cases. Layers with lower δ p532 (0.17 ± 0.01), AOT 532 (0.03 ± 0.02) and thicknesses (786 ± 212 m) have shown high dust mixing ratio, while the 535 properties of the least or no mixed dust layers (δ p532 =0.32±0.02, AOT 532 =0.32±0.23 and thickness=3158±1605 m) are in a good agreement with literature findings for pure Saharan dust cases (Tesche et al., 2009;Papayannis et al., 2009;Ansmann et al., 2012;Mona et al., 2012;Groß et al., 2011;. Lidar stand-alone microphysical retrievals of Reff, RRI and IRI are also differentiated by the level of mixing. Three Schemes have been implemented to evaluate the DRF during the selected dust outbreaks: the model 540 mass concentrations by BSC-DREAM8b (Scheme A), the mass concentrations calculated from lidar measurements (Scheme B) and the lidar extinction profiles along with the mean AOT532 values (Scheme C). The DRF variations are strong and result from significant changes in the lidar retrieved extensive optical properties (α 532 , β 532 , AOT 532 ) or model mass estimations from the BSC-DREAM8b. A mean underestimation 31% of dust mass by BSC-DREAM8b is leading to a significantly lower SW DRF compared to results obtained from the implementation of 545 lidar profiles in LibRadtran. Additional variations in the SW range are introduced due to the variations in the available solar radiation during day that depends on the SZA. However, variations of the latter are significantly lower compared to the former. The vertical structure of a layer that gives information about the base, the thickness and the intensity (AOT) of a dust event is critically important, while additional information of its mixing state can be also significant in RF and net heating rate estimations. 550 Our findings show a much more pronounced DRF at the BOA compared to the one at the TOA due to the low altitude of the studied layers (usually 2-4 km). Taking this into account we report DRF values at the BOA ranging from -40 to -13 W m -2 at SZA 25 o , from -43 to -14 W m -2 at SZA 45 o and from -44 to -15 W m -2 at 65 o . At the TOA, the corresponding ranges per SZA are -9.5 to -1.4 W m -2 (25 o ), -16 to -3.3 W m -2 (45 o ) and -24.3 to -6.9 W m -2 (65 o ). In the LW range these values range from +1.6 to +3.5 W m -2 for the BOA and from +0.8 to +2.6 W m -555 2 for the TOA where again, the values of the latter reveal less perturbation.
The dependence of DRF on the accuracy of the vertical aerosol mass concentrations can give a wide range of possible values, however, the implemented Schemes can contribute to the characterization of the aerosols' radiative forcing effects over the Mediterranean region, being one of the most sensitive regions to climate change forcing. Our findings suggest that a further investigation of aerosols' mixing state is needed since, not only their 560 optical but also their microphysical properties and radiative forcing can strongly vary, depending on the mixing types. We recommend that the synergistic usage of observations and modeling in the next generation state of the art dust cycle models for the DRF should be intensified. The systematic use of remote sensing profiling measurements to radiative transfer models is stressed in this study, creating an essential tool that will allow us to quantify both the cooling at the surface and the heating within the atmosphere produced by different aerosol types 565 such as dust and its mixtures on a regional and a global scale.
the European Unionʼs Horizon 2020 research and innovation program under grant agreement no. 654169. The authors also acknowledge the BSC-DREAM8b model, operated by the Barcelona Supercomputing Center, the NOAA Air Resources Laboratory (ARL) for the provision of the HYSPLIT, the Google Earth and the AERONET for high-quality sun/sky photometer measurements. A.P., R.F., and C.A.P. acknowledge support by the project "PANhellenic infrastructure for Atmospheric Composition and climatE change" (MIS 5021516) which is 575 implemented under the Action "Reinforcement of the Research and Innovation Infrastructure", funded by the Operational Programme "Competitiveness, Entrepreneurship and Innovation" (NSRF 2014-2020) and co-financed by Greece and the European Union (European Regional Development Fund).