A Machine Learning Approach to Retrieving Aerosol Optical Depth Using Solar Radiation Measurements

Abstract

Aerosol optical depth (AOD) constitutes a key parameter of aerosols, providing vital information for quantifying the aerosol burden and air quality at global and regional levels. This study demonstrates a machine learning strategy for retrieving AOD under cloud-free conditions based on the synergy of machine learning algorithms (MLAs) and ground-based solar irradiance data. The performance of the proposed methodology was investigated by applying different components of solar irradiance. In particular, the use of direct instead of global irradiance as a model feature led to better performance. The MLA-based AODs were compared to reference AERONET retrievals, which encompassed RMSE values between 0.01 and 0.15, regardless of the underlying climate and aerosol environments. Among the MLAs, artificial neural networks outperformed the other algorithms in terms of RMSE at 54% of the measurement sites. The overall performance of MLA-based AODs against AERONET revealed a high coefficient of determination (R² = 0.97), MAE of 0.01, and RMSE of 0.02. Compared to satellite (MODIS) and reanalysis (MERRA-2 and CAMSRA) data, the MLA-AOD retrievals revealed the highest accuracy at all stations. The ML-AOD retrievals have the potential to expand and complement the AOD information in non-existing timeframes when solar irradiances are available.

1. Introduction

Aerosol particles, suspended in the atmosphere and originating from either natural (e.g., dust and sea salt) or anthropogenic sources (e.g., nitrate and sulfate), exert a direct effect on the Earth’s radiative balance through aerosol–radiation interaction and an indirect effect as their role as condensation nuclei in cloud formation [1]. This interaction is contingent upon the aerosol optical properties [2]. While these optical properties are crucial for evaluating the radiative effect of aerosols [3,4,5,6], they also hold significance for precise solar resource estimations on both regional and global scales [7].

Reference instruments for measuring solar radiation face limitations in terms of fine temporal resolution and worldwide spatial coverage. Consequently, alternative techniques for solar resource assessment are under scrutiny, with a primary focus on radiative transfer modeling. An accurate assessment of solar resources necessitates a comprehensive understanding of aerosol properties. Integrating these properties into radiative transfer models enables the simulation of long-term solar irradiances under cloud-free conditions [8,9,10,11,12,13,14,15,16,17].

There are two primary methods for monitoring atmospheric aerosols. The first method involves the use of ground-based instruments, such as sun photometers. AERONET (Aerosol Robotic NETwork) [18] is a widely recognized ground-based monitoring network of sun photometers that has been providing highly accurate physical and optical characteristics of aerosols at over 600 stations worldwide since 1993. The second method involves utilizing sensors aboard satellites, such as the MODerate resolution Imaging Spectroradiometer (MODIS), Multi-angle Imaging Spectro-Radiometer (MISR), and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO). In general, ground-based and satellite-based aerosol retrievals complement each other. Ground-based retrievals excel in accuracy and temporal resolution but face limitations in spatial coverage, particularly in desert areas, which are significant natural sources of suspended aerosol particles in the atmosphere [19,20,21]. On the other hand, polar-orbiting satellites provide high spatial resolution of aerosol optical depth (AOD) but are limited to offering more than two retrievals per day across a specific region. On the other hand, geostationary satellites, by utilizing their higher temporal resolution, can present almost real-time images and capture the dynamic AOD variations, which is very useful for monitoring events such as forest fires [22]. However, there are two main drawbacks when using geostationary satellites for earth observations. The first is their limited storage and high orbit altitude, which make it challenging to record images with high spatial resolution. The second drawback is that higher uncertainties are observed in AOD retrieval algorithms compared to polar orbit satellites [23]. Considering the spatial and temporal constraints of both satellite and ground-based aerosol retrievals, other techniques have recently emerged, with a predominant focus on AOD.

Various alternative techniques for AOD retrieval reported in the literature fall into four primary categories: (1) backward solving radiative transfer (RT) or clear-sky models using solar radiation measurements [24,25,26,27], (2) methodologies based on sunshine duration (SD) measurements [28,29,30], (3) image processing techniques using sky radiances from all-sky imagers [31,32,33,34], and (4) machine learning (ML) and deep learning algorithms employing various independent parameters as input features [35,36,37,38,39].

Salmon et al. [25] retrieved AOD across 16 regions worldwide characterized by different climates and aerosol properties, solving backward two clear-sky radiation models (SMARTS and SOLIS). The model utilized inputs such as broadband direct normal irradiance (DNI), meteorological data, and aerosol information such as single scattering albedo and the Ångström exponent, resulting in an average dispersion error of approximately 0.04 for the retrieved AOD compared to AERONET. Wandji Nyamsi et al. [30] conducted a benchmarking analysis, evaluating three distinct methodologies (linear-based, physical-based, and hybrid) for AOD retrieval based on SD measurements over Europe. The hybrid approach outperformed the other methods, yielding a root mean square error (RMSE) ranging from 0.04 to 0.13 compared to AERONET. Román et al. [32] retrieved AOD at three different wavelengths (465, 536, and 605 nm) using sky radiances from an all-sky imager and the GRASP (Generalized Retrieval of Atmosphere and Surface Properties) [40] code. The predicted AOD agreed well with AERONET retrievals, exhibiting a coefficient of determination of approximately 0.87. Huttunen et al. [36] employed five ML algorithms (MLAs: non-linear regression, artificial neural network, random forest, Gaussian process regression, and support vector machine) to retrieve AOD. All models included global horizontal irradiance (GHI), total column water vapor, and solar geometry as independent parameters. The predicted AOD closely matched observations (absolute bias = 0.04–0.05), with the Gaussian process regression and artificial neural network to outperform other ML algorithms, revealing coefficients of determination of 0.93 and 0.92, respectively.

The primary objective of this study is to implement an ML approach for retrieving AOD, denoted as “MLA-AOD”. The innovative aspects of this proposed approach are succinctly outlined below:

• AOD is extracted utilizing the direct beam of solar radiation or direct normal irradiance (DNI) derived from ground-based reference measurements.
• The impact of different solar irradiance components (GHI vs. DNI) on the retrieved AOD is explored, underscoring the significance of selecting the most appropriate solar irradiance component for AOD retrieval.
• The validation process includes comparisons against AERONET measurement locations with diverse climatic and aerosol characteristics.
• MLAs with different prediction mechanisms are benchmarked to identify the “optimal” approach for retrieving AOD.
• The MLA-AOD retrievals are compared with AOD from reanalysis products (MERRA-2 and CAMSRA) through a comprehensive comparative analysis, showcasing the advantages of the proposed methodology at a regional level.

The proposed MLA-AOD methodology presents an aerosol retrieval technique for cloud-free conditions that can be used to reproduce the atmospheric aerosol burden at the high temporal resolution of radiometric instruments. The MLA-AOD approach is evaluated using either pyrheliometer (DNI) or pyranometer (GHI) measurements, presenting different retrieval performances. The selection of the applied solar irradiance component relies on the instrument availability. Even though AERONET sun photometers are established and measured every 5–15 min, it is very common for them to present gaps of hours, days, or weeks or to be inactive for years. This study proposes that the MLA-AOD methodology could be used (a) to fill the AOD gaps in the AERONET sites and (b) to expand the already existing temporal capabilities, regardless of the underlying climate and aerosol conditions. Upgrading the current AOD capabilities is important for the evaluation of satellite products, aerosol data assimilation, air pollution models, AOD trend analyses, climate models, and direct aerosol radiative effect estimation.

The structure of this work is as follows: Section 2 and Section 3 present details on the datasets and the proposed methodology. Section 4 covers the results and Section 5 addresses the summary and conclusion of this study.

2. Datasets

2.1. Ground-Based Data

2.1.1. AERONET

AERONET (https://aeronet.gsfc.nasa.gov; accessed on 12 January 2024) comprises a network of sunphotometers, which are ground-based, passive remote sensing instruments [18]. It provides high-quality data on aerosol optical and physical properties, encompassing over 600 stations. AERONET’s Direct Sun Algorithm (DSA) derives AOD at seven wavelengths spanning from 340 nm to 1020 nm. DSA relies on spectral solar irradiance measurements directly on the solar disk, captured from the CIMEL Electronique CE318 multiwavelength sunphotometer. Besides AOD, DSAA also provides information about total column water vapor (TCWV), based on measurements at 940 nm. AERONET, except for DSA, applies an inverse algorithm using sky radiances in the almucantar geometry to provide other aerosol properties like the single scattering albedo, aerosol volume distribution, phase function, refractive index, etc. [41].

In this study, AOD at 550 nm (AOD_550nm) is utilized to maintain consistency with MODIS, CAMSRA, and MERRA-2 reanalysis. The AERONET AOD is interpolated at 550 nm using the Ångström law. For subsequent analyses, cloud-screened data rounded to the nearest minute are employed, meeting quality criteria established by pre-field and post-field calibrations (Level 2.0 data from Version 3; L2V3) [42]. The spectral AOD uncertainties in the visible spectrum for AERONET Version 3 are reported to be 0.01 [42].

2.1.2. BSRN

The Baseline Surface Radiation Network (BSRN, https://bsrn.awi.de; accessed on 12 January 2024), is a ground-based solar radiation monitoring network that has been measuring solar radiation at various locations worldwide since 1992 [43,44]. BSRN provides comprehensive data on the various components of incoming solar radiation, namely global horizontal irradiance (GHI), direct normal irradiance (DNI), and direct horizontal irradiance (DHI). GHI refers to the total amount of shortwave radiation reaching the ground surface. DNI is the amount of solar radiation received per unit area by a surface that is always normal to the solar rays. DHI is the total amount of solar radiation that reaches the ground and does not arrive on a direct path from the sun. The network delivers quality-assured measurements at a 1 min temporal resolution. Each BSRN radiometric station employs a first-class pyranometer and a pyrheliometer to measure GHI and DNI, respectively. The raw observations from BSRN undergo quality control [45] using the procedures outlined in the SolarData R package [46]. For subsequent analysis, only the 1 min quality-assured DNI and GHI measurements from BSRN are utilized.

2.1.3. Stations’ Characteristics

Figure 1 displays the spatial distribution of the collocated AERONET-BSRN stations (N = 26), while Table S1 provides information about the station names and data availability. The majority of the stations are situated in Europe and the US (Figure 1). The five distinct colors in Figure 1 represent the Köppen–Geiger (KG) [47] climate clusters (CCs). The KG climate classification system categorizes the world’s climate into five primary CCs based on air temperature and precipitation amount, namely (A) Equatorial, (B) Arid, (C) Warm Temperate, (D) Snow, and (E) Polar. The corresponding climate zone for each station is determined using high-resolution data (100 arc sec) via the ‘kgc’ R package [48]. Each CC is represented by at least one station, with arid and temperate climates being the most prevalent in terms of station occurrence.

It is significant to evaluate the proposed methodology in diverse environments to assess its resilience to potential changes in climate and aerosol types. Previous studies showed that the accuracy of model-based and reanalysis AOD against observations varied climatically [49,50] and regionally [50,51]. For instance, the training of an ML algorithm to retrieve AOD becomes more challenging in regions potentially affected by dust intrusions or wildfires.

2.2. Reanalysis Data

2.2.1. CAMSRA

The Copernicus Atmosphere Monitoring Service Reanalysis (CAMSRA) represents the latest global dataset for atmospheric composition, created by the European Centre for Medium-Range Weather Forecasts (ECMWF). This dataset consists of 3-D time-consistent fields of atmospheric composition, incorporating aerosols and various chemical species [52]. CAMSRA was generated through 4DVar data assimilation within CY42R1 of ECMWF’s Integrated Forecast System (IFS), in which meteorological modeling is incorporated. The 4DVar analysis applies IFS dynamics and physics to produce a sequence of atmospheric states that closely align with the available observations.

Concerning AOD_550nm, CAMSRA assimilates satellite-derived AOD_550nm retrievals obtained from MODIS Collection 6 aboard Aqua and Terra, as well as the Advanced Along-Track Scanning Radiometer (AATSR) aboard Envisat [53]. For TCWV, CAMSRA assimilates satellite-based infrared and microwave radiances from various sensors onboard satellites such as AMSR-2 (GCOM-W1), AMSRE (AQUA), GMI (GPM), and TMI (TRMM), among others. Since 2003, CAMSRA data have been provided at approximately 80 km resolution and on a 3 h basis. In this study, AOD_550nm and TCWV retrievals are acquired through CAMSRA. The validation of CAMSRA AOD_550nm against AERONET observations revealed a global Root Mean Square Error (RMSE) of 0.144 [50]. The validation of CAMSRA TCWV against AERONET retrievals revealed satisfactory performance with a coefficient of determination (R²) of 0.93 and an RMSE of 0.3 cm (Figure S1).

2.2.2. MERRA-2

The Modern-Era Retrospective Analysis for Research and Applications-2 (MERRA-2) is the latest global reanalysis dataset developed by the NASA Global Modeling and Assimilation Office (GMAO). MERRA-2 assimilates meteorological and aerosol measurements obtained from a variety of ground-based and satellite-based products [54,55,56]. This analysis focuses on AOD_550nm for intercomparison purposes. The MERRA-2 AOD_550nm reanalysis product incorporates AOD products spanning different time periods. More specifically, before 2000, AOD_550nm relies exclusively on AVHRR retrievals. Subsequently, MERRA-2 assimilates data from MODIS (onboard both Terra and Aqua, covering the period from 2000 and 2002 to the present, respectively), MISR (from 2000 to 2014), and AERONET (from 1999 to 2014). MERRA-2 data (product name: tavg1_2d_aer_Nx) are provided at a spatial resolution of 0.5° × 0.625° on an hourly basis.

The global validation of MERRA-2 AOD_550nm against AERONET data yielded an RMSE of 0.126 [50], outperforming CAMSRA. The enhanced performance of MERRA-2 AOD_550nm is attributed to its utilization of AERONET retrievals as part of its data assimilation methodology.

2.3. Satellite Data

MODIS

MODerate resolution Imaging Spectroradiometer (MODIS) is a passive sensor that has been hosted on NASA’s Terra and Aqua polar satellites since 2000 and 2002, respectively. MODIS monitors several pivotal atmospheric constituents, like aerosols [57], on a global scale. MODIS takes advantage of its wide swath (~2330 km) to provide approximate global coverage by making 14–15 orbits daily. The Aqua/Terra crosses the equator in a northbound/southbound direction at 13:30/10:30 UTC. Across dark (in the visible spectrum) surfaces over land [58,59,60] or/and ocean [61,62,63] areas, the Deep Target algorithm (DT) is applied, while over bright-reflecting land surfaces [64,65], including desert, semi-arid, and urban regions, the Deep Blue algorithm (DB) is adopted.

In this study, the AOD_550_Dark_Target_Deep_Blue_Combined product from Collection 6.1 from Aqua and Terra is applied, which is a merged dataset of DB and DT algorithms. In particular, the AOD values are merged when specific quality assurance (QA) conditions are met (QA ≥ 1 for DT ocean, QA = 3 for DT land, and QA ≥ 2 for DB land).

3. Methodology

The proposed methodology (MLA-AOD) for retrieving AOD_550nm at 1 min resolution is illustrated in Figure 2. Various data preprocessing steps are included to increase the performance of the MLAs.

The applied methodology involves two primary steps: data acquisition and preprocessing, which are extensively detailed in Section 2 and Section 3.1. Subsequently, the preprocessed dataset is divided into two subsets, namely the training and testing datasets. The splitting strategy follows a 70/30 approach without using a random-sampling approach. The first 70% of the time series represents the training set and the remainder represents the test set. Table S1 includes information about the training and testing periods for each station.

The independent (input) features undergo normalization between 0 and 1 using the Min–Max normalization method. Five different ML algorithms (see Section 3.3 and Appendix A) are evaluated to retrieve AOD. The MLA structure is determined through several internal parameters called hyperparameters. The optimal configuration for these hyperparameters for each MLA is determined through a randomized search procedure during training. The full cartesian hyperparameter space is given in Table S2. The randomized search involves a 10-fold cross-validation process, with MSE serving as the fitness function. For each random hyperparameter set, the data are randomly divided into 10 roughly equal subsamples. Then, at each iteration, one of those 10 subsamples is retained as the “validation” dataset for MLA evaluation, and the 10-1 subsamples train the ML model. The MSE scores retrieved through the 10 iterations are averaged to retrieve a total score which is characteristic for each hyperparameter set. This process is repeated 100 times. The “optimal” hyperparameter combination for each MLA is determined based on the minimum of total scores. The MLA-AOD methodology (Figure 2) is applied for each collocated AERONET-BSRN station (see Section 2.1.3).

3.1. Data Preprocessing

The raw data obtained from ground-based measurements are cloud-screened to isolate clear sky conditions. Although L2V3 AERONET retrievals are described as cloud-screened, the temporal merging of AERONET and BSRN data may introduce cloud-contaminated cases, attributed to the smaller aperture of sun photometers compared to pyrheliometers [25]. To ensure clear-sky conditions, two distinct filters were applied—one for DNI and one for GHI. Concerning DNI, within one hour around each BSRN measurement, the standard deviation is calculated. If it is above 20% of the corresponding BSRN DNI measurement, the data point is detected as cloudy and removed [25]. In addition, even though the sun’s disk could be cloudless, the rest of the sky can contain sparse clouds that affect the GHI measurements. Therefore, to avoid possible cloud contamination, the Reno and Hansen clear-sky detection methodology [66] is applied to the GHI measurements. This method is preferred here over other published methodologies because it relies on periods of time. In particular, the Reno and Hansen clear-sky detection methodology [66] is applied on time windows (here, 10 min) of 1 min temporal resolution. If at least one minute is assigned as cloudy, the whole window is flagged as cloudy and eliminated from the analysis. By applying the filters for GHI and DNI described above, the cloudy solar irradiance measurements are unlikely to be found in the used data. Subsequently, the cloudless solar irradiance measurements are synchronized with AERONET 1 min retrievals.

Moreover, the reanalysis data (TCWV and AOD_550nm) are extracted at the ground-based stations through nearest-neighbor interpolation. The 3-hourly TCWV from CAMSRA is linearly interpolated at 1 min temporal resolution, aligning with the AERONET and BSRN datasets.

For validation purposes, both MLA-AOD based and AERONET AODs are temporally aggregated to match the resolution of the reanalysis and satellite products. Given that MERRA-2 hourly data are available at half-hour intervals (e.g., 00:30, 01:30, etc.), the ground observations are averaged over ±30 min intervals, centered on the MERRA-2 timestamps. On the other hand, CAMSRA provides data with a 3 h resolution (e.g., 00:00, 03:00, etc.). In this case, the ground-based data are represented as 90 min averages, centered on the CAMSRA 3 h intervals. Regarding the satellite data, the MODIS AODs are spatially averaged if at least one measurement resides within a circle area of 25 km around the AERONET station. Then, the spatially averaged AODs are temporally averaged when at least one of the AERONET measurements falls within a time window of 2 h around the satellite overpass time.

Another significant point in validation studies involves addressing elevation differences between ground-based and reanalysis datasets [50]. To mitigate these discrepancies, the reanalysis AOD is adjusted to the station elevation using the scale-height correction of Gueymard and Yang [50].

In addition, any unrealistic AERONET AOD retrievals (AOD_550nm < 0) are eliminated for the subsequent analysis. Finally, the proposed methodology is applied using data for a 15-year period (during 2004–2019). The data period is related to the station’s data availability. Finally, for validation against reanalysis, only stations with at least a minimum of 100 (3 h) CAMSRA and 300 (1 h) MERRA-2 records are retained.

3.2. Input Parameters

The direct and global clearness indexes are calculated using Equations (1) and (2), respectively,

$${\mathrm{k}\mathrm{t}}_{\mathrm{d}\mathrm{i}\mathrm{r}}=\frac{\mathrm{D}\mathrm{N}\mathrm{I}}{{\mathrm{E}}_0}$$

(1)

$${\mathrm{k}\mathrm{t}}_{\mathrm{g}\mathrm{l}\mathrm{o}\mathrm{b}}=\frac{\mathrm{G}\mathrm{H}\mathrm{I}}{{\mathrm{E}}_0}$$

(2)

where E₀ is the received irradiance on a horizontal plane at the top of the atmosphere. The optical air mass is calculated based on Kasten and Young’s [67] formula, using Equation (3),

$$\mathrm{m}=(\mathrm{P}/{\mathrm{P}}_0)/[\cos \left(\mathrm{S}\mathrm{Z}\mathrm{A}\right)+0.50572·{\left(\mathrm{S}\mathrm{Z}\mathrm{A}+6.07995\right)}^{-1.6364}]$$

(3)

where P and P₀ are the pressure level at the station’s elevation and sea level, respectively, while SZA corresponds to the solar zenith angle. The TCWV in cm is derived through the CAMSRA (see Section 2.2.1).

The retrieval of AOD through remote sensing always requires information about the solar geometry and amount of the non-aerosol atmospheric constituents. For example, the AERONET retrieval methodology relies on Beer–Lambert–Bouguer law, considering the impact of solar geometry and TCWV [18,42]. In addition, solar geometry is considered when the actual and estimated reflectances are compared through AOD retrieval from satellite measurements [68]. In the proposed methodology, two different solar irradiance components were applied, presenting different levels of alteration with variations in aerosol burden, solar geometry, and TCWV.

Figure 3 depicts the influence of AOD_550nm on DNI and GHI under various SZA (Figure 3a) and TCWV (Figure 3b) cases. The solar irradiance simulations were performed via the RTM LibRadtran package [69,70]. DNI is apparent more sensitive than GHI as AOD_550nm increases, which is consistent with the study conducted by Ruiz-Arias and Gueymard [71], becoming a very influential parameter for feeding the MLAs. Setting a constant water content (Figure 3a), the solar irradiance drops substantially as SZA increases. For instance, for SZA of 10° (70°), the DNI decreases from ~1000 Wm⁻² (~300 Wm⁻²) to ~350 Wm⁻² (~10 Wm⁻²) at higher AOD_550nm (Figure 3a). Therefore, it is essential to include information on solar geometry in MLAs. For example, if a DNI of 600 Wm⁻² is considered (Figure 3a, under constant TCWL level), this value can be simulated from the RTM from various AOD values varying from 0.02 (for SZA = 50°) to 0.80 (for SZA = 10°). Based on Figure 3b, TCWV affects DNI to a lesser extent than SZA. For instance, if DNI equals 600 Wm⁻², without knowing TCWV, the AOD_550nm ranges from 0.4 (for TCWV = 4 cm) to 0.75 (for TCWV = 0 cm). Therefore, clear-sky solar irradiance is predominantly influenced by solar geometry, aerosols, and water vapor. Therefore, all three variables jointly provide essential information for retrieving AOD_550nm.

3.3. Machine Learning Algorithms

In this study, five machine learning models with different prediction mechanisms, including linear, distance-based, and tree-based approaches, are utilized. This setup establishes a benchmarking task aimed at evaluating and comparing the effectiveness of these ML algorithms in simulating AOD. The five MLAs used in this study are Light Gradient Boosting Machine (LGBM: tree-based), Random Forest (RF: tree-based), Multivariate Adaptive Regression Splines (MARS: linear-based), K-Nearest Neighbors (KNN: distance based), and Artificial Neural Network (ANN: linear based). For the sake of brevity, a comprehensive description of each MLA is presented in Appendix A.

4. Results

4.1. Direct Normal Irradiance vs. Global Horizontal Irradiance as Model Input

This sub-section investigates the performance of the applied methodology using various ML models in terms of the direct and global clearness index. Two distinct scenarios are examined: Scenario 1 considers the three input parameters (see Section 3.2), using the kt_dir as a solar irradiance component. Conversely, in Scenario 2, kt_dir is substituted by kt_glob. Figure 4 depicts the mean absolute error (MAE) and the relative MAE between Scenario 1 and 2, using solely the retrievals of the “optimal” MLA, compared to AERONET measurements. More specifically, the “optimal” MLA is jointly assessed using the MAE, root mean square error (RMSE), and correlation coefficient (R) (Figure S2), assuming an equal contribution for each statistical metric. The “best” MLA at each station is presented in Table S1. The linear correlation between MLA-AOD and AERONET individually for each station per KG climate class is depicted in Figures S3–S7.

The selection of solar irradiance type as an input feature significantly affects the model performance. Scenario 1 reveals the optimal model performance as depicted in Figure 4. The MAE (relative MAE) values ranges are 0.007–0.037 (7.9–35.9%) and 0.018–0.64 (18.1–66.3%) for Scenario 1 and 2, respectively, with the corresponding metrics for each station and Scenario 2 being approximately twice as high. In terms of R, Scenario 1 consistently outperforms Scenario 2, extending from 0.50 to 0.99 (Figure S2). Both scenarios show that the highest RMSEs are recorded for XIA. This can be attributed to the high aerosol variability (Figure S4) in that area due to the high air pollution levels due to local fossil fuel burning from automobiles, the transported polluted plumes from nearby regions, and the transport of dust particles [72]. Despite the relatively high RMSE values in XIA, the relative errors are lower than 10%. Similar results were observed for other stations with high AOD variability such as TAM and SBO (Figure S4).

The difference in the observed performance between the two scenarios is primarily attributed to the impact of aerosol load on GHI and DNI. As shown in Section 3.2 and previous studies [71,73], DNI is highly sensitive to aerosol burden improving the performance of MLAs when it is used as an input parameter.

The performance among MLAs is quite similar at each station, highlighting the model’s versatility (Figure 5 and Figure S2). Nevertheless, for both scenarios, ANN shows the highest number of appearances (Figure S2). For Scenario 1, the majority of the stations had the lowest MAE (Figure 5) values for ANN (12/26), followed by MARS (6/26) and LGBM (5/26). Regarding RMSE (Figure S2c), again, most of the stations recorded their lowest values for ANN (14/26), followed by MARS (8/26) and RF (4/26). Similar results are revealed for the stations in Scenario 2, with ANN being the MLA with the lowest values either for MAE or RMSE (Figure S2). Altogether, ANN and MARS are the most frequently “optimal” MLAs with the most appearances (22 out of 26 of the stations). The following sections consider the input parameters from Scenario 1 (the direct beam solar irradiance showed greater performance against GHI) and the “optimal” MLA.

An interesting finding of the presented methodology is the adequate performance of MLAs at various climates, regions, and aerosol types (Figure 4 and Figure S2). For example, at Tamanrasset (TAM) station, located in Southern Algeria, the main aerosol composition consists of mineral dust particles originating from the Sahara Desert, with a relatively high average AOD over time (here, AOD_550nm = 0.17) [74,75]. Those aerosol particles are quite absorbing. Although the applied methodology is blind regarding the aerosol size or absorptivity, it can retrieve AOD_550nm with reasonable accuracy (Scenario 1; MAE = 0.015–0.037, R = 0.95–0.99). On the other hand, Carpentras, France (CAR), is placed near industrial and urban areas with fine-mode particles as the primary aerosol size. According to Logothetis et al. [75], CAR is mainly affected by fine-absorbing (black carbon) and non-absorbing (sulfate and nitrate) particles. The applied methodology proved to be adequate for accurately retrieving AOD_550nm (Scenario 1; MAE = 0.007–0.011, R = 0.97–0.98) over CAR.

MLA-AOD retrieval performance is also investigated for stations with mixed aerosol conditions, including both fine and coarse aerosol modes. For example, Sede Boker (SBO), Israel, is frequently affected by mineral dust particles emitted by nearby desert areas and fine-mode absorbing aerosols from local sources [75]. The relevant errors in AOD_550nm range from 0.018 to 0.021 (Figure 5) and the correlation coefficient ranges from 0.94 to 0.96 for the various MLAs in Scenario 1 (Figure S2).

The overall performance of MLA-AOD retrievals, including the “optimal” MLA and Scenario 1, is investigated against AERONET in Figure 6. MLAs document a high coefficient of determination (R² = 0.97), MAE of 0.01, and RMSE of 0.02 (Figure 6), highlighting the adequate performance of MLA-AOD methodology compared to AERONET. The majority of the AOD values were lower than 0.5, with the MLA-AOD performance being adequate (see the second plot in Figure 6). For high AOD values (>0.5), the MLA-AOD retrievals have the tendency to underestimate the AERONET measurements.

4.2. The Effect of Aerosol Properties on MLAs Retrieval Performance

The relationship between the differences of the MLA-AOD and observed AOD and the AOD magnitude is investigated (Figure 7). More specifically, the error is drawn against various AOD bins based on AERONET observations (Figure 1).

Based on Figure 7a, the ML-based retrievals encompass higher deviations at higher AOD_550nm. For AOD_550nm lower than 0.4, the mean differences are lower than 0.02 regardless of the underlying climate conditions. The bin frequency at each CC for AOD_550nm < 0.4 is relatively high, including at least 1000 observations. Within high AOD_550nm ranges (>0.4), the deviations of ML estimations are relatively low (between ±0.05), highlighting the model’s high accuracy for extreme AOD_550nm values. In addition, for high AOD_550nm ranges), the MLA-AOD revealed the tendency to underestimate the AERONET AOD. As expected, the frequency of AOD_550nm at high aerosol classes is low, since high AOD_550nm conditions are found in or close to deserts and heavy pollution areas. Additionally, the station density is low in the equatorial and polar climate clusters.

The relative differences, using the AERONET AOD as a skill reference, are calculated to investigate the relationship between the AOD_550nm deviations and their magnitudes (Figure 7b). For all AOD_550nm cases, the relative differences are within 0.20 (Figure 7b). The relatively high documented variability at high AOD_550nm bins (Figure 7a) is reduced when relative differences are used (Figure 7b).

Based on the discussion above, MLAs seem to accurately reproduce the observed AOD in regions with different aerosol characteristics, even if information about the type, size, and absorptivity is absent. The dependency of the MLA-AOD retrieval performance on aerosol properties indicators of aerosol size (a_440–870nm) and absorptivity (SSA_440nm) is investigated (Figure 8) by calculating the error between the MLA-AOD and AERONET AOD within specific a_440–870nm and SSA_440nm bins based on AERONET observations. It is worth mentioning that aerosol properties are not incorporated as independent parameters in the MLA-AOD, since the data availability for SSA is limited in the AERONET sites, and also to keep the MLAs as simple as possible. Regarding aerosol size (Figure 8a), the mean differences are lower than 0.02 regardless of the a_440–870nm bin. Relatively higher variability and mean error are observed for coarse (a_440–870nm < 0.4) and fine particles (a_440–870nm > 2.0), respectively. Regarding aerosol absorptivity (Figure 8b), the mean differences are also lower than 0.02 regardless of the SSA_440nm values. A relatively higher mean error is observed for highly absorbing aerosols (SSA_440nm ≤ 0.85).

Nevertheless, the incorporation of secondary parameters (e.g., single scattering albedo) regarding aerosols into the MLA design could increase the performance of MLA. However, the limited availability of this type of aerosol information supports the decision to evaluate simple MLA with easily found independent predictors such as GHI or DNI, TCWV, and SZA.

4.3. MLA-AOD vs. MODIS

In this section, a comparative analysis between the MLA-AOD and MODIS AOD (both Terra and Aqua) is performed. As described in Section 4.1, only the “best” MLA for each station is used for intercomparison purposes. Figure 9 depicts the RMSE values for the MODIS and MLA-AOD (for Scenario 1, see Section 4.1) against AERONET data. The MLA-AOD outperform MODIS at nearly all stations (except BAR). MODIS RMSEs vary between 0.03 (BAR) and 0.18 (SBO) (Figure 9). Regarding MLA-AOD, RMSE ranges from 0.01 (CAR) to 0.06 (XIA) for the MLA AOD_550nm.

The MLA-AOD and MODIS are linearly compared against AERONET in Figure 10. MLA-AOD documents a high coefficient of determination (R² = 0.98) and RMSE of 0.02 (Figure 10), while the MODIS Aqua and Terra datasets reveal lower R² (0.70) and higher RMSE (0.1) values (Figure 10). The MAE values were 0.02 and 0.07 for MLA-AOD and MODIS, respectively. The comparison between the proposed methodology and MODIS AOD revealed the adequate performance of the MLA-AOD methodology.

Figure S8 presents the daily variability of AERONET, ML-AOD, and MODIS AOD_550nm for specific months and stations. For the TAM station (Figure S8a), all three AOD products have the same pattern for all days in May 2015. The MODIS performance in this case is adequate. Nevertheless, the MLA-AOD provides very accurate predictions, which are almost identical to AERONET measurements. For the XIA station (Figure S8b), the AOD variability in March 2015 was between 0.2 and 1.3. Again, all products follow the daily AOD pattern, but retrievals from the proposed methodology are more accurate, with MODIS having the tendency to overestimate the AERONET measurements. Similar results are observed at the other stations (Figure S8c–f), revealing that the MLA-AOD methodology captures the daily variability of AERONET measurements against MODIS with higher accuracy.

4.4. MLA-AOD vs. Reanalysis Products

Globally, the primary source of aerosol properties originates from spaceborne remote sensing instruments (e.g., MODIS in the previous section). However, aerosol properties can be obtained through aerosol reanalysis datasets at a good temporal resolution. AOD_550nm from MLAs is compared against MERRA-2 and CAMSRA products. As shown in Section 4.1, the “best” MLA for each station is selected for intercomparison purposes.

Figure 11 illustrates the results obtained for each reanalysis dataset and the applied methodology (for Scenario 1, see Section 4.1 against AERONET data. The MLA-AOD outperform MERRA and CAMSRA at all stations and climate types. Regarding MERRA-2, RMSE varies between 0.03 (SMS) and 0.3 (XIA) (orange bars in Figure 11). Accordingly, RMSE ranges from 0.01 (CAR) to 0.1 (XIA) for the MLA AOD₅₅₀ (blue bars in Figure 11). CAMSRA RMSE values lie within 0.04 (SMS) and 0.24 (XIA), providing higher deviations than MERRA-2 in 23 out of 27 stations. It is worth mentioning that MLA-AOD follows the error patterns of the reanalysis dataset but with significantly lower magnitudes. In the majority of the stations, the MLA-AOD approach reveals RMSE values 2–3 times lower than reanalysis datasets compared to AERONET retrievals.

The former findings are consistent with Gueymard and Yang’s [50] study. As they reported, the greater performance of MERRA-2 against CAMSRA can be attributed to the different assimilated remote sensing datasets. Apart from the data derived from spaceborne instruments, MERRA-2 also assimilates the AERONET measurements during 1999–2014. The reanalysis and MLA-AOD are compared against AERONET in Figure 12. All stations are used independently based on their climate type. MLAs document a high coefficient of determination (R² = 0.97) and RMSE of 0.02 (Figure 12), while the reanalysis datasets reveal lower R² (MERRA-2 = 0.72, CAMSRA = 0.70) and higher RMSE (MERRA-2 = 0.07, CAMSRA = 0.09) values (Figure 12).

Consequently, the presented findings verify that MLAs have the potential to retrieve AOD from DNI, being unaware of any other aerosol characteristics during the observation period, resulting in more accurate data than those provided by CAMSRA or MERRA-2.

5. Summary and Conclusions

In this study, a machine learning-based approach (MLA-AOD) is evaluated to estimate AOD at 550 nm using solar irradiance data from ground-based observations, water vapor column (taken from the CAMS reanalysis dataset) and solar geometry. The data are imported into five MLAs with different theoretical prediction mechanisms.

Sensitivity analysis between GHI and DNI is conducted to interpret which solar component (GHI and DNI) is the most appropriate for AOD predictions. The algorithms’ performance provides errors up to two times lower when DNI is used against GHI. Among the MLAs, the ANN provided the highest retrieval performance at most stations. The MLA-AOD methodology performance is also investigated under various climates. For low (<0.4) AOD values, all climate classes showed errors around 0.02, while for higher AOD records, the respective biases slightly increased, tending to underestimate the AERONET AOD. Overall, the MLA-AOD methodology documented a high R² = 0.97, MAE of 0.01, and RMSE of 0.02 compared to AERONET. In addition, MLA-AOD retrievals were accurate in regions with different aerosol optical characteristics, even if information about the type, size, and absorptivity was absent during the model training. The proposed methodology was compared to MODIS AOD, outperforming the satellite data at nearly all stations. The MLA-AOD retrievals were also compared to MERRA-2 and CAMSRA reanalysis products. At all sites, the MLA-AOD outperformed reanalysis data, reaching R² = 0.97 and an RMSE = 0.02.

Considering the great performance of MLAs when reproducing AOD, the results of the proposed methodology could be used to fill and extend either the existing or missing AOD time series derived from ground-based, satellite, and reanalysis datasets emerging an alternative and promising technique for AOD retrieval. Moreover, MLA-AOD values can be used in numerous scientific disciplines, including radiative transfer modeling, aerosol climatology, climate change, as well as air pollution applications.