Fine and coarse dust radiative impact during an intense Saharan dust outbreak over the Iberian Peninsula – short-wave direct radiative effect

Mineral dust has a key role in the Earth's radiative balance, and it has become significant over the Iberian Peninsula (IP), where Saharan dust outbreaks have seemed to increase in frequency and intensity. This study quantifies the dust direct radiative effect (DRE) in the short-wave range (SW) during an intense persistent springtime dust episode over the IP. In particular, the vertical distribution of dust optical properties was derived at five lidar stations, and the Global Atmospheric Model (GAME) was used for radiative transfer simulations. Moreover, this study innovates by simulating the SW DRE using two distinct methodologies. The novel approach separates the impacts of fine (Df) and coarse (Dc) dust components and calculates the total DRE as their combined sum. In contrast, the commonly used approach directly simulates the DRE for the total dust. Along the dust pathway crossing the IP, the SW DRE consistently registered a pattern of aerosol-induced cooling at both the surface, i.e the bottom of the atmosphere (BOA), and top of the atmosphere (TOA). Results reveal that the role of Df must be highlighted, as Df particles contributed nearly half of the total SW DRE at BOA and TOA, particularly for this event. DRE simulations based on the separation of dust modes for solar zenith angles < 70° show that nearly 75 % of DRE values are lower (in absolute value) with respect to those obtained by considering the total dust. That is, a less pronounced cooling effect is observed overall when the separation of dust modes is taken into account, although relative differences between approaches are not highly significant in general (−3 % and −6 % on average at BOA and TOA, respectively). However, this behaviour reverses under moderate-to-high-dust conditions when the contributions of the Dc and Df components are considered separately; i.e. dust induces a more pronounced cooling effect. This indicates the relevant role of the fine dust particles in DRE estimations, which is likely underestimated when total dust (no separation) is taken into account in relatively high-dust environments. In addition, the cooling effect is more evident at TOA than at BOA, probably due to the presence of dust at higher levels than those usual in the troposphere. This fact can introduce relevant implications for radiometric measurements from satellite instrumentation.

1 Introduction

According to the latest Intergovernmental Panel on Climate Change technical report (IPCC, 2023), not only is the Earth's global average temperature rising, but an increase in the frequency, intensity, and duration of extreme hot events and extreme precipitation has also been observed over several zones in the last few years. In this context of climate change, aerosols play a crucial role in climate forcing, through both their direct impact by interacting with radiation (aerosol–radiation interactions, ARI) and their indirect effect by serving as cloud condensation nuclei and ice nucleating particles (aerosol–cloud interactions, ACI). From 1750 to 2014, aerosols contributed −1.3 W m⁻² to the global effective radiative forcing at the top of the atmosphere (TOA), resulting in a cooling effect (Arias et al., 2021) that partially compensates for the warming effect by increasing greenhouse gas emissions.

Among aerosol particles, mineral dust is one of the most abundant and widely distributed particles in the Earth's atmosphere. Significant uncertainties remain in the estimation of the dust radiative effect, mostly due to the lack of observational constraints. Regarding the SW range, the balance between scattering and absorption is determined by the dust particle size and mineralogy, and the underlying surface's albedo determines the extent to which both processes impact the TOA radiative flux (Kok et al., 2023). In particular, the Sahara desert and Sahel regions are the major contributors to mineral dust lifting, accounting for over 50 % of global emissions (Kok et al., 2021). On one hand, the geographical proximity of the Iberian Peninsula (IP) to north Africa, as well as the persistence of certain favourable weather patterns (Russo et al., 2020; Couto et al., 2021), makes the IP one of the main pathway regions for Saharan desert dust transport towards Europe. On the other hand, Sánchez-Benítez et al. (2020) found that Iberian heatwave events are primarily linked to anomalous atmospheric circulation, typically characterized by a pronounced positive 500 hPa geopotential height anomaly (Z500) aloft. In particular, an analysis of daily weather regimes during Iberian heatwave days reveals a marked dominance of positive Z500 anomalies over western Europe, resembling the occurrence of Euro-Atlantic subtropical ridges. The frequency of this occurrence during heatwave days has been observed to double relative to climatological averages. Recent studies linked both processes, showing robust evidence for increases in maximum temperatures and the frequency of heatwaves over Europe (IPCC, 2023), which can be partially associated with dust intrusions (Sousa et al., 2019; Fernandes and Fragoso, 2021; Barriopedro et al., 2023). While there are certain indications that future global dust aerosol concentrations may increase (Allen et al., 2016; Tegen and Schepanski, 2018), this likelihood is notably contingent upon alterations in precipitation patterns and atmospheric circulation (IPCC, 2023). Consequently, the radiative consequences of such potential changes remain unclear (Allen et al., 2016; Kok et al., 2018). Since preindustrial times, the dust mass that has lifted into the atmosphere from those north African regions has increased by 46 % (Kok et al., 2023). Moreover, there is an increasing frequency of not only Saharan dust outbreaks over the IP when compared to long-term records (Sousa et al., 2019; Salvador et al., 2022) but also the number of studies reporting extreme and intense dust episodes in the region (e.g. Guerrero-Rascado et al., 2009; Preißler et al., 2011; Santos et al., 2013; Cazorla et al., 2017; Sorribas et al., 2017; Titos et al., 2017; Valenzuela et al., 2017; Córdoba-Jabonero et al., 2019, 2021; Fernández et al., 2019; López-Cayuela et al., 2023).

This study focuses on the interaction of Saharan dust and radiation, i.e. on the direct radiative effect (DRE) of dust particles. Using observational measurements, estimation of the DRE is performed, underscoring the importance of determining the radiative properties of the dust particles and accurately quantifying their direct radiative impact on the Earth–atmosphere radiative budget. This work aims to study the DRE during an exceptionally intense and long-lasting (14 d, from 25 March to 7 April 2021) Saharan dust event over the IP, allowing the study of changes in the radiative properties across both temporal and spatial dimensions. This event was well characterized in the work of López-Cayuela et al. (2023), which determined the vertical behaviour of the optical and microphysical dust properties using data obtained from five Iberian lidar stations. This study focuses on the direct radiative effect of the dust particles in the short-wave (SW) spectral range.

The use of lidar measurements should be emphasized as an added value of this study. Thanks to the application of the POLIPHON (Polarisation Lidar photometer Networking; Mamouri and Ansmann, 2014, 2017; Ansmann et al., 2019) algorithm to the polarized lidar measurements, it is possible to separate the fine and coarse dust contributions (Df and Dc, respectively) to the total DRE. In this work, the Df particle role has been highlighted in contrast to other studies, in which the potential impact of Df particles on the total radiative effect has not been underlined separately. The reason for the overlooked Df contribution is that dust intrusions are typically dominated by large particles (Dc), particularly near the source. Indeed, the DRE of total dust (DD) has been extensively studied in both SW and long-wave (LW) spectral ranges (e.g. Sicard et al., 2014a, 2016; Barragan et al., 2016; Granados-Muñoz et al., 2019), but only a few studies have addressed this issue by separating the dust components. For instance, Sicard et al. (2014b) reported that Dc particles primarily affect LW solar radiation, while Df is primarily responsible for SW radiative modulation. Furthermore, the impact of Df particles in the SW spectral range of the solar radiation is greater at the top of the atmosphere (TOA) than at the surface, i.e. the bottom of the atmosphere (BOA) (Córdoba-Jabonero et al., 2021), while the contribution of the Dc particles to the LW solar spectral range at the TOA is less than at the BOA (Sicard et al., 2022). Hence, it is essential to study Df and Dc particles separately when estimating individual radiative effects.

It should be highlighted that in dust climate models, the separation can be particularly helpful for analysing DRE per size mode. On one hand, when dust is treated as a singular entity (total dust), the contribution of fine dust can be masked, as the radiative properties input into radiative transfer models tend to be skewed towards the characteristics of coarse dust. On the other hand, recent studies have found evidence that Dc is inadequately represented in global atmospheric models (Adebiyi and Kok, 2020). In particular, the mass of coarse dust in the atmosphere could be about 4 times greater than simulated in current AeroCom climate models. Since Dc warms by absorbing both SW and LW radiation (Kok et al., 2017), the underestimation of Dc by climate models indicates that the net DRE could reflect more warming than has been previously estimated. Overall, the atmosphere appears to contain approximately 40 % more dust (both Df and Dc) than what is simulated by AeroCom models, accounting for about 80 % of the total particulate mass load in the atmosphere (Textor et al., 2007; Adebiyi and Kok, 2020). Thus, Df and Dc distinction is particularly relevant for data assimilation processes in climate models, where accurately representing the radiative properties of different dust components is critical to the improvement of model performance and predictions. As an added value, this study introduces the novelty of simulating SW dust DRE using two distinct approaches. On the one hand, the contribution of the Df and Dc components is simulated separately, and then the total dust DRE is estimated as their sum (DD = Df + Dc). On the other hand, the DRE is directly simulated for the total dust component.

The paper is organized as follows. The radiative transfer model applied and the parameterizations used in terms of SW range are described in Sect. 2. An overview of the dust outbreak can be found in Sect. 3. In Sect. 4 the results and discussion are shown. Finally, the main conclusions of this study are found in Sect. 5.

2 Methodology

2.1 Radiative transfer model: GAME

GAME (Global Atmospheric Model) is the radiative transfer model (RTM) used in this study to compute the solar radiation fluxes. In this work, they are simulated in the SW spectral range (0.2 to 4.0 µm, with a wavenumber resolution of 400 cm⁻¹ from 0.2 to 0.7 µm and 100 cm⁻¹ from 0.7 to 4.0 µm). Hereafter, for simplicity, the “SW” notation is omitted. GAME was developed by Dubuisson et al. (1996, 2004) and employs the discrete ordinates method (DISORT; Stamnes et al., 1988) to compute the solar flux at the boundary of the homogenous and plane atmospheric layers. The SW radiation flux is calculated for 18 layers of the atmosphere ranging from ground level to an altitude of 20 km.

This RTM accounts for the absorption of gases, specifically H₂O, CO₂, O₃, N₂O, CO, CH₄, and N₂, using the correlated k distribution (Lacis and Oinas, 1991). Additional information on the calculation of the gas transmission functions can be found in Dubuisson et al. (2004) and Sicard et al. (2014b). The parameterization of gas absorption is based on pressure, temperature, and relative humidity profiles. In particular, the data from the Global Data Assimilation System (GDAS) provided by the National Oceanic and Atmospheric Administration (NOAA) are used (ftp://arl.noaa.gov/archives/gdas1/; last access: 28 February 2023). Additionally, the model requires input from the spectrally resolved surface albedo (SA), which was obtained from MODIS (Moderate Resolution Imaging Spectrometer; https://modis.gsfc.nasa.gov; last access: 28 March 2023).

The input parameters, radiative properties as introduced in GAME for the SW spectral range, and the databases used are summarized in Table 1.

Table 1 Input parameters and radiative properties for the GAME model and databases in the SW spectral range. Note that z denotes the vertical dependence. DRE^(I) and DRE^(II) denote the particular approach used for DRE simulation, considering either dust component separation (Dc, Df) or total dust (DD, no separation), respectively.

2.1.1 AERONET properties used as input in GAME

In GAME, aerosols are included by comprehensively parameterizing their radiative properties in terms of the spectral dust extinction coefficient (α), single-scattering albedo (ω), and asymmetry factor (g). This section focuses on ω and g taken from NASA's AERONET (AErosol RObotic NETwork; http://aeronet.gsfc.nasa.gov, last access: 18 February 2023) Version 3 Level 2.0 data inversion products. At the time of this writing, few or no AERONET Level 2.0 data were available for all the stations. In such cases, AERONET Level 1.5 data were used instead. The AERONET g and ω are given with an uncertainty of ±0.03 and ±0.04, respectively, for dust particles (Dubovik et al., 2000, 2002).

Both ω and g values are assumed to be vertically constant and equal to the columnar values of ω and g. In addition, they are interpolated at the model wavelengths up to the largest SW wavelength provided by AERONET (1020 nm) and assumed to be constant beyond. Furthermore, they are hourly averaged and interpolated in the case that those hourly values are missing. Also note that g is given separately for the fine, coarse, and total modes.

2.1.2 Lidar-derived extinction used as input in GAME

Since the focus of this work is on obtaining the DRE while considering the separation between Df and Dc, continuous hourly averaged lidar measurements together with the POLIPHON method (Mamouri and Ansmann, 2017; Córdoba-Jabonero et al., 2018) were applied to obtain the vertical profiles of the dust extinction coefficient at 532 nm (α⁵³²) instead of using the aerosol optical depth given by AERONET. The POLIPHON approach relies entirely on empirically obtained dust input parameters, specifically the dust depolarization ratio and dust lidar ratio, without the need for a particle shape model. This method has undergone extensive validation and has been applied in several works (Mamouri and Ansmann, 2017; Córdoba-Jabonero et al., 2018, 2021; Ansmann et al., 2019; Couto et al., 2021; Salgueiro et al., 2021; Sicard et al., 2022). In the first step of POLIPHON, dust particles are separated from the remaining aerosols. In the second step, Df and Dc are identified, enhancing the analysis of the backscatter profile. Once separated by particle modes, α⁵³² for each component is obtained by considering the specific particle lidar ratio at 532 nm for each one (Ansmann et al., 2019). The uncertainties in the calculation of α⁵³² using this methodology are 30 %–50 %, 20 %–30 %, and 15 %–25 %, for Df, Dc, and total dust, respectively (Ansmann et al., 2019).

The coarse and fine lidar-derived α⁵³² profiles were previously obtained by López-Cayuela et al. (2023) (see Sect. 3). As those profiles have different vertical resolutions depending on the lidar system used, they have been downscaled to the 18 model layers (ranging from the surface up to 20 km height) through trapezoidal numerical integration in order to homogenize all the datasets. In addition, the lidar-derived α⁵³² profiles are used together with the AERONET Ångström exponent at 440–870 nm (AE^440−870) to derive the height-resolved dust extinction at other wavelengths within the SW range. Finally, hourly dust optical depth (DOD, τ) values for both Dc and Df modes and for the total dust (DD) are computed by vertically integrating the spectrally estimated dust extinction coefficient profiles.

2.2 Dust radiative effect estimation

The aerosol direct radiative effect (ARE) is obtained for dust, i.e. the dust direct radiative effect (DRE), which can be defined at a given height (z) as follows (Córdoba-Jabonero et al., 2021; Sicard et al., 2022):

$$\begin{array}{}\text{(1)}& \text{DRE}\left(z\right)=\left[F_d^{\downarrow }\left(z\right)-F_d^{\uparrow }\left(z\right)\right]-\left[F_{\mathrm{0}}^{\downarrow }\left(z\right)-F_{\mathrm{0}}^{\uparrow }\left(z\right)\right],\end{array}$$

where F_d and F₀ denote the solar radiative flux (W m⁻²) computed by GAME, with (d) and without (0) dust occurrence. The ↓ and ↑ arrows indicate whether the fluxes are downward or upward, respectively. The first $\left(\left[F_d^{\downarrow }\left(z\right)-F_d^{\uparrow }\left(z\right)\right]\right)$ and second $\left(\left[F_{\mathrm{0}}^{\downarrow }\left(z\right)-F_{\mathrm{0}}^{\uparrow }\left(z\right)\right]\right)$ terms of Eq. (1) represent, respectively, the net flux with and without the presence of dust.

According to that definition, a negative (positive) DRE value represents a cooling (warming) effect. In this study, the DRE is shown at both TOA and BOA. The dust contribution in the overall atmospheric column is evaluated by quantifying the atmospheric radiative effect, DRE(ATM), which is defined as follows:

$$\begin{array}{}\text{(2)}& \text{DRE}\left(\text{ATM}\right)=\text{DRE}\left(\text{TOA}\right)-\text{DRE}\left(\text{BOA}\right).\end{array}$$

Note that α profiles are obtained from “hourly” lidar measurements, and the DRE values denote the “instantaneous” simulated radiative values obtained from those hourly α profiles; hence, for simplicity, both terms are removed. The fine-to-total ratio (Df $/$ DD) of DRE (ftr_DRE) is also calculated, together with a linear fitting analysis of that variable over time. The slope of this linear fitting, denoted δDRE, represents a measure of the temporal rate of the relative Df contribution to the DRE. Beside this, the daily averaged DRE is computed as follows:

$$\begin{array}{}\text{(3)}& {\text{DRE}}^{\mathrm{daily}}={\displaystyle \frac{N_d}{\mathrm{24}}}{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\text{DRE}}^{\mathrm{hourly}}}{N}},\end{array}$$

where N is the number of hourly DRE values computed for each day, and N_d is the total number of hours with a solar zenith angle (SZA) < 90° over the course of the day (i.e. N_d = 12 in this work). Moreover, the episode-averaged DRE is computed as the average of the daily averaged DRE values for the entire event. For simplicity, hourly DRE will be denoted DRE hereafter, unless otherwise specified.

The dust radiative efficiency (DREff) is defined as the ratio of the DRE to DOD⁵³² values (in W m⁻² τ⁻¹) throughout the day. Hence, the daily DREff is obtained from the slope of the linear fitting of DRE values with respect to DOD⁵³² and forced to zero for each day.

2.3 Comparison between different approaches

The dust-induced DRE is calculated using two different approaches.

• 1.

  The previous one can be used, i.e. considering the contribution of the Df and Dc particles separately in the DRE, following the expression below:

  $$\begin{array}{}\text{(4)}& {\text{DRE}}^{\left(\mathrm{I}\right)}={\text{DRE}}^{\mathrm{DD}}={\text{DRE}}^{\mathrm{Df}}+{\text{DRE}}^{\mathrm{Dc}}.\end{array}$$
• 2.

  DRE is obtained from the contribution of the total dust (no separation):

  $$\begin{array}{}\text{(5)}& {\text{DRE}}^{\left(\mathrm{II}\right)}={\text{DRE}}^{\mathrm{total}}.\end{array}$$

This notation will be employed throughout the entire paper. That is, “DD” will refer to variables obtained by summing the Df and Dc components, and “total dust” will denote variables obtained without considering the mode separation.

The main distinction between the two approaches lies in certain radiative inputs used in the simulations. Specifically, in the first approach (DRE^(I)), both g and α are introduced as separate fine and coarse components, whereas in the second approach (DRE^(II)), these variables are introduced as their total component. All other variables are consistent in both approaches. Further details can be found in Table 1. Note that the second approach (DRE^(II)) is the most commonly used in the literature. However, the novelty of this study lies in showing that the first approach (DRE^(I)) enables dust DRE computation using the optical properties of each dust mode separately (Sicard et al., 2014b; Córdoba-Jabonero et al., 2021; Sicard et al., 2022), assuming that the dust-induced diffuse radiation from one of the modes does not interact with the other one. Moreover, the uncertainties in DRE^(I) should be considered higher than those in DRE^(II), as they are derived from the uncertainties in the calculation of the POLIPHON α⁵³², which are higher for Df and Dc modes than for total dust (see Sect. 2.1.2).

Finally, the differences found in the dust-induced DRE (ΔDRE) and the relative differences (Δ^relDRE) between the two approaches can be calculated using the following equations:

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }\text{DRE}={\text{DRE}}^{\left(\mathrm{I}\right)}-{\text{DRE}}^{\left(\mathrm{II}\right)},\end{array}$$

$$\begin{array}{}\text{(7)}& {\mathrm{\Delta }}^{\mathrm{rel}}\text{DRE}\left(\mathit{\%}\right)=\mathrm{100}{\displaystyle \frac{\left({\text{DRE}}^{\left(\mathrm{I}\right)}-{\text{DRE}}^{\left(\mathrm{II}\right)}\right)}{{\text{DRE}}^{\left(\mathrm{II}\right)}}}.\end{array}$$

Since the dust SW DRE is negative, negative (positive) ΔDRE (Δ^relDRE) values would indicate that DRE^(I) exhibits higher absolute values than DRE^(II). This suggests that a more pronounced cooling effect would be induced by dust according to the DRE simulations derived from the dust-mode separation approach (DRE^(I)) than from the total dust approach (DRE^(II)). Additionally, a statistical analysis based on the Δ^relDRE percentiles (P) will be considered to assess the significance of discrepancies between the two methodologies. In particular, the first (P25) and third (P75) quartiles will be examined.

3 Overlook of the dust outbreak crossing the Iberian Peninsula

During spring 2021, the IP experienced a significant and prolonged dust intrusion (25 March–7 April). This event was closely monitored and analysed by five stations strategically located across the IP, namely El Arenosillo/Huelva (ARN), Granada (GRA), Torrejón/Madrid (TRJ), and Barcelona (BCN) in Spain and Évora (EVO) in Portugal. Figure 1 shows the geographical position of each station in the IP. All these stations share the commonality of being dedicated to aerosol–cloud monitoring. Among several instruments, each station is equipped with an AERONET photometer (or is close to an AERONET station) and a lidar system. Those five stations share a common exposure to Saharan dust outbreaks, particularly during spring and summer months, albeit with varying frequencies (i.e. Córdoba-Jabonero et al., 2021; López-Cayuela et al., 2023). Moreover, each station exhibits a unique aerosol background. Specifically, the aerosol background at the ARN station is mostly of marine and rural origin, as ARN is located in a rural environment in the southwestern IP and is less than 1 km from the Atlantic coastline. The EVO station is located in a rural region with limited industrialization and low levels of anthropogenic aerosol concentrations (Pereira et al., 2009; Preißler et al., 2013). Both the GRA and TRJ stations are located in populated cities, and their background aerosols are of anthropogenic origin (Lyamani et al., 2012; Molero et al., 2014). Finally, the BCN station is located on the northeast coast of the IP, within a densely populated and industrialized region; thus the background aerosol load is predominantly composed of urban and marine aerosols (Sicard et al., 2011).

Figure 1 MODIS Terra-/Aqua-corrected reflectance over the Iberian Peninsula on 31 March 2021. The five lidar stations are marked by yellow pins. From southwest to northeast the stations are as follows: El Arenosillo/Huelva (ARN), Évora (EVO), Granada (GRA), Torrejón/Madrid (TRJ), and Barcelona (BCN) sites. The coordinates and altitudes are also included.

Focusing on the dust intrusion studied in this work, a deep analysis of the synoptic situation and back-trajectory analysis revealed that the dust intrusion originated in the Saharan region and traversed the IP from the southwest to the northeast. The dust event was accompanied by extended cloud cover across all five Iberian lidar stations, resulting in the unavailability of some lidar retrievals. Consequently, the dataset exhibits certain gaps, along with periods of no lidar measurements.

An extensive study of the dust outbreak can be found in López-Cayuela et al. (2023), together with the methodology to obtain the optical and mass properties derived from the polarized lidar measurements. This section aims to summarize the most significant results for the optical properties, which are essential for the development and discussion of the current study.

This dust event is particularly noteworthy due to not only the intensity of the dust outbreak but also the prevailing meteorological situation conditioning the main findings. Among those were (1) the detection of dust in the upper layers of the troposphere in the northern half of the IP (reaching up to 10 km in height) caused by the atmospheric instability that induced vertical motion of the air parcels and (2) the absence of uniform gravitational settling observed throughout the IP; that is, the Df $/$ DD proportion remained nearly constant along the dust pathway across the IP (López-Cayuela et al., 2023), which contrasts with the long-range dust events observed in central Europe.

An overview of the temporal evolution of the dust outbreak crossing the IP in terms of the extinction (α, km⁻¹) over the five Iberian lidar stations is shown in Fig. 2. In particular, the hourly mean DD α at 532 nm (${\mathit{\alpha }}_{\mathrm{DD}}^{\mathrm{532}}$) is shown. Despite some data gaps due to either no inversion possible or no lidar measurements, the transport of dust particles can be appreciated by looking at the behaviour of ${\mathit{\alpha }}_{\mathrm{DD}}^{\mathrm{532}}$ along the southwest IP (Fig. 2e) to the northeast IP (Fig. 2a) pathway of the dust intrusion during the overall dusty period (25 March–7 April 2021). In terms of the episode averaging, DD DOD⁵³² decreases as latitude increases, ranging from 0.34 at ARN to 0.14 at BCN, without significant differences in the Df $/$ DD ratio (ftr_DOD) between the stations (approximately 30 %; see Table 2).

Figure 2 Temporal evolution of the DD extinction coefficient (${\mathit{\alpha }}_{\mathrm{DD}}^{\mathrm{532}}$, km⁻¹) at the five Iberian lidar stations (from the northeast to southwest IP in order of decreasing latitude): (a) Barcelona (BCN), (b) Torrejón/Madrid (TRJ), (c) Évora (EVO), (d) Granada (GRA), and (e) El Arenosillo/Huelva (ARN). Profile gaps correspond to either no inversion available or no lidar measurements.

Table 2 Dust optical depth at 532 nm (DOD⁵³²) for the fine (Df) and coarse (Dc) dust components and DD (DD = Df + Dc) at the five lidar stations (from left to right as latitude increases). ftr_DOD denotes the Df-to-DD ratio (in %). $\stackrel{\mathrm{‾}}{X}$ indicates the average of the whole dust event (the standard deviation is also shown), and X^max refers to the maximal hourly value within the dusty period.

In this work, the hourly averaged DOD⁵³² values are considered low when DOD⁵³² < 0.20, moderately intense when 0.2 ≤ DOD⁵³² < 0.50, intense when 0.50 ≤ DOD⁵³² < 1.00, and extreme for DOD⁵³² ≥ 1.00. Thus, this episode spans the full range, highlighting high and extreme DD DOD⁵³² values on the days of maximal incidence (29 March–1 April). The mean and maximum hourly averaged DOD⁵³² values for the Df, Dc, and DD components during the entire episode are shown in Table 2.

4 Results and discussion

This section is divided into an analysis of the specific radiative dust properties introduced in the model (Sect. 4.1), the DRE estimated considering the first approach (Sect. 4.2), and the discrepancies found by estimating the DRE from the two approaches considered (Sect. 4.3). GAME calculations are performed when the Sun is above the horizon (SZA < 90°). The DRE is calculated using Eqs. (1) and (2) for the Df, Dc, and total dust components separately (see Sect. 2.2).

4.1 Radiative properties during the dust episode

4.1.1 Surface albedo (SA)

For each station considered in this study, the SA values corresponding to the dust episode reveal a symmetric distribution (regardless of the wavelength under consideration). In fact, for each wavelength, the relative differences in SA throughout the episode at the same station are less than 0.1 %. For this reason, although hourly values are considered for the simulations, Fig. 3a represents the mean value of SA as a function of the wavelength and station.

Figure 3 Spectral behaviour of (a) the mean of the MODIS surface albedo (SA) distribution; (b) the AERONET asymmetry factor (g) for fine, coarse, and total modes; and (c) the AERONET single-scattering albedo (ω) at the five Iberian lidar stations: Barcelona (BCN, blue), Torrejón/Madrid (TRJ, orange), Évora (EVO, yellow), Granada (GRA, purple), and El Arenosillo/Huelva (ARN, green). The date and time chosen are 1 April at 10:00 UTC for the BCN station and 31 March at 12:00 UTC for the rest of the sites.

However, regarding the SA values (Fig. 3a), a similar pattern across all sites is found, with SA increasing with wavelength until 870 nm, reaching the maximum values at that wavelength. Within the visible range (VIS, wavelengths < 780 nm) the values increase with latitude, showing SA⁴⁸⁰ = 0.04 at ARN and SA⁴⁸⁰ = 0.07 at BCN. This characteristic is not preserved at the near-infrared range (NIR, 780–2500 nm). EVO and TRJ show maximum SA values of 0.33–0.36 at 870 nm, which ranges from 0.26 to 0.20 for the other stations. Along the rest of the NIR range, the SA values decrease until SA${}_{max}^{\mathrm{2130}}$ = 0.13–0.18. Notably, the SA values in this study align with those found in Granados-Muñoz et al. (2019) and Córdoba-Jabonero et al. (2021) for GRA and BCN, respectively, in summertime.

Larger SA values indicate that, at a given wavelength and constant incoming radiation reaching the surface, more radiation will be reflected upward. Overall, minimum SA values are found at BCN, while maximum ones are found at the TRJ and EVO stations. The differences between the five lidar stations are lower than 5 % in the VIS range and for wavelengths greater than 1650 nm. At the beginning of the NIR range (870–1240 nm), EVO and TRJ reflect 7 %–15 % more radiation than the other stations. As shown in Sect. 3, the stations are located in a variety of places, from rural sites with different types of vegetation to urban sites. Thus, variations in SA values between stations are expected, as SA depends on the surface coverage.

4.1.2 Asymmetry factor (g) and single-scattering albedo (ω)

Conducting the same analysis for the parameters g and ω, a skewed distribution is obtained for each station and wavelength, with maximum relative differences during the dust episode of 1.5 % and 2.5 %, respectively. For this reason, and although once again the hourly averaged value of these variables is used for each simulation, Fig. 3b–c depict an hourly data snapshot from the days of highest incidence, i.e. 1 April at 10:00 UTC for the BCN station and 31 March at 12:00 UTC for the rest of the sites. The g and ω histograms for the five lidar stations can be found in the Supplement (Figs. S1–S5).

The g values decrease with the wavelength, as the forward scattering is approximately 1.3 times greater for the coarse mode (coarse g⁴⁴⁰ = 0.82–0.86) than for the fine mode (fine g⁴⁴⁰ = 0.67–0.70). As explained in Córdoba-Jabonero et al. (2021), this result implies that the solar radiation scattered to the surface is greater for the coarse mode than for the fine mode, keeping a constant aerosol optical depth (AOD) and at low SZA, independently of the wavelength. Regarding the total g, its values range between those for the fine and coarse modes (total g⁴⁴⁰ = 0.75–0.80). The variability in ω is also low between the stations, increasing with wavelength. In particular, ω⁴⁴⁰ ranges from 0.90 to 0.93, representative of weakly absorbing particles. From 675 nm to NIR, ω increases to values of 0.98–0.99 for all stations.

Figure 4 shows the temporal evolution of g⁴⁴⁰ for the fine and coarse components, together with both g⁴⁴⁰ and ω⁴⁴⁰ for the total component at the five lidar stations. Along the entire dust outbreak, these properties remain relatively constant at each station. The variability observed is minimal, with values below 5 % and 3 % for g⁴⁴⁰ and ω⁴⁴⁰, respectively. It should be noted that, in this study, the variability has been quantified as the percentage ratio between the standard deviation and the mean value along the entire episode for each station. Specific values are detailed in Table 3.

Figure 4 Overview of the temporal evolution of the radiative properties of fine, coarse, and total g and for ω at 440 nm at the five lidar stations, i.e. (a) BCN, (b) TRJ, (c) EVO, (d) GRA, and (e) ARN.

Table 3 Radiative properties at 440 nm at the five lidar stations (from left to right as latitude increases). In particular, the parameters are the fine (Df) and coarse (Dc) dust components of the asymmetry factor (g) and the total dust component for both g and single-scattering albedo (ω). The symbol $\stackrel{\mathrm{‾}}{X}$ indicates the average of the whole dust event (the standard deviation is also shown), and X^max and X^min refer to the maximum and minimum hourly values within the dusty period, respectively.

As mentioned before, the coarse component of g has the most influence on scattered radiation. The lowest values of this parameter were observed at ARN and BCN, while the highest values are evident at the GRA station. Thus, the solar radiation scattered to the surface is lower for ARN and BCN than for the GRA station. Conversely, the pattern is reversed for ω. However, the differences among the stations are marginal, displaying relative differences in the episode-average values of less than 3 % for the total component in g⁴⁴⁰ and in ω⁴⁴⁰. Regarding the separation of the dust components, the differences in the episode-averaged g⁴⁴⁰ were 1.5 % and 4 % for Df and Dc, respectively.

In summary, the results for g and ω are similar for all the stations, showing the typical values expected for dust (Dubovik et al., 2002). Indeed, the differences between the stations are less than 5 %, regardless of the variable and wavelength being examined. Thus, the differences expected in DRE between stations should mostly rely on the specific values of SA for each station (Sect. 4.1.1) and the variability in dust extinction between the stations and during the dusty episode. A summary of the evolution of the dust extinction can be found in Sect. 3, and a complete analysis can be found in López-Cayuela et al. (2023).

4.2 Dust-induced direct radiative effect

This section is devoted to describing the episode in terms of the evolution of the dust properties and DRE across the stations considered. Figure 5 shows the temporal evolution of the DRE and the ftr_DRE at ARN station as an example. To avoid an excessively lengthy paper, the graphics for the rest of the stations can be found in the Supplement (Figs. S6–S9). Figure 6 shows DRE as a function of DOD⁵³² at ARN together with the DREff values obtained for all five lidar stations at BOA and TOA. Moreover, the specific episode-averaged and maximum DRE values together with DREff, ftr_DRE, and δDRE for each station are shown in Table 4.

Figure 5 Dust direct radiative effect (DRE, W m⁻²) at BOA (purple), TOA (green), and in the atmosphere (ATM, yellow) at the ARN station as induced by the (a) fine dust (Df) and (b) coarse dust (Dc) particles. Daily mean values are also included above (keeping the same colours). (c) Fine-to-total (Df $/$ DD) ratio of the dust DRE (ftr_DRE) at BOA (purple) and TOA (green), together with their linear fits between 25 March and 7 April. The standard deviation of the slope and intercept are shown in brackets. The absolute temporal rate of the ftr_DRE ratio (δDRE) at BOA (TOA) is +0.026 % h⁻¹ (+0.023 % h⁻¹), equivalent to an increase of +0.62 % d⁻¹ (+0.55 % d⁻¹).

Table 4 Episode-averaged hourly values of the dust radiative effect (DRE, in W m⁻²) and the standard deviation on the surface (BOA) as induced by fine (Dc), coarse (Dc), and total dust (DD) at the five Iberian lidar stations. $\stackrel{\mathrm{‾}}{X}$ indicates the average of the whole event (standard deviations are also shown), and X^max refers to the maximal value. DREff (in W m⁻² τ⁻¹) denotes the DRE efficiency. ftr_DRE denotes the Df-to-total dust DRE ratio (in %), and δDRE (in % d⁻¹) is the slope of the linear fitting analysis of hourly averaged ftr_DRE in time, considering a minimum of six ftr_DRE values each day.

4.2.1 DRE at BOA

The DRE is negative at all stations during the entire dusty period, indicating dust-induced cooling (see Figs. 5a and S6a–S9a, in purple). The intensity of the episode in terms of the radiative forcing is directly correlated with the DOD⁵³², which is greater for Dc compared to Df particles. Thus, during the first part of the episode (until 1 April) the dust event was moderate, showing hourly Df (Dc) DRE values above −10 W m⁻² (−20 W m⁻²). During the rest of the episode, the dust event was weaker, showing hourly DRE values 50 % lower. On some specific days (27 and 31 March), the event was notably intense at the ARN and TRJ stations (mean AOD ∼ 0.80), reaching hourly Df (Dc) DRE values below −40 and −20 W m⁻² (−60 and −40 W m⁻²), respectively. Indeed, the maximum hourly DRE for the dust episode is found at ARN (DOD⁵³² = 1.30), while the rest of the stations exhibit maximum values of approximately 30 %–50 % lower (Table 4).

The intensity of the overall dust episode can be understood by averaging the daily DRE values (episode-averaged DRE; see Sect. 3.2). The highest episode-averaged DRE values for DD are found at the southern and central stations, decreasing by 40 % at EVO and 50 % at BCN (see Tables 2 and 4), with differences mainly due to the decrease in DOD with latitude. As expected, the differences found between stations are mainly related to the DOD decrease with latitude (see Tables 2 and 4). Indeed, Meloni et al. (2005) found that the radiative forcing at both the TOA and BOA varies almost linearly with AOD, and the strong dependence of the DRE on AOD is well documented (Prasad et al., 2007; Sicard et al., 2014b; Lolli et al., 2018; Meloni et al., 2018; Granados-Muñoz et al., 2019).

DRE values are compared with other authors' findings over the IP. Córdoba-Jabonero et al. (2021) studied a moderately intense dust episode over BCN (maximum hourly averaged DOD⁵³² of 0.60) by separating the Df and Dc components. The daily averaged DRE results are similar for equivalent daily averaged DOD⁵³² conditions. By considering other studies when the Df and Dc separation was not applied under similarly intense dust aerosol load conditions, the daily total DRE reported in those studies is almost 50 % greater than in the present study (Cachorro et al., 2008; Valenzuela et al., 2017). However, similar results are found for extreme dust aerosol conditions (Bazo et al., 2023).

This study also analyses the DREff for Df, Dc, and DD particles. The results show an excellent degree of confidence, with a correlation coefficient greater than 0.80 at all stations. By averaging those DREff values for the five lidar stations, a mean value of −138.5 ± 20.3 W m⁻² τ⁻¹ (−95.3 ± 5.5 W m⁻² τ⁻¹) is obtained for Df (Dc) particles. However, when considering the stations separately, the DREff slightly increases with latitude (i.e. from ARN to BCN), up to a maximum of 10 % for all particle modes (see Table 4 and Fig. 6). The GRA station deserves special mention, as it is located at a similar latitude as ARN and has similar dust optical properties (see Fig. 3) and Df $/$ DD proportions but shows 5 %–10 % lower DREff values for Df and Dc particles. This difference may be due to the limited data available in the study, resulting in a smaller range of DODs. Comparing our results with other studies over the IP, similar results are found, with maximum differences in DREff of 10 % (Obregón et al., 2015; Sicard et al., 2016; Granados-Muñoz et al., 2019; Córdoba-Jabonero et al., 2021).

Figure 6 Dust direct radiative effect (DRE) as a function of DOD⁵³² at ARN station (left panels), together with the DREff values obtained for each station (right panels) at (a) BOA and (b) TOA. Both are shown separately for fine dust (Df, blue circles), coarse dust (Dc, dark blue squares), and total dust (DD, red triangles). The error in DREff is shown in brackets in the left panels and by bars in the right panels. The period considered is from 25 March to 7 April 2021.

Finally, and regarding the contribution of the Df particles, the mean daily ftr_DRE at the five lidar stations is around 40 %, indicating that Df particles contribute almost half of the total DRE at BOA. Although there is a slight increase in δDRE during the dust outbreak period, it remains below +1 % d⁻¹ (Fig. 5c, Table 4).

4.2.2 DRE at TOA and in the atmosphere

This study shows that the DRE at TOA is negative at all stations, indicating a dust-induced cooling effect (see Figs. 5 and S6–S9). Similar to the results found at BOA, the first part of the dust outbreak was moderate in terms of the radiative forcing (Córdoba-Jabonero et al., 2021) and weaker for the rest of the episode, showing slightly lower values than at BOA. Once again, the greatest hourly DRE values were found at ARN and TRJ on 27 and 31 March (see Figs. 5 and S6 and Table 5). Overall, the stations exhibit maximum daily DRE values of approximately 35 %–70 % lower than at ARN. What is new is that the highest episode-averaged DRE values for DD are found at the southern and central stations, decreasing 50 % at EVO and BCN (see Table 5). In comparison with other studies in similar dust load conditions over the IP, comparable daily averaged DRE results are found for Df and Dc particles (Córdoba-Jabonero et al., 2021) but are 20 %–30 % lower for DD particles (Cachorro et al., 2008; Valenzuela et al., 2017).

Table 5 The same as Table 4 but for the top of the atmosphere (TOA).

The linear regression analysis shows a slight positive trend in the contribution of Df particles to total dust DRE during dust outbreaks, with higher δDRE values at TRJ. However, the GRA station shows negative values of almost 3 % d⁻¹ (see Table 5). The episode-averaged ftr_DRE shows similar values of approximately 45 % at all lidar stations, suggesting that Df particles can induce almost half of the total dust DRE at TOA.

As for the BOA analysis, the linear fitting between DOD and DRE shows a good degree of confidence, with a correlation coefficient greater than 0.80. By averaging those values for all five stations, a mean DREff of −107.1 ± 10.4 W m⁻² τ⁻¹ (−57.5 ± 3.3 W m⁻² τ⁻¹) for Df (Dc) particles is found, highlighting the higher variability in Df compared to Dc particles, as indicated by the standard deviation. The study of DREff allows the assessment of the radiative effect by eliminating the aerosol burden variable. However, the aerosol load is significant for not only DRE calculations but also the structure and altitude at which the dust is located. The most notable difference among the studies from the literature is the dust layering. In this work, the dust particles reached the higher troposphere (over 6 km above sea level) several times at certain stations (TRJ and BCN). Thus, DRE is more pronounced, and DREff is greater at the DOD constant for those cases than a dust layer located at lower altitudes. For this reason, our results could differ from other studies, presenting values 9 % and 30 % greater for Df and Dc DREff (Córdoba-Jabonero et al., 2021), up to 3 times greater for DD DREff (Granados-Muñoz et al., 2019), or even 15 % lower for DD particles (Sicard et al., 2016).

In contrast with the DRE at BOA and TOA, the DRE (ATM) is positive for the entire episode, indicating an atmospheric dust-induced heating, as DRE(BOA) > DRE(TOA) in absolute values (see Eq. 2). Regarding the episode-averaged DRE, the highest values are found at ARN and GRA and decrease with latitude, i.e. up to 50 % at EVO, 25 % at TRJ, and 65 % at BCN (see Table 6). Moreover, the maximum hourly averaged DRE value is found at ARN, and results are 55 %–80 % lower for the rest of the stations (Table 6). Similar results are found in terms of the atmospheric ftr_DRE ratio, showing values approximately in the 27 %–33 % range.

Table 6 The same as Table 4 but for the atmosphere (ATM).

4.3 Differences between approaches for DRE estimation in SW radiative flux simulations

4.3.1 Potential limitations in the radiative properties

As described in Sect. 2.3, a comparative analysis is performed for the dust-induced DRE obtained using two different approaches. That is, either the separated contribution of each of the two dust modes (DRE^(I)) or the contribution of the total dust overall (DRE^(II)) is considered. Note that the second approach is assumed to be the reference in the comparison (see Eqs. 6 and 7), as it is widely used by other authors. Potential limitations of this methodology should be considered, as they could affect the differences in DRE simulations found between the approaches used in this study.

First, the particular sensitivities to coarse and super-coarse dust particles should be taken into account, as they can be differentiated by passive sensors (sun/sky photometers) and active ones (lidar systems). Specifically, the properties of particles with radii exceeding 15 µm are not retrieved by the AERONET inversion algorithm. However, Ryder et al. (2019) demonstrated that dust particles with a radius greater than 15 µm contribute only 1 %–3 % to the particle extinction coefficient at 550 nm. This indicates that the size cutoff effect in the AERONET data inversion procedure has a negligible impact on the AERONET inversion products. Thus, the coarse components regarded by AERONET and lidar inversion are comparable despite the potential different sensitivities of both sensors.

It is worth remembering another important difference when applying both approaches in the DRE simulations. In DRE^(I), g and α are introduced for the fine and coarse modes separately, whereas in DRE^(II), these variables are introduced for the total dust. Oppositely, AERONET ω is provided only for the total mode (no separation; see Table 1). To assess the impact of using a single total ω instead of separate components, the parameterizations proposed by Wang et al. (2019) were applied. In that study, an efficient approach to the parameterization of dust aerosols in SW bands is provided. The differences between the coarse ω obtained using the parameterizations and total ω of AERONET are less than 1 % on average, except for the 440 nm wavelength, which rose to 3 %. Those results are aligned with AERONET's own uncertainty in ω. Thus, using total ω in both Df and Dc DRE^(I) simulations should not have a significant impact.

As far as the fine mode is concerned, the AERONET fine g value is introduced in GAME to compute the DRE related to the Df component, as shown in Sect. 2.1.1. However, the fact that the fine g is influenced by the total fine mode, i.e. both Df and background aerosols, should be taken into account. In addition, assuming a vertically uniform g for the Df component could have substantial consequences, mainly because fine g values can be strongly affected by background aerosols, which are also mostly confined to the boundary layer. Therefore, a complementary study has been conducted to study the degree of suitability of applying AERONET fine g values to DRE computation of the Df component. This way, hourly fine g values for those cases reported under dust-dominated conditions (i.e. cases where the AERONET fine-mode fraction is less than 40 %) were compared with respect to the fine g values reported for all cases during the study period (25 March to 7 April 2021). In summary, results indicated that differences in fine g between cases under dust-dominated conditions and all conditions at each station were not significantly high for those obtained during the study period and neither were the period-averaged fine g differences between stations. Specifically, the relative differences ranged between −3.4 % and +0.4 %. Therefore, the contribution of background aerosols to the fine g can be assumed to be negligible, and its values can be used for the Df component. This outcome could be, to some extent, expected since background aerosols generally exhibit very low linear depolarization ratios, as they are predominantly small and spherical particles with a minimal contribution to the g parameter. Hence, at least in this specific dust event, the assumption of a constant g value for the Df component (i.e. AERONET fine g) throughout the entire atmospheric column can be considered reliable.

4.3.2 Comparative analysis

The comparison of the DRE computation at BOA and TOA obtained from the two approaches is based on examining their absolute (ΔDRE; W m⁻²) (see Eq. 6) and relative (Δ^relDRE; %) (see Eq. 7) differences. On the one hand, Fig. 7 shows Δ^relDRE as a function of SZA from 25 March to 7 April 2021 for all five lidar stations involved in this study. The DD DOD⁵³² is also introduced, although a clear dependence is not observed. First, note that large Δ^relDRE values are observed between methodologies at high SZA (> 70°), which is expected as the intrinsic uncertainty in GAME simulations increases at these angles due to the assumption of a plane-parallel atmosphere in the model (see Sect. 2.1). Therefore, data at these SZAs (> 70°) should be discarded. Indeed, data percentiles P(75) = + 0.65 % and P(25) = −0.85 % are obtained at both BOA and TOA, that is, only 25 % of Δ^relDRE present values below −0.85 % when considering all SZAs. However, those percentiles decreased at SZA < 70°; i.e. P(75) = + 0.02 % (−0.15 %) and P(25) = −6.8 % (−10.9 %) are computed at BOA (TOA) (see Fig. 7). These results indicate that lower absolute DRE^(I) values with respect to DRE^(II) (note that SW DRE is negative) are predominantly found. Therefore, dismissing the comparison with SZA > 70°, DRE^(I) values on average are lower (in absolute units) with respect to DRE^(II) (i.e. DRE^(I) shows less-negative values than DRE^(II)), with mean Δ^relDRE around −3 ± 5 % and −6 ± 10 % on average at BOA and TOA, respectively. Therefore, DRE estimations that consider the separation of the dust Df and Dc components with respect to total dust would result in a less pronounced dust-induced cooling effect.

Figure 7 Relative differences in DRE (Δ^relDRE, in %) obtained from the two approaches as a function of SZA at (a) TOA and (b) BOA for all five lidar stations involved in this study from 25 March to 7 April 2021. The dependence on DD DOD⁵³² is shown as a colour-scaled bar. P(25) and P(75) are marked by horizontal short-dashed and dashed lines, respectively. The vertical dashed grey line denotes SZA = 70°.

On the other hand, Fig. 8 shows the absolute differences in DRE obtained from the two approaches (ΔDRE; see Eq. 6) as a function of DRE^(II) at both BOA and TOA for all five lidar stations. The dependence on the DD DOD⁵³² is also highlighted. For SZA < 70°, approximately 75 % of ΔDRE values are positive, indicating that absolute DRE^(I) values are lower than the DRE^(II) ones at both BOA and TOA. Indeed, P(25) = −0.01 (+0.03) W m⁻² at BOA (TOA) (see Fig. 8) is computed. Hence, those results indicate, overall, that the DRE estimation considering the separation of both dust components would result in a less pronounced cooling effect than that of the total dust (note the SW dust DRE is negative; i.e. dust induces a cooling effect at both BOA and TOA), as also stated previously. However, it should be noted that absolute ΔDRE tends to increase as DOD⁵³² increases; in general, ΔDRE was mostly close to zero at lower DOD (< 0.2), increased somewhat at moderate-/high-dust-load conditions (DOD > 0.50), and reached maximum values of nearly −20 W m⁻² at both BOA and TOA. In particular, when DRE^(II) > −40 W m⁻² (shown by a dashed grey line in Fig. 8), mean ΔDRE values are +0.4 ± 0.8 and +0.3 ± 0.8 W m⁻² for BOA and TOA, respectively; meanwhile for DRE^(II) ≤ −40 W m⁻², corresponding to higher dust load conditions, ΔDRE shows mean values of −2.8 ± 4.0 and −3.7 ± 4.3 W m⁻² at BOA and TOA, respectively. Indeed, this is also reflected by examining the significant quartiles (median) computed for ΔDRE when DRE^(II) < −40 W m⁻²: P(75) (P(50)) = 0.17 (−0.1) and −1.7 (−3.3) W m⁻² at BOA and TOA, respectively. These opposite results, depending on dust loading, suggest that dust induces a more pronounced cooling effect at both BOA and TOA under moderate-to-high-dust conditions when the contributions of the Dc and Df components are considered separately. Furthermore, this cooling effect is more evident at TOA than at BOA, probably due to the presence of dust at higher altitudes than those typically found in the troposphere.

Figure 8 Differences in DRE (ΔDRE; W m⁻²) obtained from the two approaches (DRE^(I)–DRE^(II); see Eq. 6) at (a) TOA and (b) BOA for all five lidar stations involved in this study from 25 March to 7 April 2021. The dependence on the DD DOD⁵³² is shown as a colour-scaled bar. P(25) is marked by a horizontal dashed line. Data for SZA > 70° are highlighted by cross symbols. The vertical dashed grey line denotes the point at which DRE^(II) = −40 W m⁻².

The differences found between DRE^(I) and DRE^(II) can be explained by the uncertainties inherent to the radiative properties (g, ω) introduced as inputs to the RTM. Andrews et al. (2006) determined with a sensitivity analysis that a 10 % decrease in total g would correspond to 13 % and 19 % reductions in BOA and TOA forcing, respectively. On the one hand, AERONET total g and ω are given with an uncertainty of ±0.03 and ±0.04 (approximately 4 % for g and ω), respectively, for dust particles (Dubovik et al., 2000, 2002). On the other hand, it is worth examining the differences between the total and coarse g, considering the fact that the contribution of Dc plays a significant role in DRE in cases with medium and high DOD. In particular, the coarse g used in Dc simulations is, on average, 8 % to 10 % greater than total g. Thus, the results found in our study agree with those obtained by Andrews et al. (2006).

Moreover, the role of Df cannot be disregarded in DRE estimations. In this work, the impact of Df particles on DD DRE has been highlighted, even in cases with medium and high DOD when the Dc contribution plausibly dominates. It is possible that the Df contribution to the total DRE could be systematically underestimated by applying the DRE^(II) approach in relatively high-dust environments.

5 Summary and conclusions

This work serves as a complement to the study conducted by López-Cayuela et al. (2023), where the vertical distribution of both the fine- and coarse-mode dust during an intense dust episode that took place between 25 March and 7 April 2021 crossing over the Iberian Peninsula was described in detail. Now, this study focuses on the estimation of temporal variations in the instantaneous and daily direct radiative effects of dust in the short-wave spectral range during the same dust event. The investigation was carried out using data obtained from five lidar stations situated in El Arenosillo/Huelva, Granada, Torrejón/Madrid, and Barcelona in Spain, as well as Évora in Portugal, whose observations were also used in that previous study. The main findings are summarized next.

Across all the stations examined in this study, the DRE remains negative throughout the dusty period at both BOA and TOA, indicating aerosol-induced cooling at these levels. Conversely, the DRE shows a positive trend, indicating an aerosol-induced heating. The results are similar to those reported by a wide range of studies that can be found in the literature.

On some specific days, DRE values were quite high, corresponding to DODs greater than 0.80. In particular, this work shows maximum hourly averaged DRE values of −54.1 W m⁻² (−96.3 W m⁻²) and −35.9 W m⁻² (−53.5 W m⁻²) induced by Df (Dc) particles at BOA and TOA, respectively. Regarding the impact of Df particles, they induce almost half of the total DRE. In particular, Df contributions of 40 % and 45 % are found at BOA and TOA, respectively. Thus, despite dust intrusions being primarily governed by the coarse dust mode, the contribution of the fine-dust mode plays a relevant role in DRE estimations.

In terms of DREff, similar values are found at every station, with variations less than 10 % between them. By comparing with other authors, this study shows DREff values slightly higher for Df particles (< 10 %) at both BOA and TOA. However, the difference increases for Dc particles, with DREff values being 15 % and 30 % greater at BOA and TOA, respectively, than those reported by other work. In terms of DD particles, the DREff results are similar to those reported at BOA and higher than those at TOA.

DRE simulations based on the separation of dust modes (after disregarding values for SZA > 70°) show that nearly 75 % of absolute DRE values are lower with respect to those obtained by considering the total dust (denoted DRE^(II) in this work). That is, the differences between the two approaches are mostly positive. This indicates that a less pronounced cooling effect would be induced overall by the DRE estimations based on the separation of the two dust components. However, the mean relative differences between approaches are not significant, i.e. −3 % and −6 % on average for SZA < 70° at BOA and TOA, respectively.

However, opposite results are obtained under moderate-to-high-dust-load conditions, with differences between approaches, as described in this work, mostly showing negative values. By separating the two dust components, DRE estimations suggest that dust induces a more pronounced cooling effect. This may highlight the important role of Df particles in DRE simulations, which are likely underestimated when total dust (no separation) is considered in relatively high-dust environments. In this study, the effect is more evident at the TOA than at the BOA, probably due to the presence of dust at higher altitude levels than those usual in the troposphere. This fact can introduce relevant implications for radiometric measurements from satellite instrumentation.

For future work, it would be relevant to perform the same study but under dusty conditions with variable ftr_DOD to comprehensively assess the role of the fine dust particles in DRE using two different conceptual approaches. This would allow the determination of how dust ageing affects the DRE, considered by looking at the proportion of Df particles relative to the total dust (DD). It is important to note that the methodology used in this work combines lidar data, which provide vertical aerosol profiles, with photometer data, which derive columnar aerosol radiative properties. One way to improve this methodology would be to use height-resolved key radiative parameters such as ω and g profiles, which would be particularly useful in cases where either background aerosols play a significant role or there is a mixture of aerosol types. The use of multiwavelength lidars instead of single-wavelength lidar systems, together with particular inversion methods, could provide the required ω and g profiles.