## 2 Materials and Methods
### The FORTH RTM  
The FORTH RTM [51] has evolved from a radiative convective model originally developed by [52]. It has been applied in global [5; 51; 65; 66], regional [55; 67; 68], and local [53] studies. The FORTH's RTM performance has been extensively evaluated against surface-based and satellite reference radiation fluxes [51; 65; 69; 70], as well as against other SW RTMs [57]. A brief evaluation of the RTM solar fluxes at the Earth's surface and at the TOA against corresponding reference data from 7 stations (Figure S2) of the BSRN ground-based network (left panels of Figure S5) and gridded CERES-SYN1 deg satellite measurements over the entire MB and over the 7 BSRN stations (right panels of Figures S5 and S6, respectively), during the three days of the study, the DAEC is provided in Figure S5. The comparison is successful since the 3-hourly matched pairs yielded overall correlation coefficients equal to 0.88 and 0.98, with relative biases equal to 0.8 and 3.5%, respectively, while the RTM TOA fluxes closely follow the intra-daily and the day-to-day variation in the corresponding CERES fluxes.  
The FORTH RTM works in the 0.2 to 10.0 \(\upmu\)m spectral interval and uses the extraterrestrial solar spectrum of [71] weighted to yield a solar constant equal to 1367 W/m\({}^{2}\). The model solves the monochromatic equation of radiation transfer for a multiple scattering/absorbing atmosphere for 118 individual wavelengths between 0.2 and 1.0 \(\upmu\)m and for eight bands between 1.0 and 10.0 \(\upmu\)m, following the delta-Eddington method [72]. It accounts for Rayleigh scattering, O\({}_{3}\) absorption in the ultraviolet and visible wavelength areas, and H\({}_{2}\)O, CH\({}_{4}\), CO\({}_{2}\), absorptions in the near-infrared, as well as aerosol and cloud scattering and absorption and surface reflection. The FORTH RTM version, which is applied here, ran on a 3-hourly temporal and 0.5\({}^{\circ}\)\(\times\) 0.625\({}^{\circ}\) spatial resolution using 3-D aerosol optical properties from MERRA-2 and cloud properties (i.e., cloud amount, optical depth, and top pressure) from ISCCP-H.  
The RTM ran one time with all aerosol types and one with all except for DAs and outputted upwelling and downwelling solar fluxes (F\({}_{\rm all}\) and F\({}_{\rm no-DA}\)) at the TOA, at the surface, and at 50 levels in the atmosphere. Subsequently, the instantaneous DDREs were calculated from the difference between the fluxes computed by taking and not taking into account DAs, according to:  
\[{\rm DDRE_{i}}={\rm F_{i,all}}-{\rm F_{i,no-DA}} \tag{1}\]  
where i stands for the TOA, surface, and atmosphere (specific atmospheric level). Based on the above-mentioned bias of model fluxes, the expected uncertainty of DDREs should be less than 7%. In addition, the 3-hourly heating rates (K/3 h) caused by the evolving dust episode and subsequent absorption of solar radiation are also estimated in order to yield the effect of the dust DREs on the thermal structure and the dynamics of the Mediterranean atmosphere. The computation is made according to:  
\[\frac{\partial{\rm T}}{\partial{\rm t}}=\frac{{\rm g}}{{\rm c_{p}}}\frac{ \Delta{\rm F}}{\Delta{\rm P}}, \tag{2}\]  
where \(\partial{\rm T}/\partial{\rm t}\) is the heating rate, g is the gravitational acceleration, \({\rm c_{p}}\) is the specific heat capacity of air at constant pressure (1006 J/KgK), \(\Delta{\rm F}\) is the \({\rm DDRE_{atm}}\), and P is the atmospheric pressure.  
Lastly, the buoyancy (vertical acceleration) induced by the DDRE was also calculated using the following equation (derived from [PERSON]'s first law and from the ideal gas law):\[\alpha=\frac{(\Gamma^{\prime}-\mathrm{T})}{\mathrm{T}}\cdot\mathrm{g} \tag{3}\]  
where \(\alpha\) is the caused acceleration in m/s\({}^{2}\); T\({}^{\prime}\) and T are the atmospheric temperature with and without the DDRE, respectively; and g is the gravitational acceleration.