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

The dust direct radiative effect (DRE) in the long-wave (LW) (DRE_LW), and the net effect (DRE_NET), is analysed during an intense and long-lasting Saharan dust intrusion over the Iberian Peninsula, complementing the study on the short-wave (SW) DRE (DRE_SW) (López-Cayuela et al., 2025). In LW, a warming effect at both the bottom-of-atmosphere (BOA) and the top-of-atmosphere (TOA) levels is induced by the fine (Df) and coarse (Dc) dust particles, while Dc being dominant. The DRE_LW-to-DRE_SW ratio for Df ranged 4 %–8 % at BOA (1 %–4 % at TOA), and for Dc it was rather higher (39 %–54 % at BOA and 20 %–50 % at TOA). DRE_NET was consistently negative (net cooling) at both levels, and the derived atmospheric DRE_NET was positive (net warming). The Df contribution to DRE_NET was 12 % (LW) and 30 % (SW). The SW aerosol heating rate (AHR) peaked at higher altitudes, inducing warming within the dust layer, than LW AHR (weaker cooling). Consequently, a net warming inside the dust layer was found, with potential cooling below and above. While SW dominates the net atmospheric warming, LW cooling partially mitigates it. As a novelty of this study, two methodologies for estimating DRE in both LW and net spectral ranges are compared. Differences in DRE between a classical approach considering total dust and an approach separating fine and coarse modes are analysed. DRE_LW (and DRE_NET) is underestimated (overestimated) by using the dust-mode separation approach in comparison to the classical one (no separation) when fine radii are lesser (greater) than a particular threshold (e.g. 0.1 µm), revealing the particle size impact in DRE_LW. The dust-induced net effect is primarily driven by SW and modulated by LW. The classical (no separation) approach overestimates DRE_NET, with mean (standard deviation) relative differences of −5 % (7 %) at BOA and −9 % (13 %) at TOA. Moreover, under moderate-to-high dust, separating Df and Dc contributions yields a weaker (stronger) net cooling at BOA (TOA).

1 Introduction

The latest report about airborne dust from the World Meteorological Organization (WMO Bulletin, 2023) reveals that the global surface dust concentration has seen a slight increase in 2022 compared to 2021. This is attributed to increased emissions from several dust-active sources, such as West-Central Africa. Among the most affected regions, receiving a dust influx much greater than the climatological mean, the Iberian Peninsula is prominently featured. Particularly, the anomaly of the annual mean surface dust concentration in 2022 (relative to the 1981–2010 mean) shows increased values of 5–20 µg m⁻³. Those results agree with several different studies indicating that since the pre-industrial Era there has been a 46 % increase in the mass of dust lifted into the atmosphere in these North African regions (Kok et al., 2023). Notably, there is not only a rising frequency of Saharan dust episodes in the Iberian Peninsula compared to long-term historical data (Sousa et al.,2019), but also an increasing number of reports on extreme and highly intense episodes (e.g. Guerrero-Rascado et al., 2008, 2009; Preißler et al., 2011; Cazorla et al., 2017; Córdoba-Jabonero et al., 2019; Fernández et al., 2019; López-Cayuela et al., 2023; Papanikolaou et al., 2024). Those results gain significance since desert dust aerosols affect Earth's energy balance. Thus, variations in the atmospheric dust loading may induce substantial changes in the radiative forcing of the climate system (Mahowald et al., 2010).

The aerosol radiative effect in the short-wave (SW) spectral range related to desert dust intrusions over the Iberian Peninsula has been widely investigated during the last years (e.g. Cachorro et al., 2008; Obregón et al., 2015; Sicard et al., 2016; Valenzuela et al., 2017; Granados-Muñoz et al., 2019; Córdoba-Jabonero et al., 2021a; Bazo et al., 2023, López-Cayuela et al., 2025). However, part of the literature often overlooked the aerosol radiative effects on the long-wave (LW) spectral range. This omission was primarily attributed to the intricate challenges associated with precisely quantifying the optical characteristics within this spectral domain (Roger et al., 2006; Mallet et al., 2008; Sicard et al., 2012). Moreover, the radiative forcing attributed to most aerosol categories (in particular, fine particles like pollution and smoke), is generally less pronounced in the LW range compared to their effects in the SW range. However, an exception arises with large and light-scattering particles (like mineral dust), which have been shown to possess a significant radiative forcing effect in the LW domain (e.g. Fouquart et al., 1987; di Sarra et al., 2011; Sicard et al., 2014a, 2022), highlighting again its climatic importance.

Recent research has demonstrated that separating the dust population into fine and coarse modes provides a useful framework to analyse size-dependent radiative effects. Sicard et al. (2014b) found that a clear distinction between dust modes is suitable for the reliable estimation of LW radiative forcing, particularly in the presence of large particles, which are common during mineral dust outbreaks, and mostly with intense dust incidence. This requirement is further supported by Adebiyi and Kok (2020), who discovered that the atmosphere burden of coarse dust is approximately four times larger than that simulated by current climate models. Consequently, an inadequate representation of coarse particles can lead to substantial errors in modelled dust–climate interactions. The findings also highlight the contrasting radiative effects associated with the two size modes, with coarse dust inducing a net warming at the top of the atmosphere (TOA) and fine dust contributing to cooling it, indicating their fundamentally different roles in the Earth's radiative budget. In addition, Kok et al. (2017) demonstrated that the dust direct radiative forcing is highly sensitive to the atmospheric dust size distribution, with climate models systematically underestimating the coarse-mode dust while overestimating the fine-mode fraction. Their analytical framework highlights the importance of considering the size-resolved dust mass and distribution, given that key radiative properties of dust (such as single-scattering albedo, asymmetry parameter, and extinction efficiency) are strongly dependent on particle size. It should be highlighted that the two dust modes differ not only in abundance and lifetime but also in their optical characteristics and radiative impacts. Thus, separating the dust population into fine and coarse modes provides a useful framework to analyse size-dependent radiative effects.

The present paper focuses on the assessment of the direct radiative effect (DRE) of dust particles in the LW range as well as their net effect. The event in study is an exceptionally intense and long-lasting Saharan dust event that crossed the Iberian Peninsula from 25 March–7 April 2021. This event was well characterized in López-Cayuela et al. (2023), which analysed the vertical behaviour of the optical and microphysical dust properties using polarization-sensitive lidar measurements from five Iberian lidar stations, separating the dust into its fine (Df) and coarse (Dc) components following the methodology of Mamouri and Ansmann (2014, 2017). The DRE analysis in the SW range can be found in López-Cayuela et al. (2025).

Although several works investigated the LW radiative effects associated with desert dust outbreaks over the Mediterranean basin (e.g. di Sarra et al., 2011; Perrone et al., 2012; Antón et al., 2014; Bazo et al., 2023), only a few studies have addressed the separation of both Df and Dc components (e.g. Sicard et al., 2014b, 2022). Two main conclusions highlight from these studies: (i) the quasi-linearity of LW radiative forcing at the bottom-of-atmosphere (BOA) and TOA with the Aerosol Optical Depth (AOD), and (ii) the high dependency of LW radiative forcing on the coarse-mode dust. Thus, the aim of the present study is to investigate whether this quasi-linearity holds for high AOD values (>0.50), in addition to distinguish the contribution of Df and Dc components to the DRE in the LW range. Indeed, this study introduces the novelty of simulating the LW dust DRE using two different approaches, used in López-Cayuela et al. (2025) for the SW range: (i) by simulating the contribution of Df and Dc components separately, and then estimating the total dust DRE as their sum (as $\text{DD}=\text{Df}+\text{Dc}$), and (ii) directly simulating DRE for the total dust component as a whole.

From a radiative transfer perspective, the most rigorous way to estimate DRE would be to integrate the optical properties over the full particle size distribution. In this work, however, the separation between Df and Dc is not introduced as a radiative-transfer requirement but arises naturally from the observational framework used. The polarization lidar methodology allows the retrieval of vertically resolved extinction profiles separately for fine and coarse dust particles (Tesche et al., 2009; Mamouri and Ansmann, 2014, 2017), as stated before. Treating these two modes independently therefore provides a direct way to propagate the observationally constrained dust components into the radiative transfer calculations. In addition, separating the Df and Dc modes offers a useful diagnostic framework for interpreting the size-dependent radiative effects of mineral dust. Df tends to dominate SW scattering, while Dc plays a more important role in LW emission and absorption, leading to potentially different or even compensating contributions to the DRE. Analysing both modes separately therefore allows for a better understanding of how different parts of the size distribution contribute to the overall radiative impact. Finally, this separated-mode approach also enables a direct comparison considering the total dust as a whole. In this study, the DRE obtained by summing the contributions from Df and Dc particles is compared with the DRE computed using the total dust population, allowing for assessing the sensitivity of the radiative effect to the size-mode resolved representation of dust.

The paper is organized as follows. The radiative transfer model and the parametrizations used in terms of the LW range are described in Sect. 2. The results and discussion are shown in Sect. 3. Finally, the main conclusions of this study are found in Sect. 4.

2 Methodology

2.1 Monitoring stations and lidar measurements

Lidar measurements were performed for vertical dust monitoring in five Iberian lidar stations: El Arenosillo/Huelva (ARN), Granada (GRA), Torrejón/Madrid (TRJ), and Barcelona (BCN) in Spain, and Évora (EVO) in Portugal. A more detailed description of those stations and lidar systems can be found in López-Cayuela et al. (2023). Briefly, polarized Micro-Pulse lidars operated at ARN, TRJ, and BCN sites. These single-wavelength elastic lidars (532 nm) operate continuously ($\mathrm{24}/\mathrm{7}$) with high pulse repetition frequency and low-energy laser emission, and include polarization capabilities (Campbell et al., 2002; Welton and Campbell, 2002; Flynn et al., 2007; Córdoba-Jabonero et al., 2018, 2021b; Welton et al., 2018). In addition, multi-wavelength Raman lidars were operative at EVO and GRA stations, which are part of EARLINET (European Aerosol Research Network; Pappalardo et al., 2014). These lidar systems use high-energy Nd:YAG lasers and provide elastic, Raman, and polarization-sensitive measurements at several wavelengths. In this work, only the common elastic 532 nm channel has been used. For the reader's convenience, information related to the location at the Iberian Peninsula and the period of the dust intrusion observed in each station can be found in Table S1 and Fig. S1 in the Supplement.

In order to obtain the Df and Dc contribution separately, the POLIPHON method (Polarisation Lidar photometer Networking method; Mamouri and Ansmann, 2014, 2017; Ansmann et al., 2019) was applied. In the first step of POLIPHON approach, background aerosols are separated from dust particles (Df+Dc). In the second step, Df and Dc are identified, and their specific backscatter coefficient profile is discriminated. The extinction coefficient (α⁵³²) for each component is then determined by considering the typical particle lidar ratio at 532 nm for each one (Ansmann et al., 2019). The uncertainties in the α⁵³² calculation by using this method are 30 %–50 %, 20 %–30 %, and 15 %–25 % for Df, Dc, and DD ($=\text{Dc}+\text{Df}$) dust, respectively (Ansmann et al., 2019).

2.2 Radiative transfer model: GAME. The MIE and LW modules

The GAME (Global Atmospheric Model; Dubuisson et al., 1996, 2004) has been used in increasing studies because of its significant advantage, i.e. the ability to fully represent the aerosol scattering and absorption in the LW region. Moreover, the model's moderate spectral resolution accounts for the spectral variations in aerosol properties, particularly in the infrared window. An extended description of the LW module of GAME can be found in Sicard et al. (2014a).

GAME calculates spectrally integrated upward and downward radiative fluxes in 40 plane and homogeneous layers from 0–100 km. Regarding the spectral limits, GAME employs 200–2500 cm⁻¹ (i.e. wavelength: 4.0–50.0 µm) with a fixed resolution of 20 cm⁻¹ (115 points). Moreover, this radiative transfer model considers thermal emission, absorption and scattering as well as their interplay employing the discrete ordinates method (DISORT, Stamnes et al., 1988). In the framework of GAME, an explicit account is taken for the absorption of gases, including H₂O, CO₂, O₃, N₂O, CO, CH₄, and N₂, using the correlated k-distribution as proposed by Lacis and Oinas (1991). Detailed insights into the computation of gas transmission functions can be found in Dubuisson et al. (2004) and Sicard et al. (2014a). The parameterization of gas absorption is based on pressure, temperature, and relative humidity profiles. Notably, these profiles are sourced from the Global Data Assimilation System (GDAS), provided by the National Oceanic and Atmospheric Administration (NOAA, last access: 28 March 2025).

The land surface temperature (LST) is a variable needed in the LW module. In this work, LST is provided by the Copernicus Land Service (https://land.copernicus.eu/global/products/lst, last access: 28 March 2025). Particularly, the hourly LST V2 dataset is used, which has uncertainties of less than 0.5 %. Moreover, the Earth's surface is assumed Lambertian, with a constant surface albedo of 0.017 in the LW spectral range. This value was determined by Sicard et al. (2014a) in Barcelona, based on the Clouds and Earth's Radiant Energy System (CERES) measurements in the spectral range of 8.1–11.8 µm, and averaged over the spring and summer seasons during five years. This same value is used at the five stations of this study, in the basis of the work of Zhou et al. (2013), which showed that the LW surface albedo remains relatively stable across the European continent.

Information on the aerosol shape, refractive index, size distribution, and density is required for an accurate calculation of their radiative properties. Yang et al. (2007) demonstrated that the non-sphericity effect of dust particles is negligible at thermal infrared wavelengths. Therefore, it is reasonable to assume that mineral dust is “spherical” in the LW range, and, hence, a Mie code can be applied for analysis. The spectral refractive index (both real and imaginary components) is identical to that reported in Sicard et al. (2014a) and was derived from measurements of long-range transported mineral dust collected in western Germany (Volz, 1983). The data that present the refractive index as a function of wavelength was obtained from Krekov (1993). The spectral variation of both the real and imaginary parts of the refractive index is illustrated in Fig. 1 of Sicard et al. (2014a). It should be highlighted that the refractive index used in the simulations (Volz, 1983), although assumed no varying, could be a source of uncertainty. Di Biagio et al. (2014, 2017) investigated the variability of the refractive index of mineral dust in LW as a function of its mineralogical composition and size distribution using in situ measurements. That study suggested that while a constant real refractive index can be probably assumed for dust from different sources, a varying complex refractive index should be used both at global and regional scale. They reported that for Saharan dust sampled at various sites, the real refractive index ranged from 1.3–2.0, and the complex refractive index ranged from 0.3–0.9 at a wavelength of 10 µm. The refractive index reported by Volz (1983), which has been used in the GAME simulations, is within those intervals of values for both the real and complex refractive index.

Figure 1 (a) Hourly land surface temperature (LST, in °C), where the red dots represent the cases coincident with lidar measurements; AERONET geometric median radius (r_g, in µm) and standard deviation (σ_g) for the (b) fine and (c) coarse modes, where the dashed lines represent the linear fitting of r_g over time; Episode-averaged values of (d) the Mie-derived normalized spectral extinction (${\mathit{\alpha }}_{\mathrm{LW}}\left(\text{Mie}\right)/{\mathit{\alpha }}_{\mathrm{532}}\left(\text{Mie}\right)$) (see Eq. 3), (e) asymmetry factor (g_LW), and (f) single scattering albedo (ω_LW), for the fine (blue), coarse (red) and total (yellow) modes. All the panels refer to El Arenosillo/Huelva (ARN) station; for the rest of stations, see the Supplement.

Moreover, the geometric median radius (r_g), and its standard deviation (σ_g), of the lognormal distribution are also needed in the Mie code. Those parameters are obtained for both the coarse and fine modes using column-integrated AERONET (Aerosol Robotic NETwork; http://aeronet.gsfc.nasa.gov, last access: 28 March 2025) Version 3 Level 2.0 data inversion products. AERONET provides the volume median radius (r_v) and its corresponding standard deviation (σ_v); hence, the following expressions were applied to determine both r_g and σ_g:

$$\begin{array}{}\text{(1)}& r_{\mathrm{g}}=r_{\mathrm{v}}{\mathrm{e}}^{-\mathrm{3}(\ln {\mathit{\sigma }}_{\mathrm{g}}{)}^{\mathrm{2}}},\end{array}$$

with σ_g=σ_v. These data were hourly averaged (and interpolated if missing). The column-integrated number concentration (N) is also derived. The AERONET column-integrated volume concentration (v_c), together with the r_v and σ_v, is used to calculate N as follows:

$$\begin{array}{}\text{(2)}& N=\mathrm{3}{\displaystyle \frac{v_{\mathrm{c}}}{\mathrm{4}\mathit{\pi }{r}_{\mathrm{v}}^{\mathrm{3}}\sqrt{\mathrm{2}\mathit{\pi }}\log {\mathit{\sigma }}_{\mathrm{v}}}}.\end{array}$$

The Mie module is capable of computing the spectral single scattering albedo (ω_LW), the asymmetry factor (g_LW) and the normalized extinction coefficient (${\mathit{\alpha }}_{\text{LW}}/{\mathit{\alpha }}_{\mathrm{532}}$) for each atmospheric layer. Then, the estimated extinction coefficient in the LW range, α_LW, distinguishing between Df and Dc modes, is calculated as follows:

$$\begin{array}{}\text{(3)}& {\mathit{\alpha }}_{\text{LW}}^i\left(\text{estimated}\right)={\mathit{\alpha }}_{\mathrm{532}}^i\left(\text{measured}\right){\displaystyle \frac{{\mathit{\alpha }}_{\text{LW, Mie}}^i}{{\mathit{\alpha }}_{\mathrm{532},\text{Mie}}^i}}\end{array}$$

where the upper-index i refers to total dust, Dc and Df, and α₅₃₂(measured) is the extinction coefficient at 532 nm as provided in López-Cayuela et al. (2023). Table 1 shows the input parameters used in the LW spectral range module as well as the data source.

Table 1 Input parameters and radiative properties for the GAME model in the LW spectral range. Note that z denotes the vertical dependence.

2.3 Dust radiative effect and heating rate estimation

The dust-induced DRE, simulated either at the BOA or the TOA, is defined as in López-Cayuela et al. (2025) (see there Eq. 1). In particular, the atmospheric DRE (DRE^ATM) is computed as the difference between the DRE at TOA (DRE^TOA) and that at BOA (DRE^BOA), that is, ${\text{DRE}}^{\text{ATM}}={\text{DRE}}^{\text{TOA}}-{\text{DRE}}^{\text{BOA}}$. In general, all those quantities are denoted as DRE${}_j^i$, where i stands for TOA and BOA, and j is the spectral band where DRE is calculated, i.e. j=LW, and net (SW+LW). All SW magnitudes were previously obtained in López-Cayuela et al. (2025). Both hourly DRE${}_j^i$ for SW and LW were computed for solar zenith angles (SZA)<90°, since GAME calculates those fluxes only during daytime.

As in López-Cayuela et al. (2025), both the hourly- and daily-averaged DRE${}_{\text{LW}}^i$ is calculated. In the SW range, the daily averages were computed as the mean (over 24 h) of the number of daytime hourly values, as SW fluxes during night-time are zero, unlike those in the LW range. Therefore, night-time hourly LW fluxes were assumed to be equal to the mean value of the daytime LW ones, and hence the daily DRE_LW was obtained from averaging those day-time and night-time-derived hourly (over 24 h) DRE_LW values. Similar procedure has been applied by other authors (di Sarra et al., 2011; Meloni et al., 2015; Sicard et al., 2022). Under this assumption, some uncertainty may still arise from two main factors affecting the DRE: differences in DOD values, and variations in the downward radiation flux during both day- and night-time. Regarding the first factor, the episode analysed was highly cloudy, resulting in numerous observational gaps during both day and night. However, recent studies by Tindan et al. (2023, 2025) investigating diurnal differences in dust aerosols across the dust belt have shown that such variations are insignificant over the Iberian Peninsula. With respect to the second factor, Granados-Muñoz et al. (2019) reported that the downward LW radiation flux displays a moderate daily variability (approximately 13 %), which would slightly modify the contrast in radiative forcing between daytime and nighttime conditions. Consequently, assuming night-time hourly LW fluxes to be equal to the mean daytime LW flux does not significantly affect the results of the present study.

Moreover, the fine-to-total (Df/DD) ratio (ftr) of the hourly-averaged DRE_j(ftr_DRE_j, j=LW, NET) is computed. In addition, a linear fitting analysis of this variable is performed over time, thus obtaining the slope of this linear fitting (δDRE_j), which serves as an indicator of the temporal rate of the relative contribution of Df particles to the DRE_j. The dust radiative efficiency (DREff_j) is also obtained from the slope of the linear fitting (forced to zero) of DRE values as a function of the dust optical depth at 532 nm (DOD⁵³²) along the event.

Following the same methodology as in López-Cayuela et al. (2025), differences in the dust-induced DRE (ΔDRE) as obtained from the two approaches are computed as follows:

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

where DRE${}_j^{\left(\mathrm{I}\right)}$ is the contribution to DRE of Df and Dc particles in each spectral range (i.e., j=LW, NET), that is,

$$\begin{array}{}\text{(5)}& {\text{DRE}}_j^{\left(\mathrm{I}\right)}={\text{DRE}}_j^{\text{DD}}={\text{DRE}}_j^{\text{Df}}+{\text{DRE}}_j^{\text{Dc}},\end{array}$$

and DRE${}_j^{\left(\text{II}\right)}$ is the contribution of the total dust as a whole, that is,

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

Moreover, the relative differences (Δ^relDRE) between the two approaches were calculated as:

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

As in López-Cayuela et al. (2025), the classical approach (i.e., without dust component separation) is adopted as the reference. Accordingly, throughout this manuscript, cases are described in which the component-separated DRE either overestimates or underestimates this classical approximation. Furthermore, a statistical analysis based on the relevant percentiles (P), e.g. P(25), P(50) (i.e. median), and P(75), of both ΔDRE and Δ^relDRE datasets was conducted to evaluate the significance of the discrepancies between the two methodologies.

Finally, it should be noted that aerosols predominantly exhibit a net cooling effect resulting from negative radiative forcing estimates, due to their inherent capacity to scatter solar radiation. However, certain aerosol types such as mineral dust are also able to absorb radiation to a greater or lesser degree, even likely leading to an opposite effect. Consequently, dust can induce heating in specific atmospheric layers, despite the potential net cooling effect observed for the overall atmospheric column (Pilewskie, 2007). The aerosol heating rate (AHR, K d⁻¹) is defined as the radiatively aerosol-induced rate of the temperature change in time ($\mathrm{\Delta }T/\mathrm{\Delta }t$) within a layer of the atmosphere. For a plane-parallel geometry, it can be expressed as follows:

$$\begin{array}{}\text{(8)}& \text{AHR}\left(z\right)={\displaystyle \frac{\mathrm{\Delta }T\left(z\right)}{\mathrm{\Delta }t}}=-{\displaystyle \frac{\textcolor{red}{}}{c_{\text{pd}}}}{\displaystyle \frac{\mathrm{\Delta }F\left(z\right)}{\mathrm{\Delta }p\left(z\right)}},\end{array}$$

where is the gravity acceleration (9.81 m s⁻¹), c_pd is the specific heat of dry air at constant pressure (p) (1005 kJ kg⁻¹ K⁻¹), Δp is the difference of the atmospheric pressure between two layers (Δp>0), and ΔF(z) represents the corresponding vertical difference in the flux (F(z); for Δz<0), which is defined as

$$\begin{array}{}\text{(9)}& F\left(z\right)=\left(F_{\mathrm{d}}^{\downarrow }\left(z\right)-F_{\mathrm{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⁻²) as computed by GAME, with and without dust presence, respectively. The arrows indicate whether the fluxes are downward (↓) or upward (↑).

3 Results and discussion

Comprehensive details describing the dust outbreak in overall, and regarding the methodology applied to derive dust optical and microphysical properties from polarized lidar measurements, are reported in López-Cayuela et al. (2023). The DRE analysis in the SW range can be found in López-Cayuela et al. (2025).

3.1 Dust radiative and microphysical properties in the LW range

Figure 1a shows the hourly LST at ARN station, for instance, during the dust outbreak period, where red dots represent the coincident values with lidar measurements when the DRE_LW can be calculated. Data from all the five Iberian lidar stations can be found in Fig. S2 in the Supplement. Consequently, compared to the SW study (López-Cayuela et al., 2025), 18 %–45 % fewer DRE_LW data were available for analysis. The lidar stations mostly affected by this lack of LST data were ARN and EVO. Regarding the LST results, except for a few days with missing data, the diurnal LST cycle is nicely visible at the stations. The maximum values ranged from approximately 28 °C (ARN, EVO, BCN) to 32 °C (GRA and TRJ) without significant changes over time (less than 0.02 °C). The maximum night/day difference ranged from approximately 18 °C (BCN) to 30 °C (TRJ).

In Fig. 1b and c, r_g and σ_g are represented as a function of time, and split into the fine and coarse modes, respectively, during the period for ARN station (Fig. S3 in the Supplement shows the same for the rest of the Iberian lidar stations). Those values are obtained from the AERONET r_v and σ_v (see Eq. 1). A linear fitting of those values over time is also performed. The episode-averaged value and the slope of the linear fitting (γ) are shown in Table 2. For the fine mode, the mean r_g (σ_g) value over the dust episode ranged from 0.076–0.093 µm (from 0.613–0.624 µm) at the southern stations (ARN, GRA and EVO). For TRJ and BCN, those values were lower, ranging from 0.059–0.067 µm (0.552–0.575 µm). It means that the fine particles were, on average, 10 %–30 % smaller at TRJ and BCN than at the southern stations. Regarding the γ values found, they were positive for GRA, TRJ and BCN, and negative for ARN and EVO. However, that increase/decrease over time was not significant, as it was less than 1 % µm d⁻¹ at every station. Therefore, by examining each station individually, the size of the fine particles did not vary considerably throughout the episode. These results are consistent with those obtained by Sicard et al. (2022) for BCN, showing a value of 0.7 % µm d⁻¹ during a summer Saharan dust outbreak in 2019.

Table 2 Episode-averaged median radius (r_g, µm) and standard deviation (σ_g, µm) at the five lidar stations: Barcelona (BCN), Torrejón/Madrid (TRJ), Évora (EVO), Granada (GRA) and El Arenosillo/Huelva (ARN) for the fine and coarse modes. The slope of each linear fitting (γ, % µm d⁻¹) and its standard error (in brackets) is also shown.

For the coarse mode, the lowest (highest) mean r_g over the dust episode was found at ARN (TRJ), showing a value of 0.471 (0.878). Regarding the σ_g, the lowest (highest) value was found at the same stations, that is 0.585 µm (0.653 µm). For the rest of the stations, r_g (σ_g) ranged from 0.529–0.584 µm (from 0.592–0.642 µm). Indeed, a value of 6 was found for the episode-averaged coarse-to-fine r_g ratio at the southern stations, being higher than 10 at TRJ and BCN. Similarly to the fine mode, that increase/decrease over time was not significant either (lower than 2 % µm d⁻¹) except for BCN, reaching almost 7 % µm d⁻¹. This increase in the geometric radius of the coarse particles at BCN has been observed previously where the coarse dust r_g increased during a summer Saharan dust outbreak in 2019 at a rate of $+\mathrm{9}\mathrm{\%}\mathrm{\mu }\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$ (Sicard et al., 2022). A possible explanation for this phenomenon was provided in the aforementioned study. Briefly, when the transport of mineral dust occurs over polluted regions with high humidity conditions, not only anthropogenic inorganic acids can be adsorbed onto the dust surface, forming hygroscopic salt compounds that coat the dust particles (Abdelkader et al., 2015; Athanasopoulou et al., 2016), but the formation of secondary pollutants is also enhanced (Querol et al., 2019; Xu et al., 2020).

Figure 1d–f show the episode-averaged optical properties introduced in the GAME Mie module (see Sect. 2.1), i.e. ${\mathit{\alpha }}_{\text{LW}}/{\mathit{\alpha }}_{\mathrm{532}},g_{\text{LW}}$ and ω_LW, for instance, at the ARN station. Note that the LW range stands for 4–50 µm in this work, and those quantities are separated into the fine, coarse and total modes for the five lidar stations. As observed, the most sensitive spectral window for radiative forcing is between 8 and 13 µm. For that reason, the analysis will be performed by averaging the properties in that spectral window, denoted here by δLW. All those quantities were derived separately for the fine and coarse dust contributions, as well as for the overall bimodal distribution (total dust). It should be mentioned that values for ${\mathit{\alpha }}_{\mathit{\delta }\text{LW}}^{\text{coarse}}/{\mathit{\alpha }}_{\mathrm{532}}^{\text{coarse}}$ are slightly higher than those for ${\mathit{\alpha }}_{\mathit{\delta }\text{LW}}^{\text{total}}/{\mathit{\alpha }}_{\mathrm{532}}^{\text{total}}$. Those ratios depicted in Fig. 1d are normalized against the corresponding α₅₃₂ for coarse (Dc) and total dust, respectively. Thus, the ${\mathit{\alpha }}_{\text{LW}}^{\text{coarse}}/{\mathit{\alpha }}_{\mathrm{532}}^{\text{coarse}}$ ratio is slightly higher due to a slightly smaller ${\mathit{\alpha }}_{\mathrm{532}}^{\text{coarse}}$ value as compared to ${\mathit{\alpha }}_{\mathrm{532}}^{\text{total}}$. Furthermore, due to the low sensitivity of fine mode in LW, ${\mathit{\alpha }}_{\mathit{\delta }\text{LW}}^{\text{fine}}/{\mathit{\alpha }}_{\mathrm{532}}^{\text{fine}}$ is an order of magnitude lower than${\mathit{\alpha }}_{\mathit{\delta }\text{LW}}^{\text{coarse}}/{\mathit{\alpha }}_{\mathrm{532}}^{\text{coarse}}$, with mean differences in the ratio of less than 6 % accounting for all the lidar stations on average along the dust episode.

Regarding the episode-averaged $g_{\mathit{\delta }\text{LW}}^{\text{fine}}$ and ${\mathit{\omega }}_{\mathit{\delta }\text{LW}}^{\text{fine}}$ at BCN station, they were from 2–10 times greater with respect to those values found in the other stations. However, noted that those values for the fine mode were remarkably lower than those for the coarse mode. Indeed, the episode-averaged fine-to-coarse ratio of g_δLW and ω_δLW is approximately 25 % and 5 %, respectively. Furthermore, the $g_{\mathit{\delta }\text{LW}}^{\text{coarse}}$ and ${\mathit{\omega }}_{\mathit{\delta }\text{LW}}^{\text{coarse}}$ properties were either rather similar or slightly higher than those corresponding to the total dust. That finding agrees with other studies (e.g. Sicard et al., 2014b). The coarse-to-total ratio was 1.0 and 1.2 on average for g_LW and ω_LW, respectively. However, in some cases, ${\mathit{\omega }}_{\text{LW}}^{\text{coarse}}$ was 1.7–3.0 times higher. The coarse-to-total ${\mathit{\alpha }}_{\text{LW}}/{\mathit{\alpha }}_{\mathrm{532}}$ ratio was 1.5–2.0 on average reaching up to values of 3.0–7.0 at several times during the dust event. This indicates that the separated coarse dust likely produces greater extinction in the LW range than that for the total dust (considering the bimodal distribution in overall). This hypothesis will be examined in Sect. 3.4, where the comparison between the two approaches considered in this study (see Sect. 2.2) will be addressed.

3.2 Long-wave dust direct radiative effect (DRE_LW)

A detailed description of the dust incidence of the Saharan intrusion by crossing the Iberian Peninsula is provided in López-Cayuela et al. (2023). In addition, the temporal evolution of the DOD⁵³² for the five lidar stations is shown in Fig. S4, where the particular days with high aerosol loads (i.e., hourly DOD⁵³²>0.5) are also indicated, occurring between 27 March and 1 April 2021 at various stations.

3.2.1 DRE_LW at BOA

Figure 2 shows the hourly dust direct radiative effect in the LW range (DRE_LW, W m⁻²) at BOA, TOA and ATM as induced by the Df and Dc particles at the ARN station, shown as an example. The rest of the stations can be found in Figs. S5–S8 in the Supplement along with the daily DRE_LW values for all stations. Table 3 shows the episode-averaged DRE_LW as induced by Dc, Df and DD at the five Iberian lidar stations, and at BOA and TOA. At all stations, DRE_LW is positive at BOA for both Df and Dc particles, indicating the expected dust-induced warming. During the most intense part of the episode (from 26 March–1 April for all the stations, except for BCN, where lasting until 3 April), hourly Dc (Df) DRE_LW values below +10 (+1) W m⁻² were mostly found. Note that 27 and 31 March stand out at ARN and TRJ (when DOD⁵³²>0.80) showing hourly DRE_LW values of ∼ +20 (∼ +2) W m⁻² for Dc (Df) particles. For the other stations, maximum hourly DRE_LW values are 40 %–60 % lower than those in ARN (see Table 3). During the rest of the dust event, Dc (Df) DRE_LW was lower than +5 (+1) W m⁻² at all stations. Regarding the daily DRE_LW (Fig. 2c and d), for days with daily DD DOD⁵³²≥0.50, DRE_LW values remained below +20 (+2) W m⁻² for Dc (Df) particles. For days with daily DD DOD⁵³²<0.50, DRE_LW values were below +10 (+0.8) W m⁻² for Dc (Df) particles.

Table 3 Episode-averaged dust direct radiative effect in the LW range (DRE_LW, in W m⁻²), and the standard deviation (SD, in brackets), at the BOA and TOA (and ATM) as induced by fine (Dc), coarse (Dc) and total dust (DD) at the five Iberian lidar stations. $\stackrel{\mathrm{‾}}{X}$ indicates the mean value for the whole event (SD are also shown), and X^max refers to the maximal value. The DREff (in $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$) denotes the DRE efficiency. The ftr_DRE denotes the hourly Df-to-DD DRE ratio (in %,), showing also the mean, median, minimum (min), and maximum (max) values; δDRE (in % d⁻¹) is the slope of the linear fitting analysis of the hourly ftr_DRE values along time.

Figure 2 Dust direct radiative effect in the long-wave range (DRE_LW, W m⁻²) at BOA (purple), TOA (green) and in the atmosphere (ATM, yellow) at the ARN station, for instance, as induced by: (a) fine dust (Df), and (b) coarse dust (Dc) particles. Daily mean values are also included at the top (as marked by the same colours). Daily mean DRE_LW values at (c) BOA and (d) TOA for Df (right) and Dc (left) particles at the five lidar stations.

Looking at the entire episode, the DRE_LW averaged values ($\stackrel{\mathrm{‾}}{{\text{DRE}}_{\text{LW}}}$) ranged from +2.6 to +6.5 W m⁻² for Dc particles, and from +0.2 to +0.4 W m⁻² for Df particles (Table 3). As noted in López-Cayuela et al. (2023), aerosol load exhibited very high variability across all stations during the study period, resulting in considerably large the standard deviation (SD; Table 3). The contribution of fine particles to DD DRE_LW is one order of magnitude lower than that of the Dc particles. The slope of the linear fitting of DRE_LW over time (δDRE_LW) showed values that ranged from −2.57 (BCN) to +0.65 (EVO) % d⁻¹ (Table 3). The mean ftr_DRE_LW values over the entire episode ranged from 6 % to 13 %, although the maximum hourly ftr_DRE_LW values varied between 10 % (GRA) and 41 % (BCN). Indeed, 10 %–15 % of the hourly ftr_DRE_LW values exceed the mean (±SD) values for the whole dust episode. These results indicate that Dc particles are the primary contributors to the DD DRE_LW. As expected, Df DREff_LW is also much smaller than Dc DREff_LW at all stations (Table 3). Specifically, Df DREff_LW is close to $+\mathrm{5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$, whereas Dc DREff_LW ranges from +24 to $+\mathrm{37}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$.

Performing a comparative analysis of the DRE_LW obtained in this work with previous studies can be challenging, as DRE_LW depends on multiple factors (e.g. fine-to-coarse r_g ratio, DOD, ftr_DOD, LST), and may vary significantly from one dust event to another. Moreover, DRE_LW is highly sensitive to the height reached by the dust layer (Dufresne et al., 2002; Sicard et al., 2022). Therefore, although the results of the present study are generally consistent with several previous works on mineral dust in the infrared range over Iberian Peninsula stations (e.g. Sicard et al., 2014b, 2022; Granados-Muñoz et al., 2019; Bazo et al., 2023), some differences can be expected. In particular, Sicard et al. (2022) studied the LW direct radiative effect during a summer dust outbreak in 2019 over BCN, using the same approach as in this study (separately considering the contribution of the Df and Dc modes). That work reported DREff_LW values of $+\mathrm{44.3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$ for Dc particles, and $+\mathrm{5.3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$ for Df particles. In the present study, similar Df DREff_LW values were obtained, whereas Dc DREff_LW were almost 50 % lower. This discrepancy likely arises from the strong sensitivity of LW radiative forcing simulations to aerosol load, coarse-mode particle radius, refractive index, vertical distribution, LST, and surface albedo (Sicard et al., 2014a). For example, substantially higher LST values were observed in Sicard et al. (2022), largely attributed to a concurrent heatwave during the Saharan dust outbreak, with night-time LST values greater than 15 °C and maximum daytime LST of 45 °C.

3.2.2 DRE_LW at TOA and in ATM

Similar to the BOA analysis, DRE_LW at TOA is positive, representing also a dust-induced heating (Fig. 2a and b). However, the magnitude is much lower (vs. DRE_LW at BOA). In particular, the maximum hourly DRE_LW values were 2–3 times lower at TOA than at BOA, and they were found at TRJ and BCN, where the dust plume reached higher altitudes (>6 km) (López-Cayuela et al., 2023). Those specific maximum DRE_LW values were +10.0 and +6.4 (+0.5 and +0.3) W m⁻², respectively, at TRJ and BCN for Dc (Df) particles, that is, between 2–5 times greater than for the rest of the stations (Table 3). The daily DRE_LW at TOA (Fig. 2c) shows values < +6.0 (+0.3) W m⁻² for Dc (Df) particles for high daily DD DOD⁵³² (≥ 0.50). For low and moderate daily DD DOD⁵³² (< 0.50), DRE_LW decreased to values lower than +2.0 (+0.1) W m⁻² for Dc (Df) particles. Regarding the mean DRE_LW values as averaged over the whole event, the Dc (Df) $\stackrel{\mathrm{‾}}{{\text{DRE}}_{\text{LW}}}$ values ranged from +0.9 to +2.3 (from +0.04 to +0.10) W m⁻² for all the stations (Table 3). As shown in Sect. 3.2.1, aerosol load exhibited a very high variability across all stations throughout the study period. Consequently, the SD is considerably high (Table 3).

The mean ftr_DRE_LW for the entire period varied approximately from +4 % to +8 % at all the stations except BCN, and slightly increased or decreased over time, depending of those stations (Table 3), but no significant impact was observed (δDRE_LW varied from around −0.5 to $+\mathrm{0.4}\mathrm{\%}{\mathrm{d}}^{-\mathrm{1}}$). However, at the BCN station, ftr_DRE_LW values of around +12 % with $\mathit{\delta }{\text{DRE}}_{\text{LW}}\sim -\mathrm{2.8}\mathrm{\%}{\mathrm{d}}^{-\mathrm{1}}$ were found. In addition, note that ftr_DRE_LW reached maximum values up to 41 %, with the Df contribution rather relevant. Indeed, between 8 % and 15 % of the hourly ftr_DRE_LW exceed the corresponding episode-averaged values (± standard deviation). Since similar results were found at BOA, this agrees with the findings reported by other authors supporting that DRE_LW is primarily dominated by the contribution of Dc particles (e.g. Sicard et al., 2022).

As expected, DREff_LW is much smaller at TOA than at BOA (Table 3). In particular, DREff_LW for Dc (Df) particles is overall nearly 4 (7) times lower at TOA (vs. at BOA) over the southern stations (ARN, GRA and EVO), and 2 (3) times lower at the rest. Additionally, the DREff_LW at TOA for the fine dust component is much smaller than for the coarse dust (Table 3). For all the stations, Df DREff_LW is lower than $+\mathrm{2}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$, and Dc DREff_LW ranges from $+\mathrm{5}-\mathrm{7}$ (at the southern stations) to $+\mathrm{12}-\mathrm{14}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$ (on TRJ and BCN). Thus, Dc DREff_LW is around 3–10 times greater than Df DREff_LW.

By comparing with other studies (e.g. Granados-Muñoz et al., 2019; Sicard et al., 2022), and regarding the southern stations, similar results for DRE_LW and DREff_LW are found. However, those parameters are 2–4 times greater at TRJ and BCN with respect to those previous studies. The difference in Dc DRE_LW may be explained by the finding of Sicard et al. (2014b), who reported that Dc DRE_LW exhibit little variations when the aerosol optical depth is kept constant. Indeed, Sicard et al. (2014b) and Dufresne et al. (2002) demonstrated that DRE_LW is highly dependent on the dust layer heights. In comparison with the heights reached by the dust intrusion as reported by Sicard et al. (2022), the observed differences in DRE_LW and DREff_LW could be based on this fact. Indeed, the dust plumes reached higher altitudes, especially at TRJ and BCN, during the dust outbreak examined in this study (>6 km height; López-Cayuela et al., 2023).

Finally, results on the dust radiative effect in the atmospheric column are reported in Table 3. For instance, hourly DRE_LW values in ATM at ARN station are shown in Fig. 2a and b for illustration. The rest of the stations can be found in Figs. S4–S8. By examining the DRE_LW at BOA and TOA, the DRE_LW at ATM is negative during the entire episode at all stations, as DRE_LW is lower at TOA than BOA, thus indicating a generalised atmospheric dust-induced cooling. The minimum hourly DRE_LW at ATM (i.e. the most negative) values are found at ARN and TRJ: −15.8 (−1.6) and −12.7 (−0.9) W m⁻² for Dc (Df) particles, respectively. For the rest of the stations the hourly DRE_LW at ATM minima (i.e. the most negative values) are 50 % lower (vs. ARN and TRJ stations; see Table 3). Regarding the episode-averaged estimates, atmospheric Dc (Df) DRE_LW ranged from −1.3 to −4.2 (from −0.1 to −0.3) W m⁻².

3.3 Dust net direct radiative effect (DRE_NET)

3.3.1 Relationship between DRE_LW and DRE_SW

It is known that the LW range is dominated by the Dc particles, whereas the Df particles induce a more pronounced effect in the SW range. In this work, ftr_DRE is less than 12 % in the LW range (Sect. 3.2), and 45 % in the SW range (López-Cayuela et al., 2025). This fact is illustrated in Fig. 3, where the DRE_LW with respect to DRE_SW ratio (DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$, in absolute value) is represented, giving an estimation of the percentage of radiative effect that the LW component represents compared to the SW one. For all the stations, the DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ for Df particles ranges on average 4 %–8 %, at both BOA and ATM. At TOA, the magnitude is lower, showing values of 1 %–4 %. Additionally, the Dc DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ ranges 39 %–54 % at BOA, and 20 %–50 % at both TOA and ATM. Particular mention should be made on the case of TRJ, where DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ for Dc reached values of 76 %. Other studies also report higher DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ values at BOA than at TOA. Particularly for desert dust outbreaks in the Mediterranean basin, daily DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ for total dust of 49 %–52 % and 26 %–35 % were found at BOA and TOA, respectively (di Sarra et al., 2011; Meloni et al., 2015). Sicard et al. (2022), which also performed the study by separating both Df and Dc components, found greater values of DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ for coarse dust (67 % at BOA, and 60 % at TOA). The reason could be attributed to multiple factors thus simulations of LW radiative forcing have demonstrated significant sensitivity to several key parameters, including aerosol load, coarse-mode particle radius, refractive index, vertical aerosol distribution, LST, and surface albedo (Sicard et al., 2014a). Therefore, discrepancies in this variable compared to other studies may be attributed to significant differences in the key parameters described above (see Sect. 3.2.1). Moreover, it is worthy to highlight, as Granados-Muñoz et al. (2019) pointed out, that the results at TOA might not be directly comparable to previous studies due to discrepancies in vertical resolutions within the GAME model for the SW and LW ranges above 4 km, potentially resulting in numerical artefacts in the derived outcomes.

Figure 3 Daily DRE${}_{\text{LW}}/{\text{DRE}}_{\mathrm{SW}}$ ratio (%; in absolute units) at the five Iberian lidar stations, at BOA (purple), TOA (green), and in ATM (yellow), corresponding to: (a) Df particles, and (b) Dc particles. The episode-averaged values are also shown on the right of each panel.

3.3.2 DRE_NET at BOA

By looking at the results, overall, DRE_NET is negative at BOA for all the stations, indicating a dust-induced net cooling effect. Figures 4a and b shows the hourly Df and Dc DRE_NET, respectively, at TOA, BOA and ATM in ARN, as an example. Results for the rest of stations are shown in Figs. S9–S12 in the Supplement. Moreover, the daily Df and Dc DRE_NET at BOA for all the stations considered in this study is shown in Fig. 4c. Table 4 shows the episode-averaged dust radiative effect in the net range ($\stackrel{\mathrm{‾}}{{\text{DRE}}_{\mathrm{NET}}}$, in W m⁻²) at BOA, TOA and ATM as induced by Df, Dc and DD at the five Iberian lidar stations.

Figure 4 The same as Fig. 2, but for the net DRE (DRE_NET).

Table 4 The same as Table 3, but for the episode-averaged dust net direct radiative effect (DRE_NET).

Dc (Df) $\stackrel{\mathrm{‾}}{{\text{DRE}}_{\mathrm{NET}}}$ ranges from −2.0 to $-\mathrm{6.0}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (from −2.9 to $-\mathrm{5.7}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). Note that those values are rather similar for Dc and Df particles, i.e. both dust components produce on average a similar net cooling at BOA. The daily DRE_NET (Fig. 4c) showed values from −13.1 to $-\mathrm{15.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (from −14.2 to $-\mathrm{20.6}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) for Dc (Df) particles during days with high daily DD DOD⁵³² (≥0.50). During days with moderate and low daily DD DOD⁵³² (<0.50), DRE_NET is always lower than $-\mathrm{4.8}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ ($-\mathrm{10.8}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) for Dc (Df). These results slightly agree with a few findings in the literature, reporting daily DD DRE_NET ranged from −14.6 to $-\mathrm{64.0}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Di Sarra et al., 2011; Meloni et al., 2015; Valenzuela et al., 2017). By definition, the DRE_NET is the sum of their SW and LW components. Therefore, those observed differences might be related to the varying balance between the DRE_SW (negative) and DRE_LW (positive). The maximum hourly DRE_NET was found at ARN station, showing values of −50.4 (Dc) and −43.0 (Df) W m⁻². For the rest of the stations, the maximum DRE_NET are lower (in absolute value) than those at ARN (20 %–60 % and 30 %–65 % for Dc and Df, respectively).

The impact of fine particles to DD DRE_NET is mainly due to their dominating contribution in the SW (vs. LW) range, as ftr_DRE_SW was estimated to be around 40 % for all the stations (López-Cayuela et al., 2025), meanwhile their LW contribution is between 6 % (TRJ) and 13 % (BCN) only (ftr_DRE_LW, see Table 3). Indeed, the ftr_DRE_NET values are 45 %–50 %, close to those obtained in the SW for all the stations (see Table 4 in López-Cayuela et al., 2025).

Moreover, DREff_NET values at BOA for the Df particles ranged from −128 to $-\mathrm{175}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$, which is approximately twice the Dc DREff_NET (see Table 4). In addition, the DD DREff_NET showed values from −78 to $-\mathrm{114}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$. Moreover, Granados-Muñoz et al. (2019) found values of DD DREff_NET approximately 1.5 times lower at GRA station than those reported in this work. Sicard et al. (2022) reported Dc and DD DREff_NET values approximately 2 and 1.5 times greater, respectively, at BCN station than those found in this work. Differences could be attributed to the radiative balance in DRE between the LW and SW ranges. In this work, the Dc and DD DREff_NET is reduced by a factor of 1.2 and 1.4, respectively, by counting on the LW contribution, with respect to the Dc and DD DREff_SW (López-Cayuela et al., 2025). Those reducing factors agree with the findings of Granados-Muñoz et al. (2019) and Sicard et al. (2022), that is, the DD DREff_NET is 1.1–1.6 times lower (vs. DD DREff_SW), although the Dc DREff_NET is reduced by a slightly higher factor of 2.5.

3.3.3 DRE_NET at TOA and in ATM

Similarly to the BOA, the DRE_NET is negative at TOA (Fig. 4d), indicating a dust-induced net cooling effect. In addition, DRE_NET values are 20 %–30 % lower (in absolute units), overall, representing a less pronounced net cooling at TOA with respect to that at BOA. Regarding DRE_NET on average for the entire episode ($\stackrel{\mathrm{‾}}{{\text{DRE}}_{\mathrm{NET}}}$, Table 4, Dc (Df) DRE_NET values range from −2.2 to $-\mathrm{4.9}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (from −2.4 to $-\mathrm{4.6}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). As it was for $\stackrel{\mathrm{‾}}{{\text{DRE}}_{\mathrm{NET}}}$ at BOA, note that those values are rather similar for Dc and Df particles, i.e. both dust components produce, on average, a similar net cooling at TOA.

The mean ftr_DRE_NET for the entire period was approximately ranging from 49 % to 58 % between stations (Table 4), with no significant temporal change observed ($\mathit{\delta }{\text{DRE}}_{\mathrm{NET}}=\mathrm{0.3}-\mathrm{0.4}\mathrm{\%}{\mathrm{d}}^{-\mathrm{1}}$) at the southern stations, meanwhile slightly higher δDRE_NET values of approximately 3–6 % d⁻¹ were found for TRJ and BCN.

The DD DREff_NET presented nearly 22 %–34 % smaller values at TOA (between −54 and $-\mathrm{75}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$) than those at BOA (see Table 3). In particular, DREff_NET for Dc particles ranged from −37 to $-\mathrm{56}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$, which are around half of the Df DREff_NET (i.e., from −91 to $-\mathrm{123}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$). It is important to note (as highlighted in Sect. 3.3.2) that the results obtained at TOA may not be directly comparable with those of previous studies due to differences in the vertical resolution of the GAME model in the SW and LW spectral ranges above 4 km, which could introduce numerical artefacts in results (Granados-Muñoz et al., 2019).

Table 5 Mean (SD), maximal (Max) and minimal (Min) values together the percentiles P(25), P(50) and P(75) of ΔDRE_LW (W m⁻²) and Δ^relDRE_LW (%) at BOA and TOA. SD stands for the standard deviation.

Finally, a dust-induced atmospheric net warming effect can be derived as DRE_NET is positive in ATM. Those results on the dust direct radiative effect in the atmospheric column are reported in Table 4. For illustration, the hourly DRE_NET values at ATM at ARN station are shown for Df and Dc particles, respectively, in Fig. 4a and b (results for the rest of stations are shown in Figs. S9–S12). The maximum hourly DRE_NET values at ATM are found at ARN and TRJ, showing values of +19.4 (+14.0) and +19.2 (+11.6) W m⁻² for Dc (Df) particles, respectively. For the rest of the stations, the hourly DRE_NET maxima at ATM ranged from +5.0 to +8.4 W m⁻² and from +8.3 to +9.2 W m⁻² for Df and Dc particles, respectively (see Table 4).

3.4 Differences in DRE_LW and DRE_NET as estimated using different approaches

Following the approach applied by López-Cayuela et al. (2025) for the SW range, the differences in DRE_LW and DRE_NET at all the stations were examined using the two approaches described in Sect. 2.2.

Relative differences in the LW range (Δ^relDRE_LW; see Eq. 7) with respect to the classical approach (DRE${}_{\text{LW}}^{\left(\text{II}\right)}$, see Eq. 5) are shown in Fig. 5 as a function of SZA, highlighting the dependence on DD DOD⁵³². The entire dataset was considered, covering the period from 25 March to 7 April 2021 at all five Iberian lidar stations. As discussed in López-Cayuela et al. (2025), the significant Δ^relDRE_SW values found for SZA>70° are attributed to the intrinsic uncertainty in GAME simulations arising from the assumption of a plane-parallel atmosphere, and, hence, these values should be discarded. However, no clear correlation was observed between Δ^relDRE_LW and SZA. At BOA (TOA), mean Δ^relDRE_LW values of approximately +8.5 % (+6.5 %) were obtained, although relatively large SD values were observed (∼ 25 %–27 %, see Table 5). Indeed, comparable Δ^relDRE_LW values are found for SZA<70° (see Table 5). Moreover, no clear relationship is evident between Δ^relDRE_LW and DD DOD⁵³². An analysis of percentiles further revealed consistent patterns at both levels: P(75) around +16 %, P(50) in the interval of +0.8 % to +1.2 %, and P(25) close to −10 % independently on SZA (see Table 5 and Fig. 5). These results indicate that larger absolute DRE${}_{\text{LW}}^{\left(\mathrm{I}\right)}$ values relative to DRE${}_{\text{LW}}^{\left(\text{II}\right)}$ are predominantly derived when the full dataset is considered. Specifically, 75 % of the Δ^relDRE_LW values are higher than around −10 %, with only 25 % falling between −10 % to +1 %.

Figure 5 Relative differences in DRE_LW (Δ^relDRE_LW, in %) as obtained from the two approaches (Eq. 7) as a function of SZA at: (a) TOA, and (b) BOA, for all five lidar stations from 25 March–7 April 2021. The dependence on DD DOD⁵³² is shown as a colour-scaled bar at the top.

Regarding the absolute differences in DRE_LW (ΔDRE_LW; see Eq. 4), those computed from the full dataset were found to be approximately 3–4 times larger at BOA than TOA, with mean (SD) values of +0.3 (1.3) and +0.1 (0.5) W m⁻², respectively, which are close to zero. Maximum (minimum) ΔDRE_LW values of +9.7 (−2.1) and +2.6 (−1.2) W m⁻² were reached at BOA and TOA, respectively.

However, when the dependence of ΔDRE_LW on r_g is examined for the fine and coarse dust, a differentiated behaviour can be observed. Figure 6 displays ΔDRE_LW as a function of DRE${}_{\text{LW}}^{\left(\text{II}\right)}$ at both BOA and TOA, highlighting the dependence on fine r_g. It was found that, as fine r_g increases ΔDRE_LW shifts from negative to positive values; the same behaviour is observed depending on the coarse r_g. The inflexion point (ΔDRE_LW=0) was estimated at thresholds of approximately 0.1 µm for fine r_g (or 0.7 µm for coarse r_g; data not shown). As this study focuses on the relevance of fine particles, reference will be made to the threshold related to fine r_g throughout this section.

Figure 6 Differences in DRE_LW as obtained from the two approaches ($\mathrm{\Delta }{\text{DRE}}_{\text{LW}}={\text{DRE}}_{\text{LW}}^{\left(\mathrm{I}\right)}-{\text{DRE}}_{\text{LW}}^{\left(\text{II}\right)}$, Eq. 4; W m⁻²) as a function of DRE${}_{\text{LW}}^{\left(\text{II}\right)}$ at: (a) TOA, and (b) BOA, for all five lidar stations from 25 March–7 April 2021. The dependence on the fine r_g (µm) is shown as a colour-scaled bar at the top.

For cases with fine r_g<0.1 µm (i.e. for rather small fine dust particles), the use of the dust-mode separation approach resulted in negative ΔDRE_LW at both BOA and TOA (see Fig. 6). This reveals an underestimation of DRE_LW values for separated dust components, leading to a less pronounced warming effect. Conversely, when fine r_g≥0.1 µm, ΔDRE_LW tended to be positive, resulting in an overestimation of DRE_LW values with respect to the traditional method, and hence in a more pronounced dust-induced warming effect. Those results are aligned with Sicard et al. (2014b), who demonstrated that the radiative forcing produced by aerosols whose size distribution is dominated by the coarse mode is higher than the estimated by the classical approach. In terms of mean values, the largest differences are found for size distributions dominated by finer particles, for which ΔDRE_LW exhibited mean (SD) values of −0.04 (0.58) and −0.03 (0.22) W m⁻² at BOA and TOA, respectively. In contrast, for cases with fine r_g ≥ 0.1 µm, ΔDRE_LW presented mean (SD) values of +3.1 (2.5) and +0.8 (0.8) W m⁻² at BOA and TOA, respectively (see Table 5).

By looking at the main percentiles P(75), P(50) and P(25), as computed from the statistical analysis of ΔDRE_LW (see Table 5), the data distribution is nearly symmetrical, with median values closely matching those mean ones for both fine r_g intervals. This same pattern depending on fine r_g is observed at both BOA and TOA, though finding lower ΔDRE_LW at TOA. Indeed, P(25) values indicate that 75 % of the ΔDRE_LW values are above +2.0 and +0.6 W m⁻² at BOA and TOA, respectively, for cases with fine r_g < 0.1 µm, but close to zero for the remaining cases. These discrepancies observed in dependence of the size interval at both BOA and TOA further emphasizes the critical role of particle size in modulating the vertical distribution and net effect of dust radiative forcing.

The same analysis has been performed for the differences in DRE_NET. Figure 7 shows Δ^relDRE_NET as a function of SZA for all five lidar stations and the whole dataset. In this case, a clear dependence on SZA is observed, originating from the same effect as in the SW range (López-Cayuela et al., 2025). Specifically, larger differences are found for SZA>70°, although the effect is less pronounced than in the SW range (López-Cayuela et al., 2025), as it is modulated by the contribution of the LW range to the net radiative balance. At TOA and for SZA>70°, Δ^relDRE_NET values are mostly positive, ranging from approximately −10 % to +65 %, with a mean (SD) value of +14.0 (20.0) % (see Fig. 7a). For the same SZA range, Δ^relDRE_NET at BOA showed values that ranged from around −35 % to +90 %, and with a mean (SD) value of +12.7 (22.7) % (see Fig. 7b). As explained in López-Cayuela et al. (2025), the significant Δ^relDRE_NET found for SZA>70° are associated to the intrinsic uncertainty in GAME simulations resulting from the model assumption of a plane-parallel atmosphere, and hence those values should be discarded.

Figure 7 The same as Fig. 5, but for DRE_NET (Δ^relDRE_NET, in %). The vertical dashed black line denotes SZA=70°.

Thus, once disregarding values of Δ^relDRE_NET for SZA>70°, Δ^relDRE_NET are mostly negative at both BOA and TOA, showing mean (SD) values of −4.8 (6.6) % and −8.5 (13.0) %, respectively. This is also corroborated by looking at the percentiles P(25) and P(75), which show values, respectively, of −9.6 % and −1.6 % at BOA, and −13.4 % and −1.2 % at TOA (see Table 6). Indeed, those results indicate that 75 % of Δ^relDRE_NET values are below around −2 % and −1 % at BOA and TOA, respectively, with minima of −18.2 % and −80 %. This represents, as DRE_NET is negative, a less pronounced net cooling at both BOA and TOA when the Df and Dc contribution is separately (vs. total dust) accounted for and showing larger differences at TOA (vs. BOA).

Table 6 The same as Table 5, but for ΔDRE_NET (W m⁻²) and Δ^relDRE_NET (%).

Finally, Fig. 8 shows the differences in DRE_NET (ΔDRE_NET) obtained from the two approaches at both BOA and TOA for all five lidar stations. It should be noted that absolute ΔDRE_NET tend to increase as DOD⁵³² increases. In general, ΔDRE_NET were mostly close to zero at lower DOD (<0.2), and increased somewhat at moderate/high-dust-load conditions (DOD > 0.50). Indeed, those differences in DRE_NET reached minimum (maximum) values of −6.4 (+6.4) and −10.4 (+2.3) W m⁻², respectively, at both BOA and TOA for SZA < 70°, showing positive (and close to zero) mean (SD) values of +0.5 (1.5) and +0.2 (1.4) W m⁻². This can be corroborated by looking at the percentiles: 75 % of ΔDRE_NET values are mostly positive (i.e., $\mathrm{P}\left(\mathrm{25}\right)=+\mathrm{0.2}$ and +0.1 W m⁻² at BOA and TOA, respectively). Table 6 shows all those values. Moreover, it should be noted that a differentiated behaviour is observed around a DRE_NET threshold of $-\mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. In particular, when DRE${}_{\mathrm{NET}}^{\left(\text{II}\right)}>-\mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, similar low mean ΔDRE_NET values are obtained at BOA ($+\mathrm{0.5}\pm \mathrm{2.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) and TOA ($+\mathrm{0.5}\pm \mathrm{0.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). However, for DRE${}_{\mathrm{NET}}^{\left(\text{II}\right)}\le -\mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, corresponding to higher dust load conditions, ΔDRE_NET show positive mean values at BOA ($+\mathrm{0.6}\pm \mathrm{0.8}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) and negative mean values at TOA ($-\mathrm{2.0}\pm \mathrm{2.6}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$). Moreover, as shown in Table 6, 75 % of ΔDRE_NET are above $+\mathrm{0.3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ at BOA, and below $-\mathrm{0.3}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ at TOA. Overall, these results would indicate a less pronounced net cooling at BOA in contrast of a more pronounced net cooling at TOA when the separation contribution of the Df and Dc particles (vs. total dust) is regarded under high dusty conditions. This highlights a potential modulation of the dust impact in the atmosphere, which could potentially be able to produce an atmospheric net cooling (contrary to what stated in Sect. 3.3). However, those final remarks should be carefully regarded as only 8 % of those examined DRE_NET profiles correspond to DOD⁵³² values greater than 0.5.

Figure 8 Differences in DRE_NET as obtained from the two approaches ($\mathrm{\Delta }{\text{DRE}}_{\mathrm{NET}}={\text{DRE}}_{\mathrm{NET}}^{\left(\mathrm{I}\right)}-{\text{DRE}}_{\mathrm{NET}}^{\left(\text{II}\right)}$; Eq. 4; W m⁻²) at (a) TOA, and (b) BOA, for all five lidar stations from 25 March–7 April 2021. The dependence on the DD DOD⁵³² is shown as a colour-scaled bar at the top. Data for SZA>70° are highlighted by cross symbols. The vertical dashed line indicates DRE${}_{\mathrm{NET}}^{\left(\text{II}\right)}=-\mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$.

3.5 Aerosol heating rate

The vertical aerosol heating rate (AHR) has been computed in the SW and LW range (see Sect. 2.2, Eq. 8) for all the dust components (DD, Df, Dc). Maxima of the hourly AHR (AHR^max, K d⁻¹) during the entire dust episode at each lidar station, along with the episode-averaged values and their corresponding altitudes, are shown in Table 7. In the SW range, the AHR_SW is predominantly positive, with maximum values within the dust layer, indicating a warming effect in the atmosphere. On the contrary, near the surface, AHR_SW are mostly negative (cooling effect). Since the fine-to-total AHR ratio in the SW (ftr_AHR_SW) remains nearly constant for all stations (around 30 % within the dust layer; see Table 7), the discussion primarily focused on DD AHR_SW (for clarity, Fig. S13 in the Supplement shows the AHR_SW at the five Iberian lidar stations). As stated in several works, the AHR_SW is linked to the vertical distribution of the dust extinction, and its magnitude increases with the DOD (Perrone et al., 2012; Meloni et al., 2015; Peris-Ferrús et al., 2017). An extensive study of the vertical dust extinction distribution can be found in López-Cayuela et al. (2023, 2025).

Table 7 Maxima of the hourly aerosol heating rates (${\text{AHR}}_i^{max}$, K d⁻¹) found for the entire episode at each lidar stations. Values for the fine (Dc), coarse (Dc) and total dust (DD) are shown. The episode-averaged of those ${\text{AHR}}_i^{max}$ ($\stackrel{\mathrm{‾}}{{\text{AHR}}_i}$, K d⁻¹) and their corresponding heights ($\stackrel{\mathrm{‾}}{{\text{z\_AHR}}_i}$, km), and the fine-to-total AHR ratio (ftr_AHR_i) is also shown. The index i stands for the SW, LW and NET ranges. The standard deviation (SD) is shown in brackets.

To summarize, the dust plume initially appeared below 3 km at the southern stations (ARN, GRA, EVO) on 25–26 March. On 27 March, enhanced atmospheric instability lifted the plume up to 6 km. The strongest intrusion occurred on 29–31 March, with dust extending from the surface to ∼ 7 km height. From 1 April, the plume weakened and descended to ∼ 3 km (see Fig. 2 in López-Cayuela et al., 2025). Regarding the AHR_SW, the maxima varied from $\sim \mathrm{0.3}-\mathrm{1.0}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ at the beginning of the episode, peaked $\sim \mathrm{3}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ on 31 March, and decreased to $\sim \mathrm{0.1}-\mathrm{0.3}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ at the end of the episode, at altitudes of 3–5 km. At the central station (TRJ), the dust plume was initially detected below 4 km on 26 March, ascending to 10 km later that day. In the following days, plume top heights fluctuated between 6–8 km height, occasionally reaching 10 km. The strongest intrusion also occurred on 29–31 March. From 1 April, the plume subsided from 8–4 km with reduced intensity (see Fig. 2 in López-Cayuela et al., 2025). The maxima AHR_SW varied at altitudes of 4–6 km from $\sim \mathrm{0.4}-\mathrm{0.7}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ at the beginning of the dust outbreak, peaking $\sim \mathrm{2}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ on 31 March, to $\sim \mathrm{0.1}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ at the end of the episode. At BCN, the dust plume was persistently stratified and less intense than at the other stations, although plume tops occasionally reached 10 km. On 28 March, dust was confined to 2–3 km. On 29 March, a more complex structure formed with two distinct layers (at 2–3 and 9–10 km) in the morning and three layers (at 1–2, 4–7, 8–10 km) later in the day. Similar stratification persisted in the following days, peaking on 1 April. Afterwards, the dust intrusion gradually weakened until 6 April (see Fig. 2 in López-Cayuela et al., 2025). The maxima AHR_SW varied at altitudes of 2–6 km, from $\sim \mathrm{0.1}-\mathrm{0.3}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ until 1 April, when peaked to $\sim \mathrm{0.7}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$, to finally decreased to $\sim \mathrm{0.1}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ at the end of the episode. Averaging over the entire dusty period, the maximum hourly DD AHR_SW values peaked at mean altitudes between 2.0 (BCN) and 5.3 (GRA) km, with magnitudes ranging from +0.20 to $+\mathrm{0.50}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ (see Table 7).

Regarding the LW range, AHR_LW is predominantly negative. As the fine mode contribution to the DD AHR_LW is low (ftr_AHR_LW<16 %; see Table 7), the discussion will focus on DD AHR_LW (the vertical distribution of DD AHR_LW can be found on Fig. S14 in the Supplement). The absorption of SW radiation by the dust layer led to the emission of LW radiation in all directions, resulting in a negative AHR_LW (cooling effect). Indeed, examining the vertical AHR structure, the AHR_LW profiling usually peaks at lower altitudes than those for AHR_SW, with maxima located below the dust layer. This behaviour was also found in previous studies (e.g. Sicard et al., 2014a). As stated in Meloni et al. (2015), when the dust intrusion is structured in multiple layers, the sign of AHR_LW can switch from negative to positive below the most dust-loaded layer, depending on the absorption of the lowermost layer. Generally, when the extinction coefficient indicates a significant aerosol load in the lowermost atmosphere, AHR_LW remains predominantly negative, reflecting to thermal emissions within the layer itself. Conversely, AHR_LW tends to be positive when the aerosol extinction near the surface is negligible, as the primary dust layer induces a LW heating effect beneath it, attributable to the absorption of local thermal radiation emitted by the dust layer (Meloni et al., 2015). The maximum (negative) values of the hourly DD AHR_LW occurred on days of strongest dust incidence, at altitudes below the main dust layer (López-Cayuela et al., 2023), ranging from −0.12 (GRA) to −0.85 (TRJ) K d⁻¹. The period-averaged values of the maximum (negative) hourly DD AHR_LW for each station and spectral range are shown in Table 7, ranging from −0.12 K d⁻¹ at 1.6 km in TRJ to −0.04 K d⁻¹ at 1.1–1.7 km for the rest of stations.

Note that AHR results for both SW and LW ranges are consistent with previous studies on mineral dust, which reported hourly AHR_SW and AHR_LW values ranging from +0.30 to $+\mathrm{3.80}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$, and from −0.30 to $-\mathrm{0.70}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$, respectively (Sicard et al., 2014a; Meloni et al., 2015; Peris-Ferrús et al., 2017; Valenzuela et al., 2017; Bazo et al., 2023).

Figure 9 (Left) Vertical distribution of the net aerosol heating rates (AHR_NET, K d⁻¹) corresponding to dust (DD) particles at the five Iberian lidar stations (from North-East to South-West, by decreasing latitude): (a) BCN, (b) TRJ, (c) EVO, (d) GRA and (e) ARN stations. (Right) An example of an hourly-averaged AHR profile for the SW (red) and LW (blue) range at each station. These specific profiles are marked on the left panel between white arrows. Solid, dashed and dotted lines refer to DD, Df and Dc AHR_NET, respectively.

Figure 9 shows the vertical AHR_NET (K d⁻¹) for DD particles at the five Iberian lidar stations along the dust event, together with examples of hourly AHR_SW and AHR_LW profiling at each site. As stated for the AHR in the SW range, ftr_AHR_NET is also nearly constant across all stations (nearly 30 % within the dust layer; see Table 7). AHR_NET profiles indicate negative values (net cooling) in the lowermost atmosphere during the dusty period, as both LW and SW contributions are negative, which is consistent with the negative DRE_NET at BOA (Sect. 3.3.1). In contrast, positive AHR_NET (net warming) dominates within the dust layer, where the (positive) AHR_SW prevails. Figure S15 in the Supplement shows the vertical distribution of the LW-to-net AHR ratio (${\text{AHR}}_{\text{LW}}/{\text{AHR}}_{\mathrm{NET}}$; %). It can be seen that the LW contribution to AHR_NET is generally <10 % inside the dust layer, but nearly all of it occurs below, and to a lesser extent above, the layer (see Fig. S15). Overall, AHR_NET is positive in the most dust-loaded layer and negative below and above it, typically between 2–4 km height and occasionally above 8 km, which is consistent with López-Cayuela et al. (2023). Although AHR_SW dominates, AHR_LW remains relevant as it modulates the strength of the net effect. As expected, the maxima values are found on day of maxima DOD, within the dust layer, and reaching values from $+\mathrm{1.83}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ (ARN) to +0.59 K d⁻¹ (EVO). The dust period-averaged DD AHR_NET showed values ranging from $+\mathrm{0.22}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ around 5 km (BCN) to $+\mathrm{0.44}\mathrm{K}{\mathrm{d}}^{-\mathrm{1}}$ around 3 km (TRJ) (see Table 7).

4 Summary and conclusions

This work is complementary to the research conducted by López-Cayuela et al. (2023, 2025), with the aim of introducing thus the closure study about the vertical radiative impact of an intense and long-lasting Saharan dust outbreak over the Iberian Peninsula in springtime 2021. In this work, the temporal variation of the DRE in the LW range and net DRE was estimated, separating the Df and Dc contributions. For that purpose, lidar observations in five Iberian stations (namely ARN, GRA, TRJ, and BCN in Spain, and EVO in Portugal) were used. The key findings are summarised below.

The availability of DRE_LW computations was reduced by 18 %–45 % compared to the DRE_SW reported in the accompanying article due mainly to limited LST measurements, particularly at ARN and EVO. Despite some data gaps, the diurnal LST cycle was clearly observed, with maximum values between 28 and 32 °C and minor temporal variability (<0.02 °C). The fine geometric median radius r_g and its standard deviation σ_g were larger at the southern Iberian Peninsula stations (ARN, GRA, EVO) than at TRJ (central Iberian Peninsula) and BCN (north-eastern Iberian Peninsula), indicating 10 %–30 % smaller fine particles in the latter stations. Temporal trends in r_g were negligible across stations ($<\mathrm{1}\mathrm{\%}\mathrm{\mu }\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$), indicating that fine particle size remained stable during the dust outbreak. Although many stations showed stable coarse particle sizes, BCN experienced a significant increase ($\sim \mathrm{7}\mathrm{\%}\mathrm{\mu }\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$), which is consistent with previous studies. This may be due to dust interaction with anthropogenic pollutants and humid conditions enhancing particle growth. The coarse mode dominated the LW extinction and scattering processes, as indicated by higher g and ω values, and their coarse-to-total ratios exceeding unity, reaching up to 3–7 for the extinction. The dominant role of coarse dust in the LW spectral range (particularly 8–13 µm) implies that LW radiative forcing estimations should be mainly affected by coarse-mode contributions. These findings align with previous literature and were validated in Sect. 3.1.

On the one hand, a dust-induced warming at BOA was consistently observed for both fine and coarse dust particles, with Dc contributing the most. During the highest incidence of the dust episode, maximum hourly DRE_LW values at BOA reached up to $+\mathrm{20}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Dc) and $+\mathrm{2}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ (Df), especially at ARN and TRJ stations. The daily DRE_LW values were significantly lower when the dust optical depth at 532 nm (DOD⁵³²) was below 0.50, indicating a direct dependency on dust concentration. Episode-averaged DRE_LW values ranged from +2.6 to $+\mathrm{6.5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for Dc and from +0.2 to $+\mathrm{0.4}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for Df particles. This shows that Df contributes an order of magnitude less to the LW radiative forcing. The relative contribution of Df (ftr_DRE_LW) was 6 %–13 % on average, but it could reach up to 41 % during some times, particularly at BCN. Despite these peaks, the Dc component remains the dominant driver of the total dust DRE_LW. Regarding the radiative efficiency, DREff_LW values were much higher for Dc (from +24 to $+\mathrm{37}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$) than Df ($\sim +\mathrm{5}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathit{\tau }}^{-\mathrm{1}}$), reinforcing the greater radiative impact of coarse dust.

This study also provides a comprehensive assessment of the DRE_NET, which is consistently negative at both BOA and TOA, thereby reflecting a net cooling effect induced by dust. In contrast, a positive DRE_NET was observed within the atmospheric column, suggesting a dust-driven net warming at this level. Regarding the effect of the Df particles at both BOA and TOA, they contributed a maximum of 12 % and 30 % to the DRE_NET in the LW and SW ranges, respectively. The Df DRE_NET was similar to the Df DRE_SW, since the Df DRE_LW is nearly negligible. Indeed, ftr_DRE_NET was approximately 45 %–50 % at BOA, and 50 %–60 % at TOA. In addition, the DREff_NET corresponding to Df particles was almost half of Dc DREff_NET at both BOA and TOA. Among all stations, ARN exhibited the highest hourly DRE_NET magnitudes, highlighting the influence of local atmospheric and surface conditions on radiative forcing. The inclusion of the LW component was found to decrease the net radiative efficiency, with reduction factors ranging from 1.2–2.5, depending on the dust mode fraction. This underscores the importance of accounting for the SW–LW balance when quantifying the net radiative impact of dust. All these findings contribute to a better understanding of the vertical distribution of dust radiative effects and their implications for regional climate over the Iberian Peninsula.

Concerning the vertical distribution of the aerosol heating rates, the peak of AHR_SW profiling occurred at a higher altitude than that of the AHR_LW one. Moreover, the AHR_SW was predominantly positive (warming effect within the dust layer) unlike the AHR_LW, which was negative (cooling effect). Hence, in the dust layer, AHR_NET displayed a warming effect. On the contrary, below (until 2–4 km) and above (beyond 8 km) the dust layer, AHR_NET indicates a cooling effect. Moreover, the ftr_AHR_NET in the dust layer at all the stations is on average nearly 30 %, with the contribution of the AHR in the LW range reaching a maximum of 10 % to AHR_NET. Despite AHR_SW dominates leading to a predominant atmospheric warming effect of dust, the relevance of the opposite (cooling) AHR_LW effect relies on the potential atmospheric modulation by reducing the strength of the warming net impact.

On the other hand, as a novelty in this work, two methodologies for estimating DRE in each spectral range were examined in detail. For that purpose, differences in DRE (ΔDRE) were analysed by comparing the contribution of the two dust components computed separately (Dc, Df; approach 1) with the classical estimation (approach 2), which considers total dust as a whole. Results revealed that the classical approach underestimated the DRE_LW, with mean (SD) relative differences of approximately +9 % (25 %) and 8 % (26 %) at both BOA and TOA, respectively. Moreover, for cases with rather small fine dust particles, the use of dust-mode separation approach resulted in negative LW differences at both BOA and TOA, revealing an underestimation of DRE_LW values when using the separated dust components approach. This leads to a less pronounced warming effect unlike when fine dust radii exceed 0.1 µm, resulting in a DRE_LW overestimation with respect to the traditional method and, consequently, to a more pronounced dust-induced warming effect. These discrepancies observed in dependence of the size interval at both BOA and TOA further emphasizes the critical role of particle size in modulating the vertical distribution of the DRE_LW and then affecting the net dust radiative forcing.

Indeed, the dust-induced net effect is strongly affected by the SW range but modulated by the LW range. On average, DRE_NET is overestimated by the classical approach, with mean relative differences (SD) of −5 % (7 %) at BOA and −9 % (13 %) at TOA. In addition, under moderate-to-high dust conditions, the net cooling at BOA is less pronounced, whereas is stronger at TOA, when the separation contribution of the Df and Dc particles (compared to total dust) are considered. These results highlight a potential modulation of the dust impact in the atmosphere, which could actually lead to net atmospheric cooling.

The literature consistently indicates that global models underestimate the burden and extent of Dc particles, leading to an underestimation of their LW warming and reduced SW cooling effectiveness. Incorporating dust observations alters both regional and global direct radiative forcing, improving agreement between simulations and observations. While SW dominates surface cooling, Dc-induced LW contributes to warming at TOA and within the atmospheric column, further modulating semi-direct cloud effects. Dust size-resolving studies show that computing DRE in bulk overestimates (in absolute terms) SW cooling; when fine and coarse dust components are treated separately, Df and Dc contributions partly offset each other, yielding a less negative radiative forcing (i.e., less pronounced cooling) at TOA, whereas the net DRE remains strongly negative at BOA. Recent estimates highlight a globally significant positive LW direct forcing during the industrial era, previously overlooked due to the lack of realistic coarse-mode representation. Therefore, separating fine and coarse dust contributions to DRE computation is crucial, affecting the TOA radiative forcing sign, BOA attenuation magnitude, atmospheric heating profiles, and model biases linked to underrepresented Dc particles in their LW effects. Thus, the dust radiative impact, and related cloud adjustments, can be wrongly estimated in both sign and magnitude when not considering this dust component separation.