# Aerosol Forcing from Ground-Based Synergies over a Decade in Barcelona, Spain

Daniel Camilo Fortunato dos Santos Oliveira, Michaël Sicard, Alejandro Rodríguez-Gómez, Adolfo Comerón, Constantino Muñoz-Porcar, Cristina Gil-Díaz, Oleg Dubovik, Yevgeny Derimian, Masahiro Momoi, Anton Lopatin

► To cite this version:

Daniel Camilo Fortunato dos Santos Oliveira, Michaël Sicard, Alejandro Rodríguez-Gómez, Adolfo Comerón, Constantino Muñoz-Porcar, et al.. Aerosol Forcing from Ground-Based Synergies over a Decade in Barcelona, Spain. Remote Sensing, 2025, 17, 10.3390/rs17081439. hal-05041767

HAL Id: hal-05041767

https://univ-reunion.hal.science/hal-05041767v1

Submitted on 22 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Aerosol Forcing from Ground-Based Synergies over a Decade in Barcelona, Spain

Daniel Camilo Fortunato dos Santos Oliveira<fnref n="1" />, Michël Sicard <fnref n="1" />,<fnref n="2" />, Alejandro Rodríguez-Gómez <fnref n="1" />, Adolfo Comerón <fnref n="1" />,*, Oleg Dubovik <fnref n="3" />, Yevgeny Derimian <fnref n="3" />, Masahiro Momoi <fnref n="4" /> and Anton Lopatin <fnref n="4" />

**Abstract:** This research aims to estimate long-term aerosol radiative effects by combining radiation and Aerosol Optical Depth (AOD) observations in Barcelona, Spain. Aerosol Radiative Forcing and Aerosol Forcing Efficiency (ARF and AFE) were estimated by combining shortwave radiation measurements from a SolRad-Net CM-21 pyranometer (level 1.5) and AERONET AOD (level 2), using the direct method. The shortwave AFE was derived from the slope between net solar radiation and AOD at 440, 675, 879, and 1020 nm, and the ARF was computed by multiplying the AFE by AOD at six solar zenith angles (20°, 30°, 40°, 50°, 60°, and 70°). Clear-sky conditions were selected from all-skies days by a quadratic fitting. The aerosol was classified to investigate the forcing contributions from each aerosol type. The aerosol classification was based on Pace and Toledano's thresholds from AOD vs. Ångström Exponent (AE). The GRASP inversions were performed by combined AOD, radiation, Degree of Linear Polarization (DoLP) by zenith angles from the polarized sun-sky-lunar photometer and the elastic signal from the UPC-ACTRIS lidar system. The long-term AFE and ARF are both negative, with an increasing tendency (in absolute value) of +24% (AFE) and +40% (ARF) in 14 years. The yearly AFE varied from -331 to -10 Wm⁻<fnref n="2" />τ⁻<fnref n="1" />, and the ARF varied from -64 to -2 Wm⁻<fnref n="2" />, associated with an AOD (440 nm) from 0.016 to 0.690. The three types of aerosols on clear-sky days are mixed aerosols (61%), desert dust (10%), and urban/industrial-biomass burning aerosols (29%). Combined with Gobbi's method, this classification clustered the aerosols into four groups by AE analysis (two coarse- and two fine-mode aerosols). Then, the contribution of the aerosol types to the ARF showed that the desert dust forcing had the largest cooling effect in Barcelona (-61.5 to -37.4 Wm⁻<fnref n="2" />), followed by urban/industrial-biomass burning aerosols (-40.4 to -20.4 Wm⁻<fnref n="2" />) and mixed aerosols (-31.8 and -24.0 Wm⁻<fnref n="2" />). Regarding the comparison among Generalized Retrieval of Atmosphere and Surface Properties (GRASP) inversions, AERONET inversions, and direct method estimations, the AFE and ARF had some differences owing to their definitions in the algorithms. The DoLP, used as GRASP input, decreased the ARF overestimation for high AOD.

**Keywords:** aerosol forcing; AOD; flux; shortwave; CM-21; CNR4; GRASP

<figure><img src="image_2.png" /></figure>

check for
updates

Academic Editor: Matthew McGill

Received: 19 March 2025

Revised: 10 April 2025

Accepted: 14 April 2025

Published: 17 April 2025

**Citation:** Oliveira, D.C.F.d.S.; Sicard, M.; Rodríguez-Gómez, A.; Comerón, A.; Muñoz-Porcar, C.; Gil-Díaz, C.; Dubovik, O.; Derimian, Y.; Momoi, M.; Lopatin, A. Aerosol Forcing from Ground-Based Synergies over a Decade in Barcelona, Spain. *Remote Sens.* **2025**, *17*, 1439. https://doi.org/10.3390/rs17081439

**Copyright:** © 2025 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license
(https://creativecommons.org/
licenses/by/4.0/).

<footnote n="1">¹ CommSensLab, Department of Signal Theory and Communications, Universitat Politècnica de Catalunya, 08034 Barcelona, Spain; daniel.camilo.fortunato@upc.edu (D.C.F.d.S.O.); michael.sicard@univ-reunion.fr (M.S.); alejandro.rodriguez.gomez@upc.edu (A.R.-G.); constantino.munoz@upc.edu (C.M.-P.); cristina.gil.diaz@upc.edu (C.G.-D.)</footnote>

<footnote n="2">² Laboratoire de l'Atmosphère et des Cyclones, Université de la Réunion, 97744 Saint-Denis, France</footnote>

<footnote n="3">³ CNRS, UMR 8518—LOA—Laboratoire d'Optique Atmosphérique, Université de Lille, 59000 Lille, France; oleg.dubovik@univ-lille.fr (O.D.); yevgeny.derimian@univ-lille.fr (Y.D.)</footnote>

<footnote n="4">⁴ Generalized Retrieval of Atmosphere and Surface Properties, GRASP-SAS, 59260 Lezennes, France; masahiro.momoi@grasp-earth.com (M.M.); anton.lopatin@grasp-sas.com (A.L.)</footnote>

<footnote>* Correspondence: comeron@tsc.upc.edu</footnote>
# 1. Introduction

The interaction between solar radiation and aerosols affects the Earth's radiation budget by direct effect, attenuating the solar radiation through the physical processes of absorbing and scattering [1]. The aerosol attenuation is spectrally dependent, i.e., the aerosol absorption decreases with wavelength, whereas the aerosol scattering increases with wavelength for mineral dust and decreases for urban pollution [2]. This interaction also affects the Earth's radiation budget by the indirect effect, modifying the cloud microphysics and its lifetime, according to [3,4], among other authors. Any alteration of the overall radiation budget, known as radiative forcing, can cause climate change [5]. This forcing, which depends on the aerosol load and type, can cool or heat the atmosphere and, consequently, the Earth's surface. For instance, a recognized cooling effect on climate is due to the shortwave-radiation scattering from tropospheric aerosols [6]. Even with advances in aerosol monitoring, the aerosol impact on the Earth's radiation has large sources of uncertainties [7] due to the variability of aerosol optical properties and the spatial-temporal distribution [8]. In the last decades, the synergies between ground-based remote sensing instruments have been widely used to understand aerosols' effects on solar radiation at the Earth's surface. Ref. [9] is an example of pioneering research on utilizing pyranometers and radiometers to estimate aerosol forcing.

Two well-known parameters to quantify the contribution of aerosol radiative effects on the climate are Aerosol Radiative Forcing (ARF) and Aerosol Forcing Efficiency (AFE). The first is the difference between the shortwave net irradiances with and without aerosols, i.e., the change in the radiation flux produced by the particles [10]. The second parameter is the change in the net radiation per unit of Aerosol Optical Depth (AOD) [11,12] at a certain wavelength. ARF and the AFE reinforce the AOD's role in the Earth's radiation budget. AOD is an aerosol columnar property representing the aerosol load in the atmospheric column and, according to [13], it enhances or reduces solar radiation at the surface (brightening and dimming effects, respectively). For example, some of the first studies using AOD from radiometers to measure aerosol forcing were developed by ref. [1,14].

Regarding aerosol forcing estimations in Spain, many studies have been conducted in the south, especially in Granada and the northeastern region. For instance, ref. [15] found integrated aerosol forcing in the whole solar spectrum in Almería; ref. [16] evaluated the surface efficiency at 670 nm and forcing in the visible spectral range (0.4–0.7 µm at 20–80°) during the 2003 heat wave with desert dust events in Granada; ref. [17] obtained diurnal efficiency values at 380 and 440 nm at Armilla and Sabinas stations, in Granada; ref. [18] analyzed aerosol forcing during the strongest desert dust intrusion mixed with smoke at El Arenosillo station and Palencia for 440 nm at 53–75°; ref. [19] calculated efficiencies in the UV erythemal range at 20–55° in Granada; at the same place, ref. [20] also calculated aerosol forcing and efficiencies for 380–870 nm at Solar Zenith Angle (SZA) < 65°; ref. [21] obtained aerosol forcing on the surface (0.31–2.80 µm at 20–80°) during African desert dust events in Granada as well; ref. [5] computed aerosol efficiency for 0.305–2.800 µm at 15° and aerosol forcing at 675 nm in Granada; ref. [22] quantified the direct aerosol radiative effect for two case studies within a dust transport episode in Granada, using the Global Atmospheric ModEl (GAME) code [23,24] with three different input data; and [25] analyzed the efficiency for 0.28–2.80 µm and 0.4–0.7 µm at 15–75° and forcing at 500 nm in Granada.

In northeastern Spain, a few studies were located in Barcelona, such as [26], who obtained instantaneous radiative forcing by the GAME code in the longwave spectral range (4–50 µm) for 11 dust outbreaks in Barcelona, and [27], who calculated mean annual efficiencies and maximum forcing in summer for 0.3–2.8 µm at 50–80° in Palma (Balearic Islands). Additionally, the aerosol radiative effects by long-term period, aerosol type, and wavelengths have not yet been studied in Barcelona. Therefore, this research is driven by
the need to investigate the contribution of aerosols to the radiation budget in Barcelona, in
northeastern Spain.

This study aims to estimate long-term aerosol radiative effects by combining radiation
measurements from pyranometers and AOD observations from photometers (ground-
based measurements). This research focuses on the global solar radiation at the Bottom of
the Atmosphere (BOA), i.e., the downward shortwave irradiance on the Earth's surface. It is
also essential to mention that this research focuses on employment of the direct method [28]
through 14 years of AOD and net flux data; thus, a deep analysis involving satellites and
other aerosol parameters is not conducted in this article.

The manuscript is organized as follows: Section 2 presents an overview of the experi-
mental site and the instrumentation used in this paper. Section 3 presents the methodology
applied to clear-sky identification, forcing estimation and GRASP retrievals. Section 4
presents the main results and discussions of the AOD and flux trends, long-term AFE and
ARF per aerosol types, and the comparison between radiative properties from AERONET
inversions, the direct method, and GRASP inversions. The conclusions are in Section 5 and
complementary information is provided in Appendices A and B.

## 2. Materials

Barcelona is located in northeast Spain in the Mediterranean region. It is a coastal
area with almost 2 million inhabitants and its atmosphere is influenced by marine aerosol,
pollen, urban and biomass-burning aerosols from local and external sources, and desert
dust [27,29–32]. The Barcelona ground-based measurement site, at the Universitat
Politècnica de Catalunya (UPC, Barcelona, Spain), is a member of the AERosol RObotic
NETwork (AERONET) [33], Solar Radiation Network (SolRad-Net), the National Aeronau-
tics and Space Administration Micro-Pulse Lidar Network (NASA MPLNET) [34], and the
Aerosols, Clouds, and Trace Gases Research Infrastructure Network/European Aerosol
Research Lidar Network (ACTRIS/EARLINET) [35].

The CommSensLab Optical Remote Sensing group at the UPC currently operates three radiometers (a Skye Instruments SKE-510 PAR-Photosynthetically Active Radiation, Llandrindod Wells, United Kingdom, a Kipp and Zonen CM-21 pyranometer, and a Kipp and Zonen CNR4 Net Radiometer, Delft, the Netherlands); a Polarized sun-sky-lunar photometer (CE318-TP9) manufactured by Cimel Electronique, Paris, France; a multi-wavelength high-power lidar system, and a NASA Polarized Micro-Pulse Lidar (P-MPL) system, Droplet Measurement Technologies and Tesscorn AeroFluid, Bangalore, India.

2.1. UPC Photometer

Over three decades, the AERONET has provided measurements of columnar AOD and
aerosol characteristics at high temporal resolution [36], employing worldwide standard-
ization of the automatic Sun photometers [33]. Additionally, the UPC CommSensLab has
maintained a Sun photometer on the rooftop of the B3 building (41.38925° N, 2.11206° E) for
20 years (2004–2023), in which the data are available at Level 2.0 (L2) and Version 3 (V3) [36]
where a pre- and post-field calibration, an automatic cloud-screening algorithm, and au-
tomatic instrument anomaly quality controls assure the high quality of the observations
and inversions. The new UPC polarized sun-sky-lunar photometer provides polarized
radiances and lunar direct measurements [32] at Level 1.5 (L1.5).

The V3 has a larger optical air mass range than the previous version. The maximum
air mass of 7.0 increases the number of solar measurements occurring in the early morning
and evening, contributing to additional AOD measurements [36] at larger SZAs. The new
quality-control changes of the V3 play a role in this research analysis due to performing
more observations in the early mornings and evenings that match the radiative fluxes.
## 2.2. UPC Pyranometers

The CNR4 Net Radiometer was installed in February–March 2021 at the UPC. It is an instrument that includes four sensors measuring the energy balance between incoming short-wave and long-wave far infrared radiation versus surface-reflected short-wave and emitted long-wave radiation [37]. The four radiometers consist of a pyranometer pair (one facing upward and the other facing downward) that measures short-wave radiation and a pyrgeometer pair (facing upward and downward) that measures long-wave radiation. Since its installation, the UPC CNR4 has been calibrated three times, as displayed in Table 1. From the first calibration in 2016 to the last one in 2024, the calibration coefficients have remained stable over time, suggesting a minimal instrumental drift and a reliable measurement performance. The third calibration (December 2024) is the current one.

<table><caption>**Table 1.** The UPC CNR4 calibrations and its calibration coefficients. Pl stands for polynomials, T for the temperature, and S for the sensors: $T = C3 \times 3 + C2 \times 2 + C1 \times + C0$ and $S1-4 = (1000/sensor \ sensitivity)x$. The sensor sensitivity (µV/Wm⁻²) is provided after the sensor calibrations. S1 and S2 are the facing-upward and facing-downward pyranometers, respectively, S3 and S4 are the facing-upward and facing-downward pyrgeometers, and x is the experimental voltage in mV.</caption><thead><tr><th>Pl</th><th>1st Calibration (November 2016)</th><th>2nd Calibration (February 2020)</th><th>3rd Calibration</th></tr></thead><tbody><tr><td>T</td><td>1.665e-8x³ - 5.576e-5x² + 0.0991x - 44.42</td><td>1.529e-8x³ - 5.168e-5x² + 0.0955x - 43.52</td><td>1</td></tr><tr><td>S1</td><td>66.1813x</td><td>66.401x</td><td>65.7462x</td></tr><tr><td>S2</td><td>77.1605x</td><td>77.639x</td><td>77.101x</td></tr><tr><td>S3</td><td>78.7402x</td><td>81.234x</td><td>80x</td></tr><tr><td>S4</td><td>78.1861x</td><td>80.58x</td><td>78.064x</td></tr></tbody></table>

The UPC also maintains a SolRad-Net Pyranometer (CM-21 model) that has provided irradiance (Wm⁻²) on a plane surface since 2009 at L1.5, which is free of any operational problems (amplifier malfunctions, a contaminated or dirty dome, and a non-level sensor as examples) and it has shown minimal sensor drift of less than 1% according to Data Quality Assessment and Data Level Designations from NASA (https://solrad-net.gsfc.nasa.gov/system_info_additional.html, accessed on 12 June 2024). In contrast to the UPC CNR4, the CM-21 has only one sensor facing upward to measure the incoming short-wave radiation. Table 2 compares the technical specifications of both radiometer models.

<table><caption>**Table 2.** Technical specifications of the UPC CNR4 and the CM-21 model.</caption><thead><tr><th>Characteristics</th><th>UPC CNR4</th><th>CM-21</th></tr></thead><tbody><tr><td>Spectral range (nm)</td><td>300–2800 (50% points)<br>180° (Upper detectors)</td><td>305–2800 (50% points <sup>1</sup>)<br>355–2200 (95% points)</td></tr><tr><td>Field of view</td><td>150° (Lower detectors) due to a lower sun<br>shield to prevent illumination at low<br>zenith angles</td><td>180°</td></tr><tr><td>Directional error (Wm<sup>-2</sup>)</td><td>&lt;20 (angles up to 80° with 1000 Wm<sup>-2</sup> beam<br>radiation). Combined zenith and azimuth error<br>from 0° to 80° with 1000 Wm<sup>-2</sup> beam radiation</td><td>&lt;± 10 (1000 Wm<sup>-2</sup> beam radiation)</td></tr></tbody></table>

## 2.3. CNR4 Validation

The instantaneous radiative fluxes from the UPC CNR4 were compared to those from the CM-21 model, as shown in Figure 1. An excellent correlation between the CM-21 and

<footnote n="1">¹ The system to measure temperature did not change from the 2nd to the 3rd calibration, so T coefficients remain the same.</footnote>

<footnote n="1">¹ “Points” are the wavelengths where the output of the instrument is 50% or 95% reduced with 100% input (https://www.kippzonen.com/FaqCategory/3/Pyranometers, accessed on 3 April 2025).</footnote>
the UPC CNR4 pyranometers facing upward is observed in the same figure, in which the linear regression (red line) has a crescent slope (0.96). Despite the Root Mean Square Error (RMSE) having a high value (69.18 Wm⁻²), the Normalized Mean Bias (NMB) shows a CNR4 non-significant overestimation (slight positive bias of 0.5%) and a coefficient of determination (r²) of 0.94, indicating that 94% of the variance of the UPC CNR4 data can be explained by the variance of the CM-21 data, i.e., the UPC CNR4 pyranometer is properly measuring the radiative fluxes. The high value of RMSE is expected due to the high dispersion of the points (dark blue dots) in the range of 600–1200 Wm⁻².

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison between the CM-21 model from the SolRad-Net and the UPC CNR4 observations from February 2021 to August 2023. The red line is the linear fitting between data, and the dashed black line is the 1:1 reference. The color bar shows the density of the data with colors on a logarithmic scale varying from 0 to 1 (100%).</figcaption></figure>

## 2.4. UPC-ACTRIS Lidar System

The UPC-ACTRIS lidar system emits elastic signals at 355, 532, and 1064 nm by an Innolas Spitlight 400 laser with a spatial resolution of 3.75 m, a divergence of 28 mrad, a pulse duration of 3.6 ns, and a repetition rate of 20 Hz. The system setup is monostatic, vertical, and biaxial.

The receiver is composed of a Celestron CGE 1400 telescope, a focal length of 3.91 m, a telescope aperture of 0.35 m (diameter), and photomultiplier tubes and avalanche photodiodes detectors. The reception channels are the same for the elastic wavelengths (355, 532, and 1064 nm), two pure-rotational Raman channels at 354 and 530 nm, two depolarization channels at 355 and 532 nm, and one vibro-rotational channel at 407 nm for the water-vapor mixing ratio (for more details, see [32]).

# 3. Methods

## 3.1. All-Skies Time Series

This research required a differentiation between all-skies and clear-sky days. The all-skies radiative fluxes were provided by the SolRad-Net at L1.5 from 2009 to 2023. The observations were selected by (i) diurnal cycles with 80% of measurements during hours of sunshine, ensuring more complete diurnal cycles, and (ii) radiative flux values at a minimum SZA ± 1° (nearest to the zenith). The data from the pandemic year (2020) are not included in this work because there are not enough measurements from the instruments in Barcelona.
## 3.2. Clear-Sky Selection

Before calculating the aerosol forcing, a selection of the cloudless sky conditions was performed by a quadratic fitting [38–41]. This simple approach identifies only clear-sky days during a whole day by a regression equation, $ax^2 + bx + c$, where a, b, and c are its coefficients. The quadratic fitting was applied to the radiative fluxes on the surface and their respective SZA cosine (x-axis), which greatly influences the radiative fluxes [42].

To execute the fitting and to guarantee a homogenous atmosphere, i.e., the symmetry between morning and afternoon fluxes, the following criteria were applied to the diurnal cycles: (i) selection of diurnal cycles with 80% of measurements during hours of sunshine, ensuring more complete diurnal cycles; (ii) elimination of SZAs greater than 75° to avoid interferences from surrounding reliefs and to avoid errors from the cosine response on the measurements [25]; (iii) application of $r^2$ greater than 0.98 for the best fit between cosines of the SZA and radiative fluxes; and (iv) application of a seasonal adaptive RMSE (i.e., winter factor of 1, spring factor of 1.18, summer factor of 2.1, and autumn factor of 1.1) by obtaining the rates between the period average of all winter solstices (maximum flux values around noon) and the period averages of all summer solstices and autumn-spring equinoxes (maximum flux values around noon). Those criteria ensure less variability in the fluxes caused by possible clouds in the morning or afternoon period. Figure 2 is an example of a day classified as a clear-sky day by the quadratic fitting, where the $r^2$ is 0.999 and the RMSE is 7.99 Wm$^{-2}$.

<figure><img src="image_4.png" /><figcaption>Figure 2. Example of a selected clear-sky day (16 April 2022) employing the quadratic fitting (QF) in red and the SolRad-Net measurements in blue as a function of (left) cos(SZA) and (right) UTC.</figcaption></figure>

## 3.3. AFE-ARF Estimation

The ARF is the difference between the net irradiance on the surface ($I_{net}$) and the net irradiance without aerosol contribution ($I_{net,0}$), i.e., the ARF is the change in the net flux due to variations in the aerosol properties regarding an aerosol-free atmosphere. The ARF can be calculated by the equation below:

$$ARF = I_{net} - I_{net,0} \qquad (1)$$

The $I_{net}$ component is computed as follows:

$$I_{net} = (1 - A)I \qquad (2)$$
where A is the surface albedo calculated as the period mean of the ratio between outcoming and incoming radiations from the UPC CNR4 by SZA and I is the measured downward irradiance on the surface. To cover the greatest amount of SZA, six angle groups were fixed (20°, 30°, 40°, 50°, 60°, and 70°) with a range of ±1° because larger angle ranges can have higher variation in irradiation. For Barcelona, A is the mean albedo over the 3-year CNR4 observations at SZA ranges (20 ± 1°, 30 ± 1°, 40 ± 1°, 50 ± 1°, 60 ± 1°, and 70 ± 1°) with 0.104, 0.104, 0.106, 0.108, 0.115, and 0.124, respectively.

As aerosol forcing is strongly affected by the aerosol load [43], it can be assessed by Aerosol Forcing Efficiency (AFE), which is the change in the ARF per unit in AOD for a certain wavelength (λ) [12]:

$$AFE(SZA, \lambda) = \frac{dARF}{dAOD} \qquad (3)$$

Equation (3) can be rewritten based on Equation (1):

$$AFE(SZA, \lambda) = \frac{d(I_{net} - I_{net,0})}{dAOD} \qquad (4)$$

This estimation avoids some errors due to the instrumental or model offsets, more common in other methods to compute the AFE; thus, the present study used the direct method determined by [28] to estimate it as the slope derived from the linear regression between the I<sub>net</sub> and the coincident AOD at an established SZA (Equation (4)) from the clear-sky condition (Section 3.2). The direct method was also well applied by [5,25,44–46], among others. Furthermore, the efficiency was calculated for the standard wavelengths of the AERONET inversions (440 nm, 675 nm, 870 nm, and 1020 nm). It is important to mention that the outliers greater than three times the median were deleted from the database to remove possible residual contamination by clouds, and then the fitting was calculated.

The advantage of this method lies in direct estimates of the AFE from observations, without further assumptions about the radiative fluxes under aerosol-free conditions [5], i.e., there is no dependency on the aerosol load under aerosol-free conditions. Hence, the direct method is derived from Equation (4) without the I<sub>net,0</sub>:

$$AFE(SZA, \lambda) = \frac{dI_{net}}{dAOD} \qquad (5)$$

Once the AFE was estimated, the ARF was derived by multiplying the AFE by each AOD [5,25,45,46] at the corresponding SZA for each wavelength (440 nm, 675 nm, 870 nm, and 1020 nm). Then, the final AFE is averaged by SZA and by year. Furthermore, the ARF was also estimated for three groups of aerosols at 440 nm: Mixed Aerosols (MAs), Desert Dust (DD), and Urban/Industrial-Biomass Burning aerosols (UI-BB). To better determine the thresholds of the 'AOD vs. Ångström exponent (AE)' classifications, two aerosol types [47,48] were adopted, which fitted what was proposed in this study, even if there is no support of chemical speciation or other aerosol properties for ensuring the aerosol types in a mixed-source urban environment. Those classifications were chosen because (i) both sites (Lampedusa Island in Italy [47] and Huelva City in Spain [48]) are strongly influenced by marine aerosols, desert dust, and local aerosols, as in Barcelona, allowing mixtures between them; (ii) both locations are in the Mediterranean region; and (iii) there is not such a complete kind of aerosol classification for Barcelona.

### 3.4. GRASP/GARRLiC Algorithm and BRDF-BPDF Calculations for Barcelona

The Generalized Retrieval of Atmosphere and Surface Properties/Generalized Aerosol Retrieval from Radiometer and the Lidar Combined data (GRASP/GARRLiC) code is a powerful unified algorithm for characterizing atmospheric and surface properties by com-
bining a variety of remote sensing observations [49]: satellites, lidar systems, photometers, ceilometers, and others. GRASP code is based on a PARASOL (Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar) retrieval algorithm developed by [50] and also relies on the heritage of AERONET retrieval advances [51–54]. The GARRLiC code [55], part of the GRASP code, inverts coincident lidar and radiometer observations and derives a united set of aerosol parameters. It assumes a bi-component aerosol model (different complex refractive indexes) as a mixture of polydisperse spheres and randomly orientated spheroids, and the microphysical properties are height-independent.

To retrieve aerosol and surface properties, the GRASP code needs radiance and polarization measurements and requires accurate models for the Bidirectional Reflectance Distribution Function (BRDF) and the Bidirectional Polarization Distribution Function (BPDF), for example, the Ross–Li model [56–58] and Maignan–Breon model [59,60], respectively. Both reflectance parameters were validated and extensively compared with the Polarization and Directionality of Earth’s Reflectances (POLDER) data by [60,61]. The Ross–Li model is a linear combination of three-surface scatterings (isotropic, volumetric, and geometric optical), while the Maignan–Breon model is a linear one-parameter BPDF.

This study obtained the BRDF-BPDF surface values from POLDER satellite climatology [62] over Barcelona. The BRDF isotropic values were interpolated for the photometer channels. In contrast, the volumetric, geometric, and BPDF surface values are considered spectrally independent in the visible and infrared regions [63]. After the parametrization, the broadband solar flux and aerosol radiative effect calculations at the BOA were conducted using the code and principles described in [64] and their implementation in the GRASP algorithm [49,50,65]. The radiative code [64] also assumed nonspherical-spherical particles. Three case studies were used to compare the outputs from the AERONET inversions, direct method (synergy between pyranometer and photometer data), and GRASP inversions. The fluxes from the pyranometer and the photometer over 14 years were also compared for 958 inversions under the clear-sky condition.

The GRASP inversions were also performed for clear-sky days with AOD, normalized radiances, Degree of Linear Polarization (DoLP), and Range-Corrected Signals (RCS) as inputs for the same combinations proposed by [32]. The combinations are shaped by data from the UPC-ACTRIS lidar system and the sun-sky-lunar photometer for five situations, D0, D0+L, D0P+L, D1+L, and D1P+L, where D is the daytime measurement, 0 refers to standard photometer channels, L refers to three RCS from the UPC-ACTRIS lidar system (355, 532, and 1064 nm), P refers to DoLP (more information in [32]), and 1 is the standard photometer plus new channels (380, 500, and 1640 nm). The AERONET inversions had less than 3% sky residuals and AOD (440 nm) ≥ 0.2 to ensure good retrievals. The lidar signals were processed by the Single Calculus Chain [66–68]. More details are found in [32]. With this synergy, more complete information about the radiation attenuation (vertically distributed and column-integrated properties) by aerosols is provided for the forcing calculations. Additionally, improvement in the aerosol property retrievals (e.g., AOD and single scattering albedo) relevant to the forcing are expected by adding polarization data [54,69–71] to the synergy [32].

To simplify the comparison among the broadband fluxes from photometer inversions and pyranometer observations, the parameters were plotted as a function of the SZA in which only three angle groups were fixed with a range of ±1° (50°, 60°, and 70°) due to the AERONET inversions that do not perform measurements at small SZAs [33,51,72,73]. The Percentage Error (PE) was calculated to indicate the magnitude of the deviation between the inversions and the observations, as in Equation (6):
$$PE = \left( \frac{\text{invertedvalue} - \text{observedvalue}}{\text{observedvalue}} \right) 100\% \qquad (6)$$

Therefore, with the methodology proposed in this section, the present research is a pioneer for calculating aerosol radiative effects integrated from 305 to 2800 nm at 20° to 70° during 14 years in Barcelona, Spain.

# 4. Results

## 4.1. Trends

The study period from 2009 to 2023 (availability of coincident AERONET and pyranometer observations) was analyzed by the Mann-Kendall (MK) test and Sen's slope. As indicated by [74], the MK test is used to identify increasing or decreasing long-term monotonic trends where the data are not serially correlated, and the Sen's slope estimator, which is less sensitive to outliers and extreme values, is based on the median of slopes calculated from all possible data pairs. The MK test considers seasonality and reduces the chances of falsely detecting a trend where one does not exist [74]. The MK test was applied with a 95% confidence interval (CI95) for the AOD and the all-skies radiative fluxes, and it was performed by [75,76] functions. Linear regressions fit the clear-sky fluxes by season with a CI95, and the slope uncertainties were also estimated with a CI95.

### 4.1.1. AOD Trends

Figure 3 shows the AOD trends estimated by the MK test and the Sen's slope, for which the monthly means were calculated by the L2-V3 AERONET database (https://aeronet.gsfc.nasa.gov/, accessed on 6 June 2024). The former period (black trend line) represents the whole L2-V3 AERONET period and the latter (red trend line) refers to the study period of this research.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time series of AOD monthly means from the AERONET photometer at L2-V3. The blue dots are the monthly means, the black line is the 2004–2023 trend, and the red line (2009–2023) is the trend of this research with a CI95. N is the sample size, yr is an abbreviation for year, and Δτ is the change in the AOD over the period considered.</figcaption></figure>
The Sen's slopes for both AOD trends (2004–2023 and 2009–2023) computed with a CI95 from the MK test are shown in Figure 3. The 2004–2023 trend (black line) has a Sen's slope of −0.0038 yr⁻¹ with a change of −0.07 (31.8%) over the whole period (19 years). However, the 2009–2023 trend (red line) has a Sen's slope of −0.0015 yr⁻¹ with a change of −0.02 (10.5%) over 15 years.

4.1.2. Trends, Seasonal Distributions, and Variations of the Radiative Fluxes

The trends in Figures 4 and 5 were calculated from monthly means of the daily solar noon (smallest SZA ±1°) for all-skies days from 2019 to 2023 and clear-sky days by seasons and the same period, respectively. The seasonal trends were calculated to better represent the temporal evolution of the radiative fluxes, avoiding a biased tendency due to years with few clear-sky days. As depicted in Figure 4, the radiative–flux uptrend (red line) for the all-skies condition has an increase of +1.22 Wm⁻² (+0.2%) and a crescent Sen’s slope of +0.085 Wm⁻²yr⁻¹.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** All-skies trend (red line) calculated by MK and Sen's slope for the monthly means of the SolRad-Net pyranometer at L15 (CM-21 model). The trend has a CI95. N is the sample size, yr is the abbreviation for year, and ΔF is the change in the radiative fluxes from 2009 to 2023.</figcaption></figure>

The quadratic fitting as a clear-sky method could identify 614 days (17.4%) from the all-skies database (Section 3.1), which made it possible to verify the radiative fluxes and AOD trends for clear-sky days by season, as shown in Figure 5. The changes in winter trends have values of +11.85 Wm⁻² (+2.4%) and -0.006 (-7.2%) with slopes of +1.18 Wm⁻²yr⁻¹ and -0.0005 yr⁻¹ for radiative fluxes and AOD, respectively. The spring trends have a flux change of -77.97 Wm⁻² (-8.6%) with a slope of -6 Wm⁻²yr⁻¹ and an AOD change of +0.02 (+14.3%) with a slope of +0.002 yr⁻¹. The summer trends present a flux change of -3.63 Wm⁻² (-0.4%) with a slope of -0.3 Wm⁻²yr⁻¹ and an AOD change of +0.03 (+15.8%) with a slope of +0.0028 yr⁻¹. The last seasonal trends present a flux increase of +82.69 Wm⁻² (+12.8%) with a slope of +6.89 Wm⁻²yr⁻¹ and the largest AOD drop of -0.04 (-23.5%) with a slope of -0.0037 yr⁻¹.

The seasonal diurnal cycles and daily variation of the radiative fluxes are presented in Figure 6. The continuous lines are the radiative–flux maximums, and the dashed lines are the radiative–flux means followed by their standard deviations (shaded areas). The blue color refers to all-skies days, and the green color refers to clear-sky days. In all seasons,
the peak flux is around noon for both sky conditions. For clear skies, the maximums of sun hour reach 683.62, 1040.05, 1050.77, and 848.52 Wm⁻² in winter, spring, summer, and autumn, respectively.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Seasonal clear-sky trends of radiative fluxes and AOD monthly means from 2009 to 2023 with a CI95. The blue dots are the flux monthly means, the orange dots are the AOD monthly means, the continuous blue and orange lines are their respective trends, N is the sample size, ΔF is the change in the radiative fluxes, Δτ is the change in the AOD, and yr is the abbreviation for year.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Diurnal cycle of radiative fluxes by seasons in Barcelona. The maximum (Max), the mean, and the standard deviation (Std) are shown for all-skies (blue lines) and clear-sky (green lines) conditions.</figcaption></figure>

The clear-sky means are higher than the all-skies means for all seasons, followed by small variations (green shaded areas) due to the non-interference of the clouds, however, the cloud disturbances are visually expressive in the large variation of the all-skies standard deviations. For spring and summer, the solid green lines (the maximums) show some
spikes due to the disturbances produced by fractional cloud cover and by the aerosol loads that are more often in those seasons than in others.

## 4.2. Long-Term of the AFE and the ARF

Figure 7 shows the net radiance fluxes of the clear-sky days from 2009 to 2023 as a function of the AOD at 440 nm and their linear fittings for each SZA group (20 ± 1°, 30 ± 1°, 40 ± 1°, 50 ± 1°, 60 ± 1°, and 70 ± 1°). The AFE-ARF long-terms at 675, 870, and 1020 nm are in Appendix A as additional information. All linear regressions (LRs) have negative slopes varying between −331 and −10 Wm⁻²τ⁻¹, i.e., Barcelona has a dominant cooling effect in the analyzed period. The slopes are the AFE according to the direct method (Section 3.3). The AOD related to the efficiencies varies from 0.016 to 0.690.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The AFE Long-term from 2009 to 2023 for six SZA groups at 440 nm. The linear regression (LR) slopes were calculated with a CI95.</figcaption></figure>

Table 3 summarizes the AFEs and their confidence interval at 95% by fixed SZAs at 440 nm. The AFE tables for 675, 870, and 1020 nm are also in Appendix A. The CI95 values decrease at higher SZAs (thicker atmosphere) due to larger data points than at lower SZAs, agreeing with [45]. For example, in 2009, the CI95 corresponded to 57.9% of the AFE at 20°, decreasing to 18.4% at 70°, and it decayed from 83.6% to 10.9% in 2021 (Table A7). The seasonal variability of net fluxes also contributes to a CI95 of the AFE, as shown in Figure 6, where summer and spring have maximum fluxes for clear-sky days. On the one hand, the number of data points contributes to a large CI95 for years with fewer observations and fewer clear-sky days (Table A8) as occurred in 2015, 2016, and 2017, representing 3.7%, 3.9%, and 2.7% of the instantaneous observations and 3.5%, 3.9%, and 4.2% of the clear-sky days, respectively. On the other hand, 2022 has the most observations (22.7%) with 72 clear-sky
days (13.2%), as shown in Table A8, and has the smallest CI95, varying from 5.5% to 29.8% of the AFE.

<table><caption>**Table 3.** The AFE on the surface is a function of the AOD at 440 nm for six SZA groups. The slope was estimated with a CI95. Y stands for years.</caption><thead><tr><th rowspan="2">Y</th><th colspan="6">AFE, Wm<sup>-2</sup> τ<sup>-1</sup></th></tr><tr><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr></thead><tbody><tr><td>2009</td><td>-158 ± 29</td><td>-175 ± 35</td><td>-195 ± 45</td><td>-153 ± 47</td><td>-159 ± 52</td><td>-145 ± 84</td></tr><tr><td>2010</td><td>-192 ± 27</td><td>-234 ± 39</td><td>-176 ± 64</td><td>-275 ± 57</td><td>-168 ± 49</td><td>-182 ± 30</td></tr><tr><td>2011</td><td>-168 ± 36</td><td>-190 ± 51</td><td>-153 ± 42</td><td>-305 ± 197</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-145 ± 26</td><td>-148 ± 40</td><td>-116 ± 43</td><td>-158 ± 37</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-224 ± 48</td><td>-331 ± 71</td><td>-270 ± 85</td><td>-170 ± 49</td><td>-86 ± 96</td><td>-91 ± 29</td></tr><tr><td>2014</td><td>-233 ± 42</td><td>-227 ± 42</td><td>-98 ± 41</td><td>-85 ± 52</td><td>-75 ± 68</td><td>-104 ± 45</td></tr><tr><td>2015</td><td>-188 ± 33</td><td>-96 ± 49</td><td>-140 ± 30</td><td>-127 ± 33</td><td>-148 ± 36</td><td>-127 ± 30</td></tr><tr><td>2016</td><td>-107 ± 46</td><td>-122 ± 58</td><td>-117 ± 63</td><td>-144 ± 52</td><td>-161 ± 43</td><td>-</td></tr><tr><td>2017</td><td>-76 ± 73</td><td>-80 ± 81</td><td>-102 ± 75</td><td>-235 ± 139</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-110 ± 29</td><td>-152 ± 41</td><td>-196 ± 37</td><td>-108 ± 42</td><td>-10 ± 67</td><td>-138 ± 48</td></tr><tr><td>2019</td><td>-179 ± 28</td><td>-218 ± 39</td><td>-242 ± 50</td><td>-319 ± 42</td><td>-243 ± 57</td><td>-190 ± 34</td></tr><tr><td>2021</td><td>-230 ± 25</td><td>-277 ± 46</td><td>-260 ± 77</td><td>-245 ± 52</td><td>-</td><td>-128 ± 107</td></tr><tr><td>2022</td><td>-202 ± 18</td><td>-260 ± 22</td><td>-208 ± 62</td><td>-210 ± 34</td><td>-174 ± 18</td><td>-163 ± 9</td></tr><tr><td>2023</td><td>-242 ± 32</td><td>-294 ± 53</td><td>-248 ± 112</td><td>-244 ± 81</td><td>-214 ± 62</td><td>-162 ± 43</td></tr></tbody></table>

Table 4 summarizes the ARF means and their standard deviations for each SZA group at 440 nm by year. The ARF tables for 675, 870, and 1020 nm are displayed in Appendix A. The ARF in Barcelona varies from -64 to -2 Wm⁻² in the period. Ref. [27] found an ARF mean value of -26.4 Wm⁻² at 50° < SZA < 60° and values ranging from -64 to -2 Wm⁻² for Mallorca Island, Spain. High variability was also found in the ARF. For instance, the standard deviations rose from 25.9% to 48.5% of the ARF mean at 20° and 70°, respectively, in 2009, and from 26.7% at 20° to 61.1% of the ARF mean at 70° in 2021 (Table A7).

<table><caption>**Table 4.** The ARF means on the surface and their standard deviation at 440 nm for six SZA groups. Y stands for years.</caption><thead><tr><th rowspan="2">Y</th><th colspan="6">ARF, Wm<sup>-2</sup></th></tr><tr><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr></thead><tbody><tr><td>2009</td><td>-33 ± 16</td><td>-41 ± 19</td><td>-41 ± 20</td><td>-36 ± 13</td><td>-36 ± 14</td><td>-27 ± 7</td></tr><tr><td>2010</td><td>-30 ± 19</td><td>-36 ± 22</td><td>-35 ± 16</td><td>-47 ± 27</td><td>-33 ± 13</td><td>-38 ± 15</td></tr><tr><td>2011</td><td>-21 ± 12</td><td>-26 ± 13</td><td>-22 ± 13</td><td>-39 ± 13</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-17 ± 13</td><td>-21 ± 14</td><td>-27 ± 10</td><td>-35 ± 18</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-25 ± 17</td><td>-35 ± 27</td><td>-41 ± 25</td><td>-27 ± 16</td><td>-14 ± 7</td><td>-21 ± 7</td></tr><tr><td>2014</td><td>-28 ± 24</td><td>-28 ± 25</td><td>-18 ± 11</td><td>-14 ± 11</td><td>-13 ± 8</td><td>-15 ± 6</td></tr><tr><td>2015</td><td>-41 ± 20</td><td>-23 ± 8</td><td>-31 ± 12</td><td>-30 ± 12</td><td>-30 ± 14</td><td>-30 ± 11</td></tr><tr><td>2016</td><td>-13 ± 7</td><td>-15 ± 8</td><td>-12 ± 7</td><td>-16 ± 10</td><td>-16 ± 11</td><td>-</td></tr><tr><td>2017</td><td>-12 ± 4</td><td>-12 ± 5</td><td>-17 ± 8</td><td>-34 ± 9</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-23 ± 11</td><td>-30 ± 15</td><td>-40 ± 20</td><td>-28 ± 10</td><td>-2 ± 1</td><td>-29 ± 11</td></tr><tr><td>2019</td><td>-36 ± 17</td><td>-44 ± 20</td><td>-48 ± 21</td><td>-64 ± 27</td><td>-57 ± 25</td><td>-42 ± 17</td></tr><tr><td>2021</td><td>-36 ± 22</td><td>-43 ± 26</td><td>-55 ± 21</td><td>-55 ± 21</td><td>-</td><td>-30 ± 8</td></tr><tr><td>2022</td><td>-28 ± 21</td><td>-34 ± 26</td><td>-38 ± 18</td><td>-42 ± 20</td><td>-35 ± 19</td><td>-35 ± 22</td></tr><tr><td>2023</td><td>-28 ± 21</td><td>-40 ± 28</td><td>-41 ± 21</td><td>-47 ± 31</td><td>-50 ± 28</td><td>-45 ± 22</td></tr></tbody></table>

The AFE has an increase, in absolute values, at intermediate SZA followed by a decrease at higher SZA. This behavior is propagated in the ARF, as displayed in the first plot of Figure 8. The AFE varies between -200 and -175 Wm⁻²τ⁻¹ for SZAs between 40°
and 70°. In the same SZA interval, the ARF decreases (in absolute value) from −37 Wm⁻² at 40° to −27 Wm⁻² at 70°. Similar angular dependence was found by [21,25,45,77,78].

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Angular dependency of the AFE and ARF for six SZA groups at 440 nm and their trends from 2009 to 2023. The linear regression slopes were calculated with a CI95.</figcaption></figure>

Figure 8, in the second plot, also presents the AFE-ARF trends at 20–70° from 2009 to 2023. The AFE has a change ($\Delta$AFE) of –37.39 Wm$^{-2}\tau^{-1}$ (–23.9%) and a slope of –2.67 Wm$^{-2}\tau^{-1}$yr$^{-1}$. The ARF change ($\Delta$ARF) is –10.27 Wm$^{-2}$ (–40.2%) with a slope of –0.73 Wm$^{-2}$yr$^{-1}$. In absolute numbers, both trends increase, indicating that the aerosol forcings throughout the analyzed period in Barcelona are getting stronger with time. Refs. [22,25] also pointed out this effect in Granada, Spain. The single scattering albedo (SSA) of the period ($0.94 \pm 0.03$) shows a considerable absorption aerosol according to Dubovik's climatology [79].

### 4.3. Contribution of the Aerosol Types to the Forcing in Barcelona

This section quantifies the contribution of the aerosol types to the cooling effect in Barcelona. As indicated in Section 3.3, the aerosol in Barcelona was classified according to [47,48] and is displayed in Figure 9. Due to the complexity in determining the aerosol-classification thresholds for AOD-AE classifications, (i) the MA includes pure marine and marine polluted [47], mixed with continental or local-pollution aerosols [48], with AOD > 0 and 1.05 < AE < 1.5, AOD < 0.2 and AE ≤ 1.5, and AOD < 0.1 and AE > 1.5; (ii) the UI-BB includes local pollution and biomass burning [47] in the coast with AOD ≥ 0.1 and AE ≥ 1.5; and (iii) the DD includes dust and some mixture of dust and biomass-burning aerosols [48] with AOD ≥ 0.2 and AE ≤ 1.05. The DD ranges were adapted from [48]. This classification shows a 60.9% predominance of MA, 28.9% of UI-BB, and 10.2% of DD in Barcelona for clear-sky days.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The AOD scatter plot at 440 nm from clear-sky days (2009–2023) according to the aerosol thresholds from [47,48] for Barcelona, Spain. The colors and shapes refer to each aerosol type and the rectangles delimit the classification regions.</figcaption></figure>

The graphical method proposed by [80] was applied to analyze the aerosol fine-mode size, the fractional contribution to the total AOD (η), and the AOD growth due to the increase in coarse particles. Their method depends on the AE difference (δAE = AE(440–675)–AE(675–870)) as a function of the AE calculated between 440 and 870 nm for a bimodal lognormal size distribution with refractive index m = 1.4–0.001i. As stated by [80], the observations of the AOD at 675 nm >0.15 were applied on all plots to avoid AOD error propagation larger than ~30%; hence, the MA aerosol was reduced by 93%, the DD by 8.2%, and the UI-BB by 86.2%. It is worth mentioning that the MA and UI-BB contain marine aerosols (AOD < 0.2) that suffer drastic reductions when applying Gobbi's method. Furthermore, the AOD utilized in this analysis is the one observed at 440 nm. The previous particularity was well performed by [27]. For the present study, the results in the graphical method are presented from the aerosol types by [48] in Figure 9 (AOD vs. AE).

Figure 10 depicts Gobbi's graphical method for the aerosol types by SZA in Barcelona. Most aerosols have positive δAE values (85.4%) with a δAE mean of 0.21 ± 0.11, indicating the presence of different aerosol modes [80,81]. The negative values of δAE have a mean of –0.13 ± 0.09. Regarding the AE, the coarse and fine modes are separated by 1.05 [48], in which the means are 0.62 ± 0.23 and 1.47 ± 0.21, respectively.

As expected, the aerosol amount increases from 17.3% at 20° to 24.3% of the total at 70° for DD, from 12.2% at 20° to 23.3% at 70° for MA, and from 7.8% at 20° to 27.8% at 70° for UI-BB at higher SZAs. The Gobbi's method and AOD-AE classification together characterized some aerosol clusters over Barcelona: (i) The first coarse-mode cluster has an increase in the aerosol size from 0.05 to 0.20 µm, mainly for DD and a few MAs, associated with a small fine fraction (η < 30%), and AE < 0.7 agreeing with [82]; however, a δAE < 0.6 is divergent from theirs. (ii) The second coarse-mode cluster has only DD with radii around 0.10 µm at lower SZAs, exceeding 0.15 µm at higher SZAs, AE between 0.7 and 1.05, δAE < 0.6 at lower SZAs and some δAE < 0 at higher SZAs, and 30% < η < 45% at lower SZAs, exceeding 50% at higher SZAs. (iii) The first fine-mode cluster is a mixture of MAs and UI-BB aerosols with radii centered in 0.10 µm, varying from moderate to high fine
fractions (35% < η < 80%) at lower SZAs, exceeding 0.15 µm, and is associated with a high fine fraction (η > 80%) at higher SZAs, δAE varying from -0.2 to 0.75, and AE > 1.05. Ref. [27] found similar values for fine-mode clusters from observations at insular sites in the western Mediterranean Basin. (iv) A second fine-mode cluster is a small cluster dominated by UI-BB at small SZAs, and it is characterized by AE > 1.05, δAE < -0.2, 75% < η < 90%, and radii around 0.13 µm that agrees with [82] for polluted and continental aerosol in some sites in the Iberian Peninsula and the Western Mediterranean region, including Barcelona. These results reinforce that the angular dependency relies on the aerosol type, as [45] pointed out.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** The AE difference (δAE) as a function of the AE, 440–870 nm, for the classification in Figure 9. The color and dot shapes refer to each aerosol type. The solid black lines refer to the aerosol sizes (radii), the dashed blue lines refer to the fine-mode fractions of the total AOD (η), and the i, ii, iii, and iv red lines are the aerosol clusters.</figcaption></figure>

Finally, the forcing contribution from each aerosol type is shown in Figure 11 and summarized in Table 5. The DD has a cooling effect with negative AFE varying from -207 to -123 Wm⁻²τ⁻¹, increasing its absolute values at 40° and decreasing it at higher SZAs.

<table><caption>**Table 5.** The AFE and ARF mean on the surface by aerosol types at 440 nm and their standard deviation for six SZA groups. AT stands for Aerosol Types, and All stands for all aerosol (first plot in Figure 8).</caption><thead><tr><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr></thead><tbody><tr><th>AT</th><th colspan="6">AFE, Wm<sup>-2</sup> τ<sup>-1</sup></th></tr><tr><th>DD</th><td>-123 ± 39</td><td>-165 ± 39</td><td>-173 ± 67</td><td>-207 ± 52</td><td>-157 ± 29</td><td>-144 ± 12</td></tr><tr><th>MA</th><td>-224 ± 13</td><td>-270 ± 20</td><td>-173 ± 27</td><td>-196 ± 28</td><td>-149 ± 23</td><td>-181 ± 19</td></tr><tr><th>UI-BB</th><td>-138 ± 21</td><td>-190 ± 32</td><td>-170 ± 38</td><td>-157 ± 29</td><td>-94 ± 23</td><td>-142 ± 24</td></tr><tr><th>All</th><td>-175 ± 52</td><td>-200 ± 76</td><td>-180 ± 66</td><td>-198 ± 74</td><td>-144 ± 69</td><td>-143 ± 32</td></tr><tr><th></th><th colspan="6">ARF, Wm<sup>-2</sup></th></tr><tr><th>DD</th><td>-37 ± 9</td><td>-50 ± 14</td><td>-50 ± 12</td><td>-62 ± 16</td><td>-49 ± 14</td><td>-50 ± 19</td></tr><tr><th>MA</th><td>-24 ± 17</td><td>-30 ± 20</td><td>-26 ± 14</td><td>-30 ± 16</td><td>-25 ± 15</td><td>-32 ± 14</td></tr><tr><th>UI-BB</th><td>-29 ± 11</td><td>-40 ± 15</td><td>-38 ± 14</td><td>-35 ± 14</td><td>-20 ± 7</td><td>-29 ± 9</td></tr><tr><th>All</th><td>-27 ± 9</td><td>-31 ± 10</td><td>-33 ± 13</td><td>-37 ± 14</td><td>-29 ± 17</td><td>-31 ± 9</td></tr></tbody></table>
<figure><img src="image_13.png" /></figure>

**Figure 11.** The AFE-ARF means by aerosol type (classification from Figure 9) for six SZA groups at 440 nm. The linear regressions (LR) slopes were calculated with a CI95. The vertical bars are the standard deviations.

The MA varies from −270 to −149 Wm⁻²τ⁻¹ and seems not to have a clear angular dependency because of the alternating drops and rises. The DD has the lowest CI95s, varying from 24% to 38% of the mean, while MA has the highest, ranging from 43.8% to 70.8% of the mean. The AFE from UI-BB ranges from −190 to −94 Wm⁻²τ⁻¹ with an inflection point at 30° also present in the efficiencies calculated by [45]. Its AOD is also quite variable, averaging 0.22 ± 0.08, leading to high CI95s (31–40% of the mean).

As depicted in Figure 11 (the last plot in the second row), the AFE angular dependencies are propagated to the ARF of the UI-BB and DD. The DD has the strongest cooling effect, and its forcing varies from −62 to −37 Wm⁻², increasing in absolute value at 40° and then decreasing at higher SZAs. The ARF of UI-BB varies from −40 to −20 Wm⁻² with an inflection point at 30° and increases (in absolute values) until 50° followed by a decrease at 70°. However, the MA has the smallest ARF at 40–70°, varying between −32 and −24 Wm⁻².

### 4.4. Radiative Fluxes and Aerosol Forcing Comparisons

The comparison between the radiative fluxes from AERONET and CM−21 pyranometer is plotted in Figure 12. The AERONET fluxes slightly overestimate those from CM−21 (on average, PE = 4 ± 1%), showing that both instruments measure comparable fluxes.

The comparisons among the radiative properties and fluxes from AERONET inversions, the direct method, and GRASP inversions are in Table 6. The AERONET fluxes slightly overestimate those from CM-21 in 1.1% and 2.7% (Cases I and II, respectively), which was expected according to [73], but it underestimates the CM-21 flux in 1.6% (Case III). The GRASP combinations (D0, D0+L, D0P+L, D1+L, and D1P+L) overestimate the CM-21 fluxes by ~7% for Cases I and II and ~1% for Case III.

The AFE-ARF comparison between GRASP combinations and the direct method shows (i) high underestimations of the AFE from the direct method in the three case studies (PE ≈ 32%, 35%, and 14% for Case I, II, and III, respectively); (ii) overestimations of the ARF for high AOD (Case I and II) and ARF underestimations by GRASP combinations for low AOD (Case III); and (iii) the GRASP inversions with DoLP also show higher AFE underestimations to the direct method for Cases I and II, i.e., PE = 38% and 37% (D0P+L
and D1P+L, respectively) in Case I, PE = 36% and 34% (D0P+L and D1P+L, respectively) in Case II, and PE = 17% (both) in Case III, but the DoLP addition keeps the ARF closer to that estimated by the direct method, mainly for Case I (high AOD).

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Radiative-flux comparison between AERONET and CM-21 as a function of the SZA from 2009 to 2023. N is the number of observations and the vertical bars are the standard deviations.</figcaption></figure>

<table><caption>**Table 6.** Radiative fluxes, AFE, and ARF comparisons for three case studies.</caption><thead><tr><th rowspan="2">Observations and Inversions</th><th colspan="3">Case I: AOD = 0.30 and SZA = 60°</th><th colspan="3">Case II: AOD = 0.25 and SZA = 60°</th><th colspan="3">Case III: AOD = 0.21 and SZA = 70°</th></tr><tr><th>Flux</th><th>AFE</th><th>ARF</th><th>Flux</th><th>AFE</th><th>ARF</th><th>Flux</th><th>AFE</th><th>ARF</th></tr></thead><tbody><tr><td>CM-21 <sup>1</sup></td><td>440</td><td>-260 ± 22</td><td>-34 ± 26</td><td>444</td><td>-260 ± 22</td><td>-34 ± 26</td><td>306</td><td>-202 ± 18</td><td>-28 ± 21</td></tr><tr><td>AERONET</td><td>445</td><td>-192</td><td>-39</td><td>456</td><td>-175</td><td>-35</td><td>301</td><td>-200</td><td>-30</td></tr><tr><td>D0</td><td>474</td><td>-187</td><td>-38</td><td>480</td><td>-166</td><td>-33</td><td>309</td><td>-187</td><td>-27</td></tr><tr><td>D0+L</td><td>472</td><td>-179</td><td>-39</td><td>476</td><td>-165</td><td>-36</td><td>310</td><td>-179</td><td>-27</td></tr><tr><td>D0P+L</td><td>474</td><td>-161</td><td>-37</td><td>476</td><td>-167</td><td>-36</td><td>308</td><td>-167</td><td>-27</td></tr><tr><td>D1+L</td><td>471</td><td>-177</td><td>-38</td><td>475</td><td>-164</td><td>-34</td><td>307</td><td>-173</td><td>-26</td></tr><tr><td>D1P+L</td><td>473</td><td>-165</td><td>-35</td><td>473</td><td>-171</td><td>-36</td><td>307</td><td>-167</td><td>-26</td></tr></tbody></table>

When the radiative fluxes, ARF, and AFE from GRASP combinations are compared with those from AERONET inversions, (i) the flux overestimations by GRASP inversions decrease the PE to ~6% and ~4% in Cases I and II, respectively; (ii) the AFE underestimations also decrease PE to ~10%, ~5%, and ~13% for Cases I, II, and III, respectively; (iii) the ARF varies according to GRASP combination, under- or overestimating the AERONET inversions (Cases I and II) and PE remains around 10% in Case III; and (iv) the GRASP inversions with DoLP also show AFE underestimations to the AERONET inversions for Cases I and III as well, i.e., PE = 16% and 14% (D0P+L and D1P+L, respectively) in Case I, PE = 17% (both) in Case III, but the DoLP addition keeps the ARF inversions closer to that estimated by AERONET inversions.

# 5. Discussion

## 5.1. AOD and Flux Trends

The downtrends in Figure 3 may be related to the air quality policies performed at Barcelona, hence, the drop in particulate matter (PM1, 2.5, and 10). However, an opposite relation was found between dust transportation and the annual AOD drop, given that dust

<footnote n="1">¹ The AFE and ARF are from 2022, the year of the case studies (Tables 3 and 4, respectively).</footnote>
frequency has increased through the years [83]. Therefore, the local aerosol sources have a dominant impact on the urban-environment AOD compared with dust events.

The first policy, started in 2007, aimed to reduce the speed limit to 30 kmh⁻¹ in 16 zones of Barcelona (https://www.barcelona.cat/mobilitat/es/barcelona-ciudad%E2%E8%88%9230, accessed on 8 August 2024) followed by the Barcelona Urban Mobility Plans of 2006–2012 (http://hdl.handle.net/11703/109492, accessed on 31 July 2024), 2013–2018 (http://hdl.handle.net/11703/85163, accessed on 31 July 2024), and 2019–2024 (http://hdl.handle.net/11703/115353, accessed on 31 July 2024) to reduce the environmental impacts and the contribution of mobility to climate change. The Low Emission Zone policy was then implemented in December 2020 (https://ajuntament.barcelona.cat/qualitataire, accessed on 8 August 2024) to protect areas where vehicles could not circulate without the environmental label from the General Directorate of Traffic.

Consequently, the AOD drop can be connected to the drop in PM1, 2.5, and 10 at Barcelona, improving the local air quality, which is supported by the following studies. Ref. [84] showed decreases in PM2.5 and 10 by applying the MK test and multi-exponential fit for the 2004–2014 trends in the urban background of Barcelona. They attributed the PM decreases to the effectiveness of pollution control measures implemented at European or regional/local levels. Ref. [85] also found a significant decrease in PM2.5 concentrations related to the success of mitigation strategies by reducing the anthropogenic source contributions (heavy oil combustion, mixed combustion, industry, and secondary sulfate). Ref. [86] quantified the composition of PM2.5 and 10 and found contributions of dust (8% and 26%, respectively), marine (<1% and 4%, respectively), and anthropogenic (73% and 54%, respectively). Ref. [87] also showed decreases in the PM1 trend (MK test) of the period from 2012 to 2018 in the urban background, which is related to the reduction of pollutant sources, the oxidation of the secondary precursors, and local meteorological phenomena. Additionally, the 2000–2019 trends of PM2.5 and 10 [88] for the western part of the European Monitoring and Evaluation Programme domain (including Barcelona) dropped for the model and observations. They also attributed the trend drops to the abatement strategies of anthropogenic gaseous precursors of secondary aerosols. Thus, the drop in the AOD could be related to the efficiency of the air quality policies in reducing aerosol concentration in Barcelona in recent decades.

The radiative fluxes (Figure 4) reaching the ground have increased over the years in Barcelona, which could be related to the AOD drop (2009–2023 trend in Figure 3) due to less radiation attenuation by aerosol effects in the last 14 years and due to the scattering caused by clouds, which interferes with global radiation at the surface.

The increase in radiative fluxes (Figure 5) is associated with a decrease in the AOD for winter and autumn from 2009 to 2023, i.e., the solar radiation that reaches the ground is getting less attenuated over the years. The AOD decay may be related to the policies implemented since 2007 by Barcelona's government. Furthermore, the radiative flux decrease is associated with an AOD increase for spring and summer due to the presence of aerosol from local sources (pollens) and from external sources (desert dust and biomass burning), which are common in these seasons [27,29–31]. The local pollution reduction policies cannot control those aerosol types.

Regarding the diurnal cycles, the all-skies maximums (Figure 6) exceed 1000 Wm⁻² in winter, reach 1200 Wm⁻² in autumn, and exceed 1200 Wm⁻² in spring and summer due to the multiple scattering produced by fractional cloud cover [89–92]. The presence of pollen, desert dust outbreaks, biomass burning from forest fires, and European pollution episodes [27,29–31] is more common in spring and summer, scattering or absorbing the radiation that reaches the ground. Furthermore, aerosol fluxes are greater in summer [11] with a factor of 1.5 compared to winter.
## 5.2. AFE and ARF in 14 Years

The results presented in Figure 7 follow [27], who found high efficiencies at high SZAs, -132.2 Wm⁻²τ⁻¹ (50° < SZA < 60°) in Mallorca Island, Spain; ref. [20], who obtained -123 ± 11 Wm⁻²τ⁻¹ at 440 nm (SZA < 65°); and ref. [93], who derived daily aerosol efficiencies varying from -188 ± 18 to -163 ± 16 Wm⁻²τ⁻¹ (SZAs from 32° to 73°) in Lampedusa Island, Italy.

The high standard deviations (Table 4) may be related to the direct method (Section 3.3), which directly depends on the aerosol load (high variability throughout the year for all aerosol types). This leads to an opposite angular dependency related to that in Table 3: at higher solar angles (thicker atmosphere), there is more AOD variability, increasing their standard deviations.

Higher absolute values of AFE (first plot of Figure 8) are found at intermediate SZAs because the atmosphere is optically thicker than at smaller SZAs, and when the solar radiation path increases in the atmosphere, the attenuation increases [94], especially at shorter wavelengths, and then the AFE decreases [25]. The AFE drop at the highest SZA is associated with the decrease in the aerosol-scattered flux because the slant path is no longer optically thin [95]. Ref. [21] associated the ARF decay at 70° with the slant path no longer being optically thin and [77] attributed the decay at larger SZAs to the upscatter flux at sunset and sunrise being larger in comparison to the one at local noon.

Regarding the contribution of the other aerosol properties to the ARF, the AERONET inversions do not allow a robust analysis of the SSA and other parameters (e.g., absorption AOD and asymmetry factor) for the corresponding clear-sky days, owing to some data being available at level 2. The relevance of the other aerosol parameters is an open issue.

Additionally, according to reports of the Ministry for Ecological Transition and the Demographic Challenge (https://www.miteco.gob.es/es/calidad-y-evaluacion-ambiental/temas/atmosfera-y-calidad-del-aire/calidad-del-aire/evaluacion-datos/fuentes-naturales/anuales.html, accessed on 3 April 2025) and ref. [83], dust episodes have increased by 5% to 15% from 2009 to 2023 in Spain and in the western Mediterranean for the 1948–2020 period, respectively. Table A8 shows the number of days by year with desert dust has increased. This fact contributes to the increase in the absolute values of the aerosol forcings (cooling effect) due to dust upscattering, mainly from 2019 (19 days).

## 5.3. AFE and ARF from Aerosol Types

The AFE values for DD (Figure 11 and Table 5) are close to -120 Wm⁻²τ⁻¹ (440 nm at 53–75°) found by [18] during a strong dust event mixed with smoke in Granada; close to those found by [45], -185 to -81.7 Wm⁻²τ⁻¹ (495.7 nm at 20–75°), in the Mediterranean island of Lampedusa; and follow those found by [22] for two case studies (-124.4 and -149.1 Wm⁻²τ⁻¹ for the first case and -169.4, -178.9, and -198.9 Wm⁻²τ⁻¹ for the second case) within a dust transport episode in Granada, Spain.

The high MA variability (Figure 11 and Table 5) can be related to the aerosol mixture that contains a large variability of single scattering albedo [25], bearing in mind that DD has the largest aerosol load (0.31 ± 0.09) and that MA is the dominant aerosol (60.9%) with the largest AOD variability (0.13 ± 0.08), which can explain their angular dependence and their CI95s. The AFE dependence on SZA for UI-BB changes according to the aerosol types, the optical thickness of the atmosphere, and its attenuation [25,45].

The ARF for DD (Figure 11 and Table 5) are in the range of those estimated by [26] for instantaneous radiative forcing at 500 nm, varying from -93.1 to -0.5 Wm⁻² during dust outbreaks in Barcelona. The DD has the largest ARF because of the largest AOD at short wavelengths.
### 5.4. Forcing Comparisons

Regarding the flux comparisons, a similar overestimation of the pyranometer by AERONET fluxes (4 ± 5.2%) was also found by [73], who attributed it to the combination of errors, such as cosine response, calibration, and linearity, in the pyranometer measurements. The overestimations of the CM-21 fluxes by GRASP ones can be relative to the cosine response (greater SZAs, smaller PE) and due to the spectral range in which the GRASP code integrates the broadband solar flux from 200 to 4000 nm [64] and the CM-21 from 305 to 2800 nm (Table 2).

As the ARF depends directly on the AOD, any AOD variation will change the ARF observed when DoLP is added to the GRASP synergies and compared with the direct method. Refs. [32,71] pointed to adjustments in the higher AOD when they were inverted with DoLP. The small over- and underestimations of the probably owing to the small AOD variation. Furthermore, the differences between the AERONET and GRASP inversions can be attributed to the AFE-ARF definitions (Appendix B and Section 3.3). According to the AERONET DOCUMENT (n.d.) and [73], the ARF is the difference between the broadband fluxes with and without aerosols at the BOA, and the AFE is the rate at which the atmosphere is forced per unit of the AOD at 550 nm (Appendix B). However, the GRASP code retrieves the aerosol net radiative effect at the BOA, which is the difference between downwelling and upwelling fluxes at a given atmospheric layer in aerosol-free and aerosol-laden conditions [64]. The formulas are in Appendix B.

## 6. Conclusions

This research presented the long-term basis aerosol radiative effects over a decade in Barcelona, Spain, for the first time. The ARF and AFE were computed by the direct method and by combining radiation measurements from the CM-21 pyranometer (L1.5) and AOD from the AERONET photometer (L2). The AFE was derived from the slope between net fluxes and AOD at 440, 675, 870, and 1020 nm, and the ARF was calculated for six SZA groups (20°, 30°, 40°, 50°, 60°, and 70°). Clear-sky conditions were selected from all-skies conditions by a quadratic fitting with seasonal adaptative criteria. Moreover, the aerosol was characterized and classified to investigate the AFE-ARF contributions from each aerosol type.

The all-skies trend of the AOD showed a decrease in the last decade due to the local government's policies to reduce air pollution in Barcelona, whereas the flux trend increased for the same period. Additionally, the clear-sky seasonality showed decays in the flux trends followed by increases in the AOD trends and vice versa. The AOD increases in spring and summer were expected because of the contributions of other aerosol types (dust and pollen, for example). All trends highlighted the influence of the aerosol load on global solar radiation that reaches the surface.

The AFE-ARF long-term revealed a cooling effect, which increases in absolute number over 14 years, 23.9% and 40.2% for AFE and ARF, respectively. The efficiencies varied from -331 and -10 Wm⁻²τ⁻¹ with an AOD varying from 0.16 to 0.69 and the ARF from -64 to -2 Wm⁻². Furthermore, the AFE-ARF angular dependency increased in the absolute number at centered SZA (40° to 60°) because of the attenuation in a thicker atmosphere.

This research showed that a CI95 of the AFE depended mostly on the SZA where the smaller angles have fewer observations, hence, a higher CI95. The seasonal variability of the CI95 was also considered. The years with fewer clear-sky days (2015, 2016, and 2017) also contributed to a high CI95. On the other hand, the ARF standard deviations were directly influenced by the aerosol load variability due to the direct method.

Regarding the aerosol types, Barcelona had three main dominant aerosols for clear-sky days: MA (60.9%), UI-BB (28.9%), and DD (10.2%). This classification combined with
Gobbi's method clustered the aerosols into four groups by AE analysis, two coarse modes and two fine modes. Finally, the contribution of the aerosol types to the forcing concluded that DD forcing has had the greatest cooling effect in Barcelona, varying from -62 to -37 Wm⁻², followed by UI-BB (-40 to -20 Wm⁻²) and MA (-32 to -24 Wm⁻²).

Comparisons between CM-21 observations and AERONET and GRASP inversions presented similar broadband-flux measurements for 958 inversions and three case studies. The AFE-ARF estimations had noticeable differences when the direct method was compared with GRASP retrievals; however, these differences decreased when GRASP outputs were compared with AERONET inversions, mostly for Cases I and II. The DoLP addition in GRASP synergy (D0P+L and D1P+L combinations) impacted the AFE by raising the underestimation of the direct method (38% with D0P+L and 37% with D1P+L in Case I, and 17% with both combinations in Case III) and reducing ARF overestimations in Case I (high AOD). Therefore, the GRASP code could retrieve the ARF satisfactorily, owing to the approximation between its ARF values and those from the direct method.

More complementary studies need to be performed in Barcelona to improve the understanding of aerosol-forcing contributions to the radiation budget, for instance, (i) investigating the influence of water vapor on the radiation budget [44,93,96] and the aerosol forcing; (ii) analyzing surface albedo and other aerosol properties relevant to estimate their radiative effects, as single scattering albedo, asymmetry factor, phase function, and fine-mode fraction [97,98] which depend on the composition, size distribution, and shape of the particles, varying with wavelength and height [10]; (iii) characterizing the chemical speciation of the aerosol types; (iv) validating the radiative fluxes from GRASP code by retrieving aerosol radiative properties from a large database; and (v) validating the forcing by satellite observations.

**Author Contributions:** Conceptualization, D.C.F.d.S.O., M.S. and A.R.-G.; Methodology, D.C.F.d.S.O., M.S., A.R.-G., C.M.-P., A.C., C.G.-D. and Y.D.; Writing—Original Draft Preparation, D.C.F.d.S.O.; Writing—Review and Editing D.C.F.d.S.O., M.S., A.R.-G., C.M.-P., A.C., C.G.-D., O.D., A.L., M.M. and Y.D.; Software, O.D., A.L., M.M. and Y.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by (i) the Spanish Ministry of Science and Innovation and Universities MCIN/AEI/10.13039/501100011033: Red de Investigación RED2022-134824-E RED ESPAÑOLA DE ENTIDADES PARTICIPANTES EN ACTRIS; PID2023-149747NB-I00, Synergetic passive and active remote sensing of aerosol vertical profiling with applications to aerosol typing, cloud characterization and energy budget (SYRSAC); Support for membership to ACTRIS-ERIC; and (ii) European commission: H2020—Research and Innovation Framework Programme (ATMO-ACCESS, grant no. 10100800, ACTRIS IMP, grant no. 871115, GRASP ACE, grant no. 778349) and Horizon Europe Framework Programme (REALISTIC, grant no. 101086690). M.S. also received funding from CNES through EECLAT, AOS, and EXTRA-SAT projects.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. This data can be found here: https://aeronet.gsfc.nasa.gov (accessed on 4 March 2025), https://solrad-net.gsfc.nasa.gov (accessed on 4 March 2025), and https://www.earlinet.org/index.php?id=328 (accessed on 4 March 2025).

**Acknowledgments:** The authors are especially grateful to the Laboratoire d'Optique Atmosphérique at the Université de Lille and the GRASP-SAS company. The authors also acknowledge the use of GRASP-processed data (https://www.grasp-open.com (accessed on 4 March 2025)) in this work.

**Conflicts of Interest:** Oleg Dubovik supported GRASP SAS with scientific consulting and held the actions of the company under "concours scientifique" convention of CNRS. Furthermore, authors Masahiro Momoi and Anton Lopatin were employed by the company GRASP-SAS. The remaining
authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Appendix A

The AFE values at 675, 870, and 1020 nm are plotted in Figures A1–A3 and summarized in Tables A1, A3 and A5, consecutively. The ARFs are summarized in Tables A2, A4 and A6 at the same wavelengths. All AFE values greater than $600 \text{ Wm}^{-2}\tau^{-1}$ (absolute value) were eliminated from this research to avoid efficiencies that are not documented in the literature. The AOD spectral dependence is propagated to the AFE-ARF spectral dependence because the direct method uses the AOD for estimating the forcing (Section 3.3). According to [44], the larger fraction of diffuse radiation concerning the total is found at larger SZAs, and when the SZA increases, radiation at the shorter wavelengths is reduced more strongly than at long wavelengths.

<figure><img src="image_15.png" /><figcaption>**Figure A1.** The AFE long-term from 2009 to 2023 for six SZA groups at 675 nm. The linear regressions (LR) slopes were calculated with a CI95.</figcaption></figure>

The AFE-ARF time series at 675 nm and 20° vary from -343 to -118 $\text{Wm}^{-2}\tau^{-1}$ (Table A1) and from -35 to -10 $\text{Wm}^{-2}$ (Table A4), respectively. The AOD varies from 0.04 to 0.63. These results follow those found by [15], whose integrated radiative forcing was 43 $\text{Wm}^{-2}$ at 670 nm and 20° in Almería, Spain; by [16], who obtained an efficiency of -73 $\text{Wm}^{-2}\tau^{-1}$ and forcing of -30 $\text{Wm}^{-2}$ at 670 nm and 20° in Granada during the 2003 heat wave; and by [5], who calculated efficiency of 279 ± 21 $\text{Wm}^{-2}\tau^{-1}$ and radiative forcing of 23 ± 7 $\text{Wm}^{-2}$ at 675 nm and 15° in Granada, Spain. The AFE and ARF vary from
−524 to −11 Wm⁻²τ⁻¹ and −53 to −2 Wm⁻², respectively, in the entire SZA range. These results are similar to those found by [25], where the AFE varied from −450 to −9 Wm⁻²τ⁻¹ and the ARF varied from −40 to −1 Wm⁻² at 675 nm in Granada. Ref. [20] also calculated an efficiency of 115 ± 11 Wm⁻²τ⁻¹ at 675 nm and SZA < 65° in Granada.

The AFE and ARF at 870 nm (Tables A3 and A4, respectively) have variations from −560 to −39 Wm⁻²τ⁻¹ and from −47 to −5 Wm⁻² associated with AOD that vary from 0.007 to 0.61. These values are similar to those found by [20,45], the first one obtained 116 ± 10 Wm⁻²τ⁻¹ and the second one obtained daily efficiencies varying from −116.6 to −56.0 Wm⁻²τ⁻¹ at 868.7 nm and an SZA of 60°, even though the efficiency values have a large variability over the years for this study. Lastly, the AFE and ARF at 1020 nm (Tables A5 and A6, respectively) show variations from −599 to −58 Wm⁻²τ⁻¹ and from −44 to −7 Wm⁻², with AOD varying from 0.003 to 0.60.

<table><caption>**Table A1.** The AFE on the surface is a function of the AOD at 675 nm for six SZA groups. The slope was estimated with a CI95. Y stands for years.</caption><thead><tr><th></th><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr><tr><th></th><th>Y</th><th colspan="6">AFE, Wm<sup>-2</sup> τ<sup>-1</sup></th></tr></thead><tbody><tr><td>2009</td><td></td><td>-295 ± 55</td><td>-323 ± 63</td><td>-355 ± 88</td><td>-274 ± 82</td><td>-300 ± 83</td><td>-343 ± 120</td></tr><tr><td>2010</td><td></td><td>-316 ± 47</td><td>-409 ± 69</td><td>-288 ± 112</td><td>-464 ± 109</td><td>-251 ± 89</td><td>-266 ± 58</td></tr><tr><td>2011</td><td></td><td>-305 ± 71</td><td>-372 ± 101</td><td>-302 ± 91</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td></td><td>-256 ± 47</td><td>-261 ± 75</td><td>-182 ± 87</td><td>-332 ± 86</td><td>-</td><td>-</td></tr><tr><td>2013</td><td></td><td>-411 ± 93</td><td>-</td><td>-524 ± 170</td><td>-330 ± 101</td><td>-203 ± 165</td><td>-216 ± 62</td></tr><tr><td>2014</td><td></td><td>-407 ± 70</td><td>-420 ± 80</td><td>-172 ± 76</td><td>-134 ± 94</td><td>-131 ± 130</td><td>-132 ± 85</td></tr><tr><td>2015</td><td></td><td>-255 ± 60</td><td>-93 ± 68</td><td>-186 ± 47</td><td>-156 ± 59</td><td>-282 ± 77</td><td>-175 ± 40</td></tr><tr><td>2016</td><td></td><td>-163 ± 73</td><td>-190 ± 86</td><td>-210 ± 104</td><td>-210 ± 79</td><td>-233 ± 62</td><td>-</td></tr><tr><td>2017</td><td></td><td>-204 ± 146</td><td>-148 ± 165</td><td>-194 ± 154</td><td>-331 ± 287</td><td>-371 ± 1845</td><td>-</td></tr><tr><td>2018</td><td></td><td>-172 ± 51</td><td>-195 ± 75</td><td>-306 ± 73</td><td>-134 ± 70</td><td>-</td><td>-137 ± 64</td></tr><tr><td>2019</td><td></td><td>-242 ± 35</td><td>-273 ± 52</td><td>-294 ± 66</td><td>-450 ± 63</td><td>-318 ± 85</td><td>-199 ± 46</td></tr><tr><td>2021</td><td></td><td>-306 ± 41</td><td>-341 ± 77</td><td>-319 ± 111</td><td>-253 ± 71</td><td>-11 ± 117</td><td>-119 ± 126</td></tr><tr><td>2022</td><td></td><td>-277 ± 30</td><td>-331 ± 37</td><td>-325 ± 102</td><td>-241 ± 56</td><td>-201 ± 26</td><td>-166 ± 15</td></tr><tr><td>2023</td><td></td><td>-298 ± 45</td><td>-406 ± 67</td><td>-344 ± 143</td><td>-275 ± 95</td><td>-232 ± 82</td><td>-172 ± 53</td></tr></tbody></table>

<table><caption>**Table A2.** The ARF means on the surface and their standard deviation at 675 nm for six SZA groups. Y stands for years.</caption><thead><tr><th></th><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr><tr><th></th><th>Y</th><th colspan="6">ARF, Wm<sup>-2</sup></th></tr></thead><tbody><tr><td>2009</td><td></td><td>-33 ± 16</td><td>-41 ± 19</td><td>-41 ± 19</td><td>-34 ± 13</td><td>-38 ± 15</td><td>-35 ± 10</td></tr><tr><td>2010</td><td></td><td>-27 ± 19</td><td>-34 ± 22</td><td>-31 ± 15</td><td>-43 ± 26</td><td>-26 ± 12</td><td>-30 ± 14</td></tr><tr><td>2011</td><td></td><td>-20 ± 11</td><td>-26 ± 13</td><td>-22 ± 12</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td></td><td>-17 ± 13</td><td>-20 ± 13</td><td>-23 ± 9</td><td>-44 ± 17</td><td>-</td><td>-</td></tr><tr><td>2013</td><td></td><td>-25 ± 16</td><td>-</td><td>-42 ± 25</td><td>-27 ± 16</td><td>-17 ± 9</td><td>-24 ± 7</td></tr><tr><td>2014</td><td></td><td>-27 ± 24</td><td>-29 ± 25</td><td>-17 ± 11</td><td>-13 ± 10</td><td>-12 ± 8</td><td>-10 ± 5</td></tr><tr><td>2015</td><td></td><td>-29 ± 17</td><td>-12 ± 6</td><td>-22 ± 12</td><td>-20 ± 11</td><td>-28 ± 14</td><td>-22 ± 11</td></tr><tr><td>2016</td><td></td><td>-11 ± 7</td><td>-12 ± 8</td><td>-11 ± 7</td><td>-13 ± 9</td><td>-12 ± 11</td><td>-</td></tr><tr><td>2017</td><td></td><td>-16 ± 5</td><td>-11 ± 4</td><td>-16 ± 7</td><td>-24 ± 7</td><td>-24 ± 3</td><td>-</td></tr><tr><td>2018</td><td></td><td>-19 ± 10</td><td>-21 ± 12</td><td>-34 ± 18</td><td>-19 ± 8</td><td>-</td><td>-16 ± 9</td></tr><tr><td>2019</td><td></td><td>-31 ± 18</td><td>-34 ± 20</td><td>-35 ± 20</td><td>-53 ± 26</td><td>-47 ± 23</td><td>-28 ± 15</td></tr><tr><td>2021</td><td></td><td>-27 ± 20</td><td>-32 ± 22</td><td>-39 ± 19</td><td>-36 ± 17</td><td>-2 ± 1</td><td>-18 ± 7</td></tr><tr><td>2022</td><td></td><td>-23 ± 19</td><td>-25 ± 23</td><td>-34 ± 18</td><td>-27 ± 17</td><td>-23 ± 17</td><td>-21 ± 19</td></tr><tr><td>2023</td><td></td><td>-21 ± 20</td><td>-33 ± 29</td><td>-33 ± 22</td><td>-37 ± 30</td><td>-36 ± 27</td><td>-33 ± 21</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Figure A2.** The AFE long-term from 2009 to 2023 for six SZA groups at 870 nm. The linear regressions (LR) slopes were calculated with a CI95.</figcaption></figure>

<table><caption>**Table A3.** The AFE on the surface is a function of the AOD at 870 nm for six SZA groups. The slope was estimated with a CI95. Y stands for years.</caption><thead><tr><th rowspan="2">Y</th><th colspan="6">AFE, Wm<sup>-2</sup> τ<sup>-1</sup></th></tr><tr><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr></thead><tbody><tr><td>2009</td><td>-426 ± 76</td><td>-448 ± 86</td><td>-503 ± 122</td><td>-366 ± 103</td><td>-417 ± 103</td><td>-437 ± 118</td></tr><tr><td>2010</td><td>-392 ± 61</td><td>-527 ± 90</td><td>-377 ± 139</td><td>-560 ± 144</td><td>-267 ± 115</td><td>-271 ± 81</td></tr><tr><td>2011</td><td>-406 ± 104</td><td>-508 ± 147</td><td>-431 ± 154</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-371 ± 67</td><td>-375 ± 109</td><td>-261 ± 140</td><td>-458 ± 184</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-534 ± 124</td><td>-</td><td>-</td><td>-420 ± 155</td><td>-270 ± 202</td><td>-407 ± 113</td></tr><tr><td>2014</td><td>-505 ± 89</td><td>-527 ± 112</td><td>-205 ± 103</td><td>-164 ± 128</td><td>-168 ± 180</td><td>-127 ± 111</td></tr><tr><td>2015</td><td>-252 ± 74</td><td>-84 ± 74</td><td>-186 ± 61</td><td>-154 ± 75</td><td>-371 ± 120</td><td>-179 ± 48</td></tr><tr><td>2016</td><td>-200 ± 89</td><td>-229 ± 98</td><td>-261 ± 130</td><td>-236 ± 93</td><td>-259 ± 73</td><td>-</td></tr><tr><td>2017</td><td>-321 ± 209</td><td>-217 ± 244</td><td>-304 ± 229</td><td>-347 ± 390</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-179 ± 66</td><td>-183 ± 95</td><td>-314 ± 103</td><td>-124 ± 82</td><td>-</td><td>-123 ± 71</td></tr><tr><td>2019</td><td>-263 ± 41</td><td>-288 ± 62</td><td>-304 ± 77</td><td>-489 ± 83</td><td>-332 ± 106</td><td>-194 ± 52</td></tr><tr><td>2021</td><td>-324 ± 50</td><td>-342 ± 95</td><td>-314 ± 127</td><td>-237 ± 80</td><td>-39 ± 122</td><td>-109 ± 129</td></tr><tr><td>2022</td><td>-297 ± 36</td><td>-344 ± 45</td><td>-378 ± 123</td><td>-233 ± 65</td><td>-199 ± 29</td><td>-161 ± 17</td></tr><tr><td>2023</td><td>-316 ± 51</td><td>-440 ± 79</td><td>-381 ± 155</td><td>-290 ± 104</td><td>-233 ± 94</td><td>-175 ± 58</td></tr></tbody></table>
<table><caption>**Table A4.** The ARF means on the surface and their standard deviation at 870 nm for six SZA groups. Y stands for years.</caption><thead><tr><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr><tr><th>Y</th><th colspan="6">ARF, Wm<sup>-2</sup></th></tr></thead><tbody><tr><td>2009</td><td>-34 ± 16</td><td>-42 ± 19</td><td>-44 ± 20</td><td>-35 ± 13</td><td>-41 ± 16</td><td>-35 ± 11</td></tr><tr><td>2010</td><td>-25 ± 18</td><td>-32 ± 22</td><td>-30 ± 16</td><td>-39 ± 25</td><td>-21 ± 11</td><td>-23 ± 12</td></tr><tr><td>2011</td><td>-19 ± 10</td><td>-26 ± 12</td><td>-23 ± 11</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-17 ± 13</td><td>-21 ± 13</td><td>-24 ± 8</td><td>-46 ± 15</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-24 ± 16</td><td>-</td><td>-</td><td>-26 ± 14</td><td>-17 ± 10</td><td>-30 ± 7</td></tr><tr><td>2014</td><td>-24 ± 24</td><td>-27 ± 24</td><td>-15 ± 10</td><td>-11 ± 9</td><td>-11 ± 7</td><td>-7 ± 4</td></tr><tr><td>2015</td><td>-21 ± 15</td><td>-8 ± 5</td><td>-16 ± 11</td><td>-15 ± 10</td><td>-26 ± 13</td><td>-17 ± 10</td></tr><tr><td>2016</td><td>-11 ± 7</td><td>-12 ± 8</td><td>-12 ± 7</td><td>-13 ± 9</td><td>-12 ± 10</td><td>-</td></tr><tr><td>2017</td><td>-18 ± 6</td><td>-12 ± 4</td><td>-18 ± 8</td><td>-18 ± 5</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-15 ± 9</td><td>-16 ± 10</td><td>-27 ± 16</td><td>-13 ± 7</td><td>-</td><td>-12 ± 8</td></tr><tr><td>2019</td><td>-27 ± 17</td><td>-29 ± 19</td><td>-30 ± 19</td><td>-47 ± 25</td><td>-40 ± 22</td><td>-24 ± 14</td></tr><tr><td>2021</td><td>-23 ± 19</td><td>-24 ± 19</td><td>-30 ± 18</td><td>-27 ± 16</td><td>-5 ± 2</td><td>-13 ± 6</td></tr><tr><td>2022</td><td>-19 ± 18</td><td>-21 ± 21</td><td>-31 ± 17</td><td>-21 ± 15</td><td>-17 ± 16</td><td>-17 ± 18</td></tr><tr><td>2023</td><td>-17 ± 19</td><td>-27 ± 29</td><td>-28 ± 22</td><td>-34 ± 30</td><td>-30 ± 25</td><td>-29 ± 20</td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Figure A3.** The AFE long-term from 2009 to 2023 for six SZA groups at 1020 nm. The linear regressions (LR) slopes were calculated with a CI95.</figcaption></figure>
<table><caption>**Table A5.** The AFE on the surface is a function of the AOD at 1020 nm for six SZA groups. The slope was estimated with a CI95. Y stands for years.</caption><thead><tr><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr><tr><th>Y</th><th colspan="6">AFE, Wm<sup>-2</sup> τ<sup>-1</sup></th></tr></thead><tbody><tr><td>2009</td><td>-512 ± 93</td><td>-513 ± 104</td><td>-585 ± 150</td><td>-421 ± 116</td><td>-487 ± 118</td><td>-461 ± 122</td></tr><tr><td>2010</td><td>-426 ± 70</td><td>-580 ± 104</td><td>-418 ± 154</td><td>-599 ± 166</td><td>-268 ± 129</td><td>-263 ± 93</td></tr><tr><td>2011</td><td>-476 ± 118</td><td>-592 ± 167</td><td>-525 ± 192</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-484 ± 83</td><td>-493 ± 128</td><td>-343 ± 183</td><td>-475 ± 251</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-589 ± 143</td><td>-</td><td>-</td><td>-453 ± 190</td><td>-303 ± 220</td><td>-588 ± 177</td></tr><tr><td>2014</td><td>-548 ± 100</td><td>-567 ± 132</td><td>-215 ± 117</td><td>-181 ± 148</td><td>-175 ± 202</td><td>-125 ± 121</td></tr><tr><td>2015</td><td>-249 ± 82</td><td>-81 ± 77</td><td>-186 ± 68</td><td>-151 ± 82</td><td>-412 ± 149</td><td>-181 ± 52</td></tr><tr><td>2016</td><td>-218 ± 98</td><td>-248 ± 104</td><td>-277 ± 149</td><td>-248 ± 101</td><td>-273 ± 80</td><td>-</td></tr><tr><td>2017</td><td>-382 ± 253</td><td>-296 ± 303</td><td>-442 ± 285</td><td>-408 ± 444</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-172 ± 74</td><td>-172 ± 104</td><td>-299 ± 119</td><td>-118 ± 87</td><td>-</td><td>-120 ± 74</td></tr><tr><td>2019</td><td>-273 ± 44</td><td>-292 ± 66</td><td>-315 ± 81</td><td>-513 ± 90</td><td>-338 ± 115</td><td>-193 ± 55</td></tr><tr><td>2021</td><td>-338 ± 55</td><td>-339 ± 104</td><td>-314 ± 134</td><td>-228 ± 85</td><td>-58 ± 126</td><td>-108 ± 131</td></tr><tr><td>2022</td><td>-305 ± 40</td><td>-348 ± 49</td><td>-392 ± 138</td><td>-224 ± 70</td><td>-197 ± 31</td><td>-160 ± 19</td></tr><tr><td>2023</td><td>-326 ± 55</td><td>-459 ± 81</td><td>-392 ± 164</td><td>-300 ± 111</td><td>-235 ± 101</td><td>-178 ± 62</td></tr></tbody></table>

<table><caption>**Table A6.** The ARF means on the surface and their standard deviation at 1020 nm for six SZA groups. Y stands for years.</caption><thead><tr><th>SZA</th><th>70 ± 1°</th><th>60 ± 1°</th><th>50 ± 1°</th><th>40 ± 1°</th><th>30 ± 1°</th><th>20 ± 1°</th></tr><tr><th>Y</th><th colspan="6">ARF, Wm<sup>-2</sup></th></tr></thead><tbody><tr><td>2009</td><td>-34 ± 16</td><td>-40 ± 18</td><td>-43 ± 19</td><td>-34 ± 13</td><td>-41 ± 16</td><td>-32 ± 11</td></tr><tr><td>2010</td><td>-23 ± 18</td><td>-30 ± 22</td><td>-28 ± 16</td><td>-36 ± 23</td><td>-18 ± 10</td><td>-19 ± 11</td></tr><tr><td>2011</td><td>-18 ± 11</td><td>-25 ± 13</td><td>-24 ± 11</td><td>-</td><td>-</td><td>-</td></tr><tr><td>2012</td><td>-18 ± 13</td><td>-21 ± 14</td><td>-25 ± 8</td><td>-39 ± 13</td><td>-</td><td>-</td></tr><tr><td>2013</td><td>-24 ± 15</td><td>-</td><td>-</td><td>-26 ± 13</td><td>-18 ± 10</td><td>-38 ± 7</td></tr><tr><td>2014</td><td>-23 ± 24</td><td>-26 ± 23</td><td>-13 ± 9</td><td>-11 ± 9</td><td>-10 ± 7</td><td>-7 ± 4</td></tr><tr><td>2015</td><td>-18 ± 14</td><td>-7 ± 5</td><td>-14 ± 10</td><td>-13 ± 9</td><td>-24 ± 13</td><td>-15 ± 10</td></tr><tr><td>2016</td><td>-10 ± 7</td><td>-11 ± 8</td><td>-11 ± 7</td><td>-11 ± 9</td><td>-11 ± 10</td><td>-</td></tr><tr><td>2017</td><td>-17 ± 6</td><td>-13 ± 5</td><td>-20 ± 9</td><td>-17 ± 6</td><td>-</td><td>-</td></tr><tr><td>2018</td><td>-12 ± 8</td><td>-12 ± 9</td><td>-21 ± 14</td><td>-11 ± 6</td><td>-</td><td>-10 ± 7</td></tr><tr><td>2019</td><td>-25 ± 17</td><td>-27 ± 18</td><td>-28 ± 18</td><td>-44 ± 24</td><td>-37 ± 21</td><td>-21 ± 14</td></tr><tr><td>2021</td><td>-21 ± 18</td><td>-21 ± 18</td><td>-26 ± 17</td><td>-23 ± 14</td><td>-7 ± 4</td><td>-12 ± 6</td></tr><tr><td>2022</td><td>-17 ± 17</td><td>-18 ± 20</td><td>-28 ± 16</td><td>-17 ± 13</td><td>-16 ± 15</td><td>-14 ± 18</td></tr><tr><td>2023</td><td>-15 ± 19</td><td>-24 ± 29</td><td>-24 ± 22</td><td>-32 ± 29</td><td>-27 ± 24</td><td>-27 ± 20</td></tr></tbody></table>

Table A7 summarizes the percentages of the CI95s to their AFEs and the percentages of the standard deviations to their ARF means. Table A8 shows the number of clear-sky days and dust days, and sample sizes of net fluxes as a function of the AOD for each SZA and year presented in Figure 7.

<table><caption>**Table A7.** Ratio between CI95s and AFEs and between standard deviations and ARF means by SZA and years at 440 nm. The ratios were calculated from Tables 3 and 4. Y stands for years.</caption><thead><tr><th>SZA</th><th>20 ± 1°</th><th>30 ± 1°</th><th>40 ± 1°</th><th>50 ± 1°</th><th>60 ± 1°</th><th>70 ± 1°</th></tr><tr><th>Y</th><th colspan="6">CL95/AFE (%), Standard Deviation/ARF (%)</th></tr></thead><tbody><tr><td>2009</td><td>57.9, 25.9</td><td>32.7, 38.9</td><td>30.7, 36.1</td><td>23.1, 48.8</td><td>0.2, 46.3</td><td>18.4, 48.5</td></tr><tr><td>2010</td><td>16.5, 39.5</td><td>29.2, 39.4</td><td>20.7, 57.5</td><td>36.4, 45.7</td><td>16.7, 61.1</td><td>14.1, 63.3</td></tr><tr><td>2011</td><td>-</td><td>-</td><td>64.6, 33.3</td><td>27.5, 59.1</td><td>26.8, 50</td><td>21.4, 57.1</td></tr></tbody></table>
<table><caption>Table A7. Cont.</caption><thead><tr><th rowspan="2">Y</th><th colspan="6">CL95/AFE (%), Standard Deviation/ARF (%)</th></tr><tr><th>20 ± 1°</th><th>30 ± 1°</th><th>40 ± 1°</th><th>50 ± 1°</th><th>60 ± 1°</th><th>70 ± 1°</th></tr></thead><tbody><tr><td>2012</td><td>-</td><td>-</td><td>23.4, 51.4</td><td>37.1, 37</td><td>27, 66.7</td><td>17.9, 76.5</td></tr><tr><td>2013</td><td>31.9, 33.3</td><td>111.6, 50</td><td>28.8, 59.3</td><td>31.5, 61</td><td>21.5, 77.1</td><td>21.4, 68</td></tr><tr><td>2014</td><td>43.3, 40</td><td>90.7, 61.5</td><td>61.2, 78.6</td><td>41.8, 61.1</td><td>18.5, 89.3</td><td>18, 85.7</td></tr><tr><td>2015</td><td>23.6, 36.7</td><td>24.3, 46.7</td><td>26, 40</td><td>21.4, 38.7</td><td>51, 34.8</td><td>17.6, 48.8</td></tr><tr><td>2016</td><td>-</td><td>26.7, 68.8</td><td>36.1, 62.5</td><td>53.8, 58.3</td><td>47.5, 53.3</td><td>43, 53.9</td></tr><tr><td>2017</td><td>-</td><td>-</td><td>59.1, 26.5</td><td>73.5, 47.1</td><td>101.3, 41.7</td><td>96.1, 33.3</td></tr><tr><td>2018</td><td>34.8, 37.9</td><td>670, 50</td><td>38.9, 35.7</td><td>18.9, 50</td><td>27, 50</td><td>26.4, 47.8</td></tr><tr><td>2019</td><td>17.9, 40.5</td><td>23.5, 43.9</td><td>13.2, 42.2</td><td>20.7, 43.8</td><td>17.9, 45.5</td><td>15.6, 47.2</td></tr><tr><td>2021</td><td>83.6, 26.7</td><td>-</td><td>21.2, 38.2</td><td>29.6, 38.2</td><td>16.6, 60.5</td><td>10.9, 61.1</td></tr><tr><td>2022</td><td>5.5, 62.9</td><td>10.3, 54.3</td><td>16.2, 47.6</td><td>29.8, 47.4</td><td>8.5, 76.5</td><td>8.9, 75</td></tr><tr><td>2023</td><td>26.5, 48.9</td><td>29, 56</td><td>33.2, 66</td><td>45.2, 51.2</td><td>18, 70</td><td>13.2, 75</td></tr></tbody></table>

<table><caption>Table A8. Number of clear-sky days, including dust days, and sample size of the net fluxes as a function of the AOD, 440 nm, for each SZA group by year, as depicted in Figure 7. Y stands for years.</caption><thead><tr><th rowspan="2">Y</th><th colspan="6">Number of Days, Sample Sizes</th><th rowspan="2">Total <sup>1</sup></th><th rowspan="2">Dust Days</th></tr><tr><th>20 ± 1°</th><th>30 ± 1°</th><th>40 ± 1°</th><th>50 ± 1°</th><th>60 ± 1°</th><th>70 ± 1°</th></tr></thead><tbody><tr><td>2009</td><td>10, 40</td><td>20, 52</td><td>23, 44</td><td>23, 74</td><td>31, 62</td><td>34, 104</td><td>35, 376</td><td>4</td></tr><tr><td>2010</td><td>13, 50</td><td>18, 37</td><td>29, 79</td><td>35, 74</td><td>52, 122</td><td>56, 161</td><td>58, 523</td><td>1</td></tr><tr><td>2011</td><td>-</td><td>-</td><td>4, 10</td><td>23, 69</td><td>32, 103</td><td>50, 148</td><td>51, 330</td><td>3</td></tr><tr><td>2012</td><td>-</td><td>-</td><td>6, 24</td><td>17, 46</td><td>33, 95</td><td>45, 146</td><td>45, 311</td><td>1</td></tr><tr><td>2013</td><td>4, 12</td><td>14, 33</td><td>23, 48</td><td>24, 45</td><td>29, 67</td><td>32, 104</td><td>35, 309</td><td>1</td></tr><tr><td>2014</td><td>7, 21</td><td>14, 21</td><td>16, 30</td><td>17, 33</td><td>25, 69</td><td>31, 103</td><td>32, 277</td><td>3</td></tr><tr><td>2015</td><td>15, 44</td><td>13, 18</td><td>16, 20</td><td>15, 27</td><td>15, 28</td><td>19, 76</td><td>19, 213</td><td>3</td></tr><tr><td>2016</td><td>-</td><td>13, 31</td><td>19, 43</td><td>18, 34</td><td>21, 42</td><td>21, 73</td><td>21, 223</td><td>2</td></tr><tr><td>2017</td><td>-</td><td>-</td><td>18, 40</td><td>21, 33</td><td>23, 41</td><td>21, 42</td><td>23, 156</td><td>-</td></tr><tr><td>2018</td><td>13, 44</td><td>20, 51</td><td>22, 37</td><td>28, 53</td><td>28, 55</td><td>27, 89</td><td>28, 329</td><td>5</td></tr><tr><td>2019</td><td>19, 71</td><td>26, 40</td><td>39, 90</td><td>45, 80</td><td>48, 87</td><td>45, 147</td><td>49, 515</td><td>19</td></tr><tr><td>2021</td><td>10, 40</td><td>-</td><td>25, 85</td><td>27, 69</td><td>34, 114</td><td>42, 164</td><td>42, 472</td><td>11</td></tr><tr><td>2022</td><td>22, 209</td><td>37, 209</td><td>40, 150</td><td>35, 84</td><td>59, 298</td><td>70, 337</td><td>72, 1287</td><td>9</td></tr><tr><td>2023</td><td>8, 28</td><td>13, 24</td><td>14, 25</td><td>19, 54</td><td>30, 86</td><td>35, 144</td><td>36, 361</td><td>6</td></tr></tbody></table>

## Appendix B

The AERONET inverts the aerosol radiative forcing (ΔF) and aerosol radiative forcing efficiency (ΔF<sup>eff</sup>) at BOA, confirming the following equations [72]:

$$ \Delta F_{BOA} = F_{BOA}^{\downarrow} - F_{BOA}^{\downarrow 0} \qquad (A1) $$

$$ \Delta F_{BOA}^{eff} = \Delta F_{BOA} / AOD(550 \text{ nm}) \qquad (A2) $$

where ↓ stands for downward flux in Wm<sup>-2</sup> and 0 stands for an aerosol-free condition in the atmosphere. However, the inverted radiative properties by the GRASP algorithm are estimated as follows:

$$ \Delta F_{BOA}^{Net} = (F_{BOA}^{\downarrow} - F_{BOA}^{\uparrow}) - (F_{BOA}^{\downarrow 0} - F_{BOA}^{\uparrow 0}) \qquad (A3) $$

where ΔF<sup>Net</sup><sub>BOA</sub> is the aerosol net radiative forcing [64]. Therefore, aerosol radiative forcing efficiency by the GRASP code is

$$ \Delta F_{GRASP,BOA}^{eff} = \Delta F_{BOA}^{Net} / AOD(550 \text{ nm}) \qquad (A4) $$

<footnote n="1">¹ The total does not include the repeated days.</footnote>
## References

1. Charlson, R.J.; Schwartz, S.E.; Hales, J.M.; Cess, R.D.; Coakley, J.A.; Hansen, J.E.; Hofmann, D.J. Climate Forcing by Anthropogenic Aerosols. *Science* **1992**, *255*, 423–430. [CrossRef] [PubMed]

2. Bergstrom, R.W.; Pilewskie, P.; Russell, P.B.; Redemann, J.; Bond, T.C.; Quinn, P.K.; Sierau, B. Spectral absorption properties of atmospheric aerosols. *Atmos. Chem. Phys.* **2007**, *7*, 5937–5943. [CrossRef]

3. Albrecht, B.A. Aerosols, Cloud Microphysics, and Fractional Cloudiness. *Science* **1989**, *245*, 1227–1230. [CrossRef] [PubMed]

4. Twomey, S. Pollution and the planetary albedo. *Atmos. Environ.* **1974**, *8*, 1251–1256. [CrossRef]

5. Foyo-Moreno, I.; Alados, I.; Antón, M.; Fernández-Gálvez, J.; Cazorla, A.; Alados-Arboledas, L. Estimating aerosol characteristics from solar irradiance measurements at an urban location in southeastern Spain. *J. Geophys. Res. Atmos.* **2014**, *119*, 1845–1859. [CrossRef]

6. Ångström, A. On the Atmospheric Transmission of Sun Radiation and on Dust in the Air. *Geogr. Ann.* **1929**, *11*, 156. [CrossRef]

7. IPCC. IPCC Summary for Policymakers. In *Climate Change 2021: The Physical Science Basis*; Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; IPCC: Geneva, Switzerland, 2021; ISBN 978-92-9169-158-6.

8. Forster, P.; Ramaswamy, V.; Artaxo, P.; Berntsen, T.; Betts, R.; Fahey, D.W.; Haywood, J.; Lean, J.; Lowe, D.C.; Raga, G.; et al. Changes in Atmospheric Constituents and in Radiative Forcing. In *Climate Change 2007: The Physical Science Basis*; Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK; New York, NY, USA, 2007.

9. Takayabu, Y.N.; Ueno, T.; Nakajima, T.; Matsui, I.; Tsushima, Y.; Aoki, K.; Sugimoto, N.; Uno, I. Estimate of the Cloud and Aerosol Effects on the Surface Radiative Flux Based on the Measurements and the Transfer Model Calculations Part I: Shortwave Forcing at Tateno, Japan. *J. Meteorol. Soc. Jpn.* **1999**, *77*, 1007–1021. [CrossRef]

10. Meloni, D.; Di Sarra, A.; Di Iorio, T.; Fiocco, G. Direct radiative forcing of Saharan dust in the Mediterranean from measurements at Lampedusa Island and MISR space-borne observations. *J. Geophys. Res.* **2004**, *109*, 2003JD003960. [CrossRef]

11. Charlson, R.J.; Langner, J.; Rodhe, H.; Leovy, C.B.; Warren, S.G. Perturbation of the northern hemisphere radiative balance by backscattering from anthropogenic sulfate aerosols. *Tellus A* **1991**, *43*, 152–163. [CrossRef]

12. Bush, B.C.; Valero, F.P.J. Surface aerosol radiative forcing at Gosan during the ACE-Asia campaign. *J. Geophys. Res.* **2003**, *118*, 8660. [CrossRef]

13. Norris, J.R.; Wild, M. Trends in aerosol radiative effects over Europe inferred from observed cloud cover, solar “dimming,” and solar “brightening”. *J. Geophys. Res.* **2007**, *112*, D08214. [CrossRef]

14. Valero, F.P.J.; Ackerman, T.P.; Gore, W.J.Y. The absorption of solar radiation by the Arctic atmosphere during the haze season and its effects on the radiation balance. *Geophys. Res. Lett.* **1984**, *11*, 465–468. [CrossRef]

15. Horvath, H.; Alados Arboledas, L.; Olmo, F.J.; Jovanović, O.; Gangl, M.; Kaller, W.; Sánchez, C.; Sauerzopf, H.; Seidl, S. Optical characteristics of the aerosol in Spain and Austria and its effect on radiative forcing. *J. Geophys. Res.* **2002**, *107*, 4386. [CrossRef]

16. Lyamani, H.; Olmo, F.J.; Alcántara, A.; Alados-Arboledas, L. Atmospheric aerosols during the 2003 heat wave in southeastern Spain II: Microphysical columnar properties and radiative forcing. *Atmos. Environ.* **2006**, *40*, 6465–6476. [CrossRef]

17. Díaz, A.M.; García, O.E.; Díaz, J.P.; Expósito, F.J.; Utrillas, M.P.; Martínez-Lozano, J.A.; Alados-Arboledas, L.; Olmo, F.J.; Lorente, J.; Cachorro, V.; et al. Aerosol radiative forcing efficiency in the UV region over southeastern Mediterranean: VELETA2002 campaign. *J. Geophys. Res.* **2007**, *112*, 2006JD007348. [CrossRef]

18. Cachorro, V.E.; Toledano, C.; Prats, N.; Sorribas, M.; Mogo, S.; Berjón, A.; Torres, B.; Rodrigo, R.; De La Rosa, J.; De Frutos, A.M. The strongest desert dust intrusion mixed with smoke over the Iberian Peninsula registered with Sun photometry. *J. Geophys. Res.* **2008**, *113*, 2007JD009582. [CrossRef]

19. Antón, M.; Gil, J.E.; Fernández-Gálvez, J.; Lyamani, H.; Valenzuela, A.; Foyo-Moreno, I.; Olmo, F.J.; Alados-Arboledas, L. Evaluation of the aerosol forcing efficiency in the UV erythemal range at Granada, Spain. *J. Geophys. Res.* **2011**, *116*, D20214. [CrossRef]

20. Antón, M.; Valenzuela, A.; Cazorla, A.; Gil, J.E.; Fernández-Gálvez, J.; Lyamani, H.; Foyo-Moreno, I.; Olmo, F.J.; Alados-Arboledas, L. Global and diffuse shortwave irradiance during a strong desert dust episode at Granada (Spain). *Atmos. Res.* **2012**, *118*, 232–239. [CrossRef]

21. Valenzuela, A.; Olmo, F.J.; Lyamani, H.; Antón, M.; Quirantes, A.; Alados-Arboledas, L. Aerosol radiative forcing during African desert dust events (2005–2010) over Southeastern Spain. *Atmos. Chem. Phys.* **2012**, *12*, 10331–10351. [CrossRef]

22. Granados-Muñoz, M.J.; Sicard, M.; Román, R.; Benavent-Oltra, J.A.; Barragán, R.; Brogniez, G.; Denjean, C.; Mallet, M.; Formenti, P.; Torres, B.; et al. Impact of mineral dust on shortwave and longwave radiation: Evaluation of different vertically resolved parameterizations in 1-D radiative transfer computations. *Atmos. Chem. Phys.* **2019**, *19*, 523–542. [CrossRef]

23. Dubuisson, P.; Giraud, V.; Chomette, O.; Chepfer, H.; Pelon, J. Fast radiative transfer modeling for infrared imaging radiometry. *J. Quant. Spectrosc. Radiat. Transf.* **2005**, *95*, 201–220. [CrossRef]
24. Dubuisson, P.; Buriez, J.C.; Fouquart, Y. High spectral resolution solar radiative transfer in absorbing and scattering media: Application to the satellite simulation. *J. Quant. Spectrosc. Radiat. Transf.* **1996**, *55*, 103–126. [CrossRef]

25. Lozano, I.L.; Sánchez-Hernández, G.; Guerrero-Rascado, J.L.; Alados, I.; Foyo-Moreno, I. Aerosol radiative effects in photosynthetically active radiation and total irradiance at a Mediterranean site from an 11-year database. *Atmos. Res.* **2021**, *255*, 105538. [CrossRef]

26. Sicard, M.; Bertolín, S.; Mallet, M.; Dubuisson, P.; Comerón, A. Estimation of mineral dust long-wave radiative forcing: Sensitivity study to particle properties and application to real cases in the region of Barcelona. *Atmos. Chem. Phys.* **2014**, *14*, 9213–9231. [CrossRef]

27. Sicard, M.; Barragan, R.; Dulac, F.; Alados-Arboledas, L.; Mallet, M. Aerosol optical, microphysical and radiative properties at regional background insular sites in the western Mediterranean. *Atmos. Chem. Phys.* **2016**, *16*, 12177–12203. [CrossRef]

28. Satheesh, S.K.; Ramanathan, V. Large differences in tropical aerosol forcing at the top of the atmosphere and Earth's surface. *Nature* **2000**, *405*, 60–63. [CrossRef]

29. Sicard, M.; Izquierdo, R.; Alarcón, M.; Belmonte, J.; Comerón, A.; Baldasano, J.M. Near-surface and columnar measurements with a micro pulse lidar of atmospheric pollen in Barcelona, Spain. *Atmos. Chem. Phys.* **2016**, *16*, 6805–6821. [CrossRef]

30. Sicard, M.; Rocadenbosch, F.; Reba, M.N.M.; Comerón, A.; Tomás, S.; García-Vízcamino, D.; Batet, O.; Barrios, R.; Kumar, D.; Baldasano, J.M. Seasonal variability of aerosol optical properties observed by means of a Raman lidar at an EARLINET site over Northeastern Spain. *Atmos. Chem. Phys.* **2011**, *11*, 175–190. [CrossRef]

31. Sicard, M.; Jorba, O.; Ho, J.J.; Izquierdo, R.; De Linares, C.; Alarcón, M.; Comerón, A.; Belmonte, J. Measurement report: Characterization of the vertical distribution of airborne Pinus pollen in the atmosphere with lidar-derived profiles—A modeling case study in the region of Barcelona, NE Spain. *Atmos. Chem. Phys.* **2021**, *21*, 17807–17832. [CrossRef]

32. dos Santos Oliveira, D.C.F.; Sicard, M.; Rodríguez-Gómez, A.; Comerón, A.; Muñoz-Porcar, C.; Gil-Díaz, C.; Lolli, S.; Dubovik, O.; Lopatin, A.; Herrera, M.E.; et al. Evaluation of the Accuracy of the Aerosol Optical and Microphysical Retrievals by the GRASP Algorithm from Combined Measurements of a Polarized Sun-Sky-Lunar Photometer and a Three-Wavelength Elastic Lidar. *Remote Sens.* **2023**, *15*, 5010. [CrossRef]

33. Holben, B.N.; Eck, T.F.; Slutsker, I.; Tanré, D.; Buis, J.P.; Setzer, A.; Vermote, E.; Reagan, J.A.; Kaufman, Y.J.; Nakajima, T.; et al. AERONET—A Federated Instrument Network and Data Archive for Aerosol Characterization. *Remote Sens. Environ.* **1998**, *66*, 1–16. [CrossRef]

34. Welton, E.J.; Campbell, J.R.; Berkoff, T.A.; Valencia, S.; Spinhirne, J.D.; Holben, B.; Tsay, S.-C. *The NASA Micro-Pulse Lidar Network (Mplnet): Co-Location of Lidars with Aeronet Sunphotometers and Related Earth Science Applications*; NASA: Washington, DC, USA, 2004. Available online: https://ams.confex.com/ams/Annual2005/techprogram/paper_86086.htm (accessed on 13 April 2025).

35. Pappalardo, G.; Amodeo, A.; Apituley, A.; Comeron, A.; Freudenthaler, V.; Linné, H.; Ansmann, A.; Bösenberg, J.; D'Amico, G.; Mattis, I.; et al. EARLINET: Towards an advanced sustainable European aerosol lidar network. *Atmos. Meas. Tech.* **2014**, *7*, 2389–2409. [CrossRef]

36. Giles, D.M.; Sinyuk, A.; Sorokin, M.G.; Schafer, J.S.; Smirnov, A.; Slutsker, I.; Eck, T.F.; Holben, B.N.; Lewis, J.R.; Campbell, J.R.; et al. Advancements in the Aerosol Robotic Network (AERONET) Version 3 database—Automated near-real-time quality control algorithm with improved cloud screening for Sun photometer aerosol optical depth (AOD) measurements. *Atmos. Meas. Tech.* **2019**, *12*, 169–209. [CrossRef]

37. CNR4 Net Radiometer: Instruction Manual; Kipp & Zonen: Delft, The Netherlands, 2014. Available online: https://www.kippzonen.com/Download/354/Manual-CNR-4-Net-Radiometer-English-V2104?ShowInfo=true (accessed on 13 April 2025).

38. Karim, S.A.A.; Singh, B.S.M. Global solar radiation modeling using polynomial fitting. *Appl. Math. Sci.* **2014**, *8*, 367–378. [CrossRef]

39. Bayrakçı, H.C.; Demircan, C.; Keçebaş, A. The development of empirical models for estimating global solar radiation on horizontal surface: A case study. *Renew. Sustain. Energy Rev.* **2018**, *81*, 2771–2782. [CrossRef]

40. Nethwadzi, L.C.; Winkler, H. Calculation of direct and diffuse solar irradiance components using a Slob Algorithm model in Gauteng conditions. In Proceedings of the SAIP2017, Stellenbosch, South Africa, 3–7 July 2017; pp. 361–366. Available online: https://events.saip.org.za/event/91/attachments/2266/3057/SAIP2017_361-366.pdf (accessed on 13 April 2025).

41. Kaplan, A.G.; Kaplan, Y.A. Developing of the new models in solar radiation estimation with curve fitting based on moving least-squares approximation. *Renew. Energy* **2020**, *146*, 2462–2471. [CrossRef]

42. Alados-Arboledas, L.; Olmo, F.J.; Alados, I.; Pérez, M. Parametric models to estimate photosynthetically active radiation in Spain. *Agric. For. Meteorol.* **2000**, *101*, 187–201. [CrossRef]

43. García, O.E.; Díaz, A.M.; Expósito, F.J.; Díaz, J.P.; Redondas, A.; Sasaki, T. Aerosol Radiative Forcing and Forcing Efficiency in the UVB for Regions Affected by Saharan and Asian Mineral Dust. *J. Atmos. Sci.* **2009**, *66*, 1033–1040. [CrossRef]

44. Di Sarra, A.; Pace, G.; Meloni, D.; De Silvestri, L.; Piacentino, S.; Monteleone, F. Surface shortwave radiative forcing of different aerosol types in the central Mediterranean. *Geophys. Res. Lett.* **2008**, *35*, 2007GL032395. [CrossRef]
45. Di Biagio, C.; Di Sarra, A.; Meloni, D.; Monteleone, F.; Piacentino, S.; Sferlazzo, D. Measurements of Mediterranean aerosol radiative forcing and influence of the single scattering albedo. *J. Geophys. Res.* **2009**, *114*, 2008JD011037. [CrossRef]

46. Di Biagio, C.; Di Sarra, A.; Meloni, D. Large atmospheric shortwave radiative forcing by Mediterranean aerosols derived from simultaneous ground-based and spaceborne observations and dependence on the aerosol type and single scattering albedo. *J. Geophys. Res.* **2010**, *115*, 2009JD012697. [CrossRef]

47. Pace, G.; di Sarra, A.; Meloni, D.; Piacentino, S.; Chamard, P. Aerosol optical properties at Lampedusa (Central Mediterranean). 1. Influence of transport and identification of different aerosol types. *Atmos. Chem. Phys.* **2006**, *6*, 697–713. [CrossRef]

48. Toledano, C.; Cachorro, V.E.; Berjon, A.; De Frutos, A.M.; Sorribas, M.; De La Morena, B.A.; Goloub, P. Aerosol optical depth and Ångström exponent climatology at El Arenosillo AERONET site (Huelva, Spain). *Q. J. R. Meteorol. Soc.* **2007**, *133*, 795–807. [CrossRef]

49. Dubovik, O.; Lapyonok, T.; Litvinov, P.; Herman, M.; Fuertes, D.; Ducos, F.; Torres, B.; Derimian, Y.; Huang, X.; Lopatin, A.; et al. GRASP: A versatile algorithm for characterizing the atmosphere. *SPIE Newsroom* **2014**, 25. [CrossRef]

50. Dubovik, O.; Herman, M.; Holdak, A.; Lapyonok, T.; Tanré, D.; Deuzé, J.L.; Ducos, F.; Sinyuk, A.; Lopatin, A. Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations. *Atmos. Meas. Tech.* **2011**, *4*, 975–1018. [CrossRef]

51. Dubovik, O.; Smirnov, A.; Holben, B.N.; King, M.D.; Kaufman, Y.J.; Eck, T.F.; Slutsker, I. Accuracy assessments of aerosol optical properties retrieved from Aerosol Robotic Network (AERONET) Sun and sky radiance measurements. *J. Geophys. Res.* **2000**, *105*, 9791–9806. [CrossRef]

52. Dubovik, O.; King, M.D. A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements. *J. Geophys. Res.* **2000**, *105*, 20673–20696. [CrossRef]

53. Dubovik, O.; Holben, B.N.; Lapyonok, T.; Sinyuk, A.; Mishchenko, M.I.; Yang, P.; Slutsker, I. Non-spherical aerosol retrieval method employing light scattering by spheroids: Non-Spherical Aerosol Retrieval Method. *Geophys. Res. Lett.* **2002**, *29*, 1415. [CrossRef]

54. Dubovik, O.; Sinyuk, A.; Lapyonok, T.; Holben, B.N.; Mishchenko, M.; Yang, P.; Eck, T.F.; Volten, H.; Muñoz, O.; Vehelmann, B.; et al. Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust. *J. Geophys. Res.* **2006**, *111*, D11208. [CrossRef]

55. Lopatin, A.; Dubovik, O.; Chaikovsky, A.; Goloub, P.; Lapyonok, T.; Tanré, D.; Litvinov, P. Enhancement of aerosol characterization using synergy of lidar and sun-photometer coincident observations: The GARRLiC algorithm. *Atmos. Meas. Tech.* **2013**, *6*, 2065–2088. [CrossRef]

56. Ross, J. Tasks for Vegetation Science. In *The Radiation Regime and Architecture of Plant Stands*; Springer: Dordrecht, The Netherlands, 1981; ISBN 978-90-6193-607-7.

57. Li, X.; Strahler, A.H. Geometric-optical bidirectional reflectance modeling of the discrete crown vegetation canopy: Effect of crown shape and mutual shadowing. *IEEE Trans. Geosci. Remote Sens.* **1992**, *30*, 276–292. [CrossRef]

58. Wanner, W.; Li, X.; Strahler, A.H. On the derivation of kernels for kernel-driven models of bidirectional reflectance. *J. Geophys. Res.* **1995**, *100*, 21077–21089. [CrossRef]

59. Nadal, F.; Breon, F.-M. Parameterization of surface polarized reflectance derived from POLDER spaceborne measurements. *IEEE Trans. Geosci. Remote Sens.* **1999**, *37*, 1709–1718. [CrossRef]

60. Maignan, F.; Bréon, F.-M.; Fédèle, E.; Bouvier, M. Polarized reflectances of natural surfaces: Spaceborne measurements and analytical modeling. *Remote Sens. Environ.* **2009**, *113*, 2642–2650. [CrossRef]

61. Maignan, F.; Bréon, F.-M.; Lacaze, R. Bidirectional reflectance of Earth targets: Evaluation of analytical models using a large set of spaceborne measurements with emphasis on the Hot Spot. *Remote Sens. Environ.* **2004**, *90*, 210–220. [CrossRef]

62. Deschamps, P.-Y.; Breon, F.-M.; Leroy, M.; Podaire, A.; Bricaud, A.; Buriez, J.-C.; Seze, G. The POLDER mission: Instrument characteristics and scientific objectives. *IEEE Trans. Geosci. Remote Sens.* **1994**, *32*, 598–615. [CrossRef]

63. Litvinov, P.; Hasekamp, O.; Cairns, B. Models for surface reflection of radiance and polarized radiance: Comparison with airborne multi-angle photopolarimetric measurements and implications for modeling top-of-atmosphere measurements. *Remote Sens. Environ.* **2011**, *115*, 781–792. [CrossRef]

64. Derimian, Y.; Dubovik, O.; Huang, X.; Lapyonok, T.; Litvinov, P.; Kostinski, A.B.; Dubuisson, P.; Ducos, F. Comprehensive tool for calculation of radiative fluxes: Illustration of shortwave aerosol radiative effect sensitivities to the details in aerosol and underlying surface characteristics. *Atmos. Chem. Phys.* **2016**, *16*, 5763–5780. [CrossRef]

65. Dubovik, O.; Fuertes, D.; Litvinov, P.; Lopatin, A.; Lapyonok, T.; Doubovik, I.; Xu, F.; Ducos, F.; Chen, C.; Torres, B.; et al. A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Constraints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications. *Front. Remote Sens.* **2021**, *2*, 706851. [CrossRef]

66. D'Amico, G.; Amodeo, A.; Mattis, I.; Freudenthaler, V.; Pappalardo, G. EARLINET Single Calculus Chain—Technical—Part 1: Pre-processing of raw lidar data. *Atmos. Meas. Tech.* **2016**, *9*, 491–507. [CrossRef]
67. D'Amico, G.; Amodeo, A.; Baars, H.; Binietoglou, I.; Freudenthaler, V.; Mattis, I.; Wandinger, U.; Pappalardo, G. EARLINET Single Calculus Chain—Overview on methodology and strategy. *Atmos. Meas. Tech.* **2015**, *8*, 4891–4916. [CrossRef]

68. Mattis, I.; D'Amico, G.; Baars, H.; Amodeo, A.; Madonna, F.; Iarlori, M. EARLINET Single Calculus Chain—technical—Part 2: Calculation of optical products. *Atmos. Meas. Tech.* **2016**, *9*, 3009–3029. [CrossRef]

69. Mishchenko, M.I.; Travis, L.D. Satellite retrieval of aerosol properties over the ocean using polarization as well as intensity of reflected sunlight. *J. Geophys. Res.* **1997**, *102*, 16989–17013. [CrossRef]

70. Li, Z.; Goloub, P.; Devaux, C.; Gu, X.; Deuzé, J.-L.; Qiao, Y.; Zhao, F. Retrieval of aerosol optical and physical properties from ground-based spectral, multi-angular, and polarized sun-photometer measurements. *Remote Sens. Environ.* **2006**, *101*, 519–533. [CrossRef]

71. Fedarenka, A.; Dubovik, O.; Goloub, P.; Li, Z.; Lapyonok, T.; Litvinov, P.; Barel, L.; Gonzalez, L.; Podvin, T.; Crozel, D. Utilization of AERONET polarimetric measurements for improving retrieval of aerosol microphysics: GSFC, Beijing and Dakar data analysis. *J. Quant. Spectrosc. Radiat. Transf.* **2016**, *179*, 72–97. [CrossRef]

72. AERONET Document. *AERONET Inversion Products (Version 3)*; NASA: Washington, DC, USA.

73. García, O.E.; Díaz, A.M.; Expósito, F.J.; Díaz, J.P.; Dubovik, O.; Dubuisson, P.; Roger, J.-C.; Eck, T.F.; Sinyuk, A.; Derimian, Y.; et al. Validation of AERONET estimates of atmospheric solar fluxes and aerosol radiative forcing by ground-based broadband measurements. *J. Geophys. Res.* **2008**, *113*, 2008JD010211. [CrossRef]

74. Lolli, S.; Sicard, M.; Amato, F.; Comeron, A.; Gil-Diaz, C.; Landi, T.C.; Munoz-Porcar, C.; Oliveira, D.C.F.S.; Dios Otin, F.; Rocadenbosch, F.; et al. Climatological assessment of the vertically resolved optical and microphysical aerosol properties by lidar measurements, sun photometer, and in situ observations over 17 years at Universitat Politècnica de Catalunya (UPC) Barcelona. *Atmos. Chem. Phys.* **2023**, *23*, 12887–12906. [CrossRef]

75. Burkey, J. Seasonal Kendall Test with Slope for Serial Dependent Data. 2024. Available online: https://www.mathworks.com/matlabcentral/fileexchange/22389-seasonal-kendall-test-with-slope-for-serial-dependent-data (accessed on 25 June 2024).

76. Burkey, J. Mann-Kendall Tau-b with Sen's Method (Enhanced). 2024. Available online: https://www.mathworks.com/matlabcentral/fileexchange/11190-mann-kendall-tau-b-with-sen-s-method-enhanced (accessed on 25 June 2024).

77. Markowicz, K.M.; Flatau, P.J.; Remiszewska, J.; Witek, M.; Reid, E.A.; Reid, J.S.; Bucholtz, A.; Holben, B. Observations and Modeling of the Surface Aerosol Radiative Forcing during UAE2. *J. Atmos. Sci.* **2008**, *65*, 2877–2891. [CrossRef]

78. Guan, H.; Schmid, B.; Bucholtz, A.; Bergstrom, R. Sensitivity of shortwave radiative flux density, forcing, and heating rate to the aerosol vertical profile. *J. Geophys. Res.* **2010**, *115*, 2009JD012907. [CrossRef]

79. Dubovik, O.; Holben, B.; Eck, T.F.; Smirnov, A.; Kaufman, Y.J.; King, M.D.; Tanré, D.; Slutsker, I. Variability of Absorption and Optical Properties of Key Aerosol Types Observed in Worldwide Locations. *J. Atmos. Sci.* **2002**, *59*, 590–608. [CrossRef]

80. Gobbi, G.P.; Kaufman, Y.J.; Koren, I.; Eck, T.F. Classification of aerosol properties derived from AERONET direct sun data. *Atmos. Chem. Phys.* **2007**, *7*, 453–458. [CrossRef]

81. Kaufman, Y.J. Aerosol optical thickness and atmospheric path radiance. *J. Geophys. Res.* **1993**, *98*, 2677–2692. [CrossRef]

82. Basart, S.; Perez, C.; Cuevas, E.; Baldasano, J.M.; Gobbi, G.P. Aerosol characterization in Northern Africa, Northeastern Atlantic, Mediterranean Basin and Middle East from direct-sun AERONET observations. *Atmos. Chem. Phys.* **2009**, *9*, 8265–8282. [CrossRef]

83. Salvador, P.; Pey, J.; Pérez, N.; Querol, X.; Artíñano, B. Increasing atmospheric dust transport towards the western Mediterranean over 1948–2020. *npj Clim Atmos Sci* **2022**, *5*, 34. [CrossRef]

84. Pandolfi, M.; Alastuey, A.; Pérez, N.; Reche, C.; Castro, I.; Shatalov, V.; Querol, X. Trends analysis of PM source contributions and chemical tracers in NE Spain during 2004–2014: A multi-exponential approach. *Atmos. Chem. Phys.* **2016**, *16*, 11787–11805. [CrossRef]

85. Veld, M.I.T.; Alastuey, A.; Pandolfi, M.; Amato, F.; Pérez, N.; Reche, C.; Via, M.; Minguillón, M.C.; Escudero, M.; Querol, X. Compositional changes of PM<sub>2.5</sub> in NE Spain during 2009–2018: A trend analysis of the chemical composition and source apportionment. *Sci. Total Environ.* **2021**, *795*, 148728. [CrossRef]

86. Querol, X. PM10 and PM<sub>2.5</sub> source apportionment in the Barcelona Metropolitan area, Catalonia, Spain. *Atmos. Environ.* **2001**, *35*, 6407–6419. [CrossRef]

87. Via, M.; Minguillón, M.C.; Reche, C.; Querol, X.; Alastuey, A. Increase in secondary organic aerosol in an urban environment. *Atmos. Chem. Phys.* **2021**, *21*, 8323–8339. [CrossRef]

88. Aas, W.; Fagerli, H.; Alastuey, A.; Cavalli, F.; Degorska, A.; Feigenspan, S.; Brenna, H.; Gliß, J.; Heinesen, D.; Hueglin, C.; et al. Trends in Air Pollution in Europe, 2000–2019. *Aerosol Air Qual. Res.* **2024**, *24*, 230237. [CrossRef]

89. Sabburg, J.; Wong, J. The effect of clouds on enhancing UVB irradiance at the Earth's surface: A one year study. *Geophys. Res. Lett.* **2000**, *27*, 3337–3340. [CrossRef]

90. Badosa, J.; Calbó, J.; Mckenzie, R.; Liley, B.; González, J.; Forgan, B.; Long, C.N. Two Methods for Retrieving UV Index for All Cloud Conditions from Sky Imager Products or Total SW Radiation Measurements. *Photochem Photobiol.* **2014**, *90*, 941–951. [CrossRef] [PubMed]
91. Lamy, K.; Portafaix, T.; Brogniez, C.; Lakkala, K.; Pitkänen, M.R.A.; Arola, A.; Forestier, J.-B.; Amelie, V.; Toihir, M.A.; Rakotoniaina, S. UV-Indien network: Ground-based measurements dedicated to the monitoring of UV radiation over the western Indian Ocean. *Earth Syst. Sci. Data* **2021**, *13*, 4275–4301. [CrossRef]

92. Lamy, K.; Portafaix, T.; Brogniez, C.; Godin-Beekmann, S.; Bencherif, H.; Morel, B.; Pazmino, A.; Metzger, J.M.; Auriol, F.; Deroo, C.; et al. Ultraviolet radiation modelling from ground-based and satellite measurements on Reunion Island, southern tropics. *Atmos. Chem. Phys.* **2018**, *18*, 227–246. [CrossRef]

93. Di Sarra, A.; Fuà, D.; Meloni, D. Estimate of surface direct radiative forcing of desert dust from atmospheric modulation of the aerosol optical depth. *Atmos. Chem. Phys.* **2013**, *13*, 5647–5654. [CrossRef]

94. Meloni, D.; Di Sarra, A.; Di Iorio, T.; Fiocco, G. Influence of the vertical profile of Saharan dust on the visible direct radiative forcing. *J. Quant. Spectrosc. Radiat. Transf.* **2005**, *93*, 397–413. [CrossRef]

95. Boucher, O.; Schwartz, S.E.; Ackerman, T.P.; Anderson, T.L.; Bergstrom, B.; Bonnel, B.; Chýlek, P.; Dahlback, A.; Fouquart, Y.; Fu, Q.; et al. Intercomparison of models representing direct shortwave radiative forcing by sulfate aerosols. *J. Geophys. Res.* **1998**, *103*, 16979–16998. [CrossRef]

96. Obregón, M.A.; Costa, M.J.; Silva, A.M.; Serrano, A. Spatial and Temporal Variation of Aerosol and Water Vapour Effects on Solar Radiation in the Mediterranean Basin during the Last Two Decades. *Remote Sens.* **2020**, *12*, 1316. [CrossRef]

97. Logothetis, S.-A.; Salamalikis, V.; Kazantzidis, A. The impact of different aerosol properties and types on direct aerosol radiative forcing and efficiency using AERONET version 3. *Atmos. Res.* **2021**, *250*, 105343. [CrossRef]

98. Chen, A.; Zhao, C.; Shen, L.; Fan, T. Influence of aerosol properties and surface albedo on radiative forcing efficiency of key aerosol types using global AERONET data. *Atmos. Res.* **2023**, *282*, 106519. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.