Multi-layer retrieval of aerosol optical depth in the troposphere using SEVIRI data: a case study of the European continent

Multi-layer aerosol optical depth (AOD) estimation with sufficient spatial and temporal resolution is crucial for effective aerosol monitoring, given the significant variations over time and space. While ground-based observations provide detailed vertical profiles, satellite data are essential for addressing the spatial and temporal gaps. This study utilizes profiles from the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) and data from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) to estimate vertical AOD values at 1.5, 3, 5, and 10 km layers. These estimations are achieved with spatial and temporal resolutions of 3 km × 3 km and 15 min, respectively, over the European troposphere. We employed machine learning models – XGBoost (XGB) and random forest (RF) – trained on SEVIRI data from 2017 to 2018 for the estimations. Validation using CALIOP AOD retrievals in 2019 confirmed the reliability of our findings, emphasizing the importance of wind speed (Ws) and wind direction (Wd) in improving AOD estimation accuracy. A comparison between seasonal and annual models revealed slight variations in accuracy, leading to the selection of annual models as the preferred approach for estimating SEVIRI multi-layer AOD values. Among the annual models, the XGB model demonstrated superior performance over the RF model at all four layers, yielding more reliable AOD estimations with R² values of 0.99, 0.97, 0.98, and 0.98 for the four layers from low- to high-altitude layers. Further validation using data from European Aerosol Research Lidar Network (EARLINET) stations across Europe in 2020 indicated that the XGB model still achieved better agreement with EARLINET AOD profiles, with R² values of 0.86, 0.80, 0.75, and 0.59 and RMSE values of 0.022, 0.012, 0.015, and 0.005. We performed a qualitative validation of multi-layer AOD estimations by comparing spatial trends with CALIOP AOD retrievals for SEVIRI pixels on four dates in 2019, showing strong agreement across varying AOD levels. Additionally, the model successfully estimated AOD at 15 min intervals for two real events – a Saharan dust plume and the Mount Etna eruption – revealing consistent physical characteristics, including long-range transport in the upper layers and a gradual increase in AOD from lower to higher tropospheric layers during volcanic events. The results demonstrate that the proposed method facilitates comprehensive monitoring of AOD behavior throughout the four vertical layers of the troposphere, offering important insights into the dynamics of aerosol occurrence.

1 Introduction

Aerosols are recognized as significant contributors to air pollution, climate change, and the modification of solar and thermal infrared radiation absorption and scattering (Hyslop, 2009; Pope et al., 2019; Li et al., 2022). Understanding aerosol behavior in the troposphere is vital for enhancing atmospheric models and refining monitoring techniques. Aerosol optical depth (AOD) is a critical parameter for quantitatively estimating aerosol concentration and its optical properties. Recent research emphasizes the importance of multi-layer retrieval of AOD in reducing uncertainties associated with aerosol characterization (Wang et al., 2018; Rogozovsky et al., 2021; Gupta et al., 2021; Rogozovsky et al., 2023). Additionally, an in-depth investigation of the multi-layer distribution of aerosol properties within the troposphere is essential for elucidating aerosol transport mechanisms, facilitating source identification, and improving atmospheric dynamics models. This understanding ultimately enhances the accuracy of simulations related to long-range aerosol transport (Chen et al., 2023). Vertical AOD retrieval can be conducted through ground-based observations or inferred from remote sensing data. Ground-based lidar networks, like the European Aerosol Research Lidar Network (EARLINET), provide detailed insights into aerosol characteristics by offering vertical profiles of optical properties, enabling high-resolution, multi-layer AOD retrieval through precise quantification of aerosol loading across distinct atmospheric layers (Bösenberg et al., 2001, 2003). While these observations offer detailed vertical information, their sparse nature necessitates supplementation with satellite observations. Satellite lidar remote sensing emerges as the primary method for capturing global temporal and spatial variations in aerosol profiles. The Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP), on board the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) satellite launched in 2006, offers distributions of aerosols and clouds, along with their geometrical and optical properties. Multi-layer AOD values are retrieved using level 2 aerosol extinction profiles at both 532 and 1064 nm, where the aerosol extinction profiles are determined from backscatter measurements (Winker et al., 2004, 2006, 2007). However, the CALIOP sensor encounters challenges in achieving adequate spatial and temporal coverage, with limitations in daily and global resolution (16 d temporal resolution and 5 km profile distance).

Recent advancements have sought to overcome these limitations through the use of passive satellite sensors with varying temporal resolutions, such as the Tropospheric Monitoring Instrument (TROPOMI), which provides near-daily global coverage with a spatial resolution of 3.5 × 7 km (improved to 3.5 × 5.5 km in 2019) and was launched in 2017 on the Sentinel-5P satellite (Veefkind et al., 2012); the Earth Polychromatic Imaging Camera (EPIC), offering a continuous daytime view every 60 to 100 min with a spatial resolution of about 8 × 8 km since its launch on 11 February 2015, on board the Deep Space Climate Observatory (DSCOVR) satellite (Marshak and Knyazikhin, 2017); the Global Ozone Monitoring Experiment-2 (GOME-2) on the Meteorological Operational satellite program (MetOp-C), with a 3 d revisit cycle and a spatial resolution of approximately 40 × 40 km since 2018; and the Moderate Resolution Imaging Spectroradiometer (MODIS), on board Terra (launched in 1999) and Aqua (launched in 2002), providing daily global coverage with spatial resolutions ranging from 0.25 to 1 km (Lyapustin et al., 2011).

Relevant research focuses on various methods specifically aimed at retrieving aerosol layer height (ALH) rather than AOD at different altitudes. One prominent method, oxygen (O₂) A- and B-band absorption spectroscopy, utilizes the differential absorption of sunlight by O₂ molecules at different altitudes (Zeng et al., 2018; Xu et al., 2017, 2019). Elevated aerosol layers scatter sunlight back to space, shortening the atmospheric path length and decreasing O₂ absorption. By analyzing spectral characteristics in the O₂ A and B bands, researchers infer ALH. However, retrieval sensitivity is enhanced over darker surfaces and higher AOD, making it challenging over bright surfaces or under low aerosol loading. For instance, Nanda et al. (2020) employed TROPOMI observations with an optimal estimation scheme in the O₂ A band, assuming a uniformly distributed aerosol layer. Similarly, the algorithm developed using EPIC/DSCOVR data leverages atmospheric window bands and differential optical absorption spectroscopy (DOAS) ratios, integrating MODIS and GOME-2 surface reflectance data. For retrievals over vegetated areas, the algorithm favors the O₂ B band due to its lower surface reflectance (Xu et al., 2019). Another study combined O₂ A- and B-band data from the Scanning Imaging Absorption Spectrometer for Atmospheric Chartography (SCIAMACHY) and GOME-2 for enhanced ALH sensitivity, especially near boundary layers (Hollstein and Fischer, 2014).

An additional retrieval method, stereoscopic techniques – employed by the Multi-angle Imaging SpectroRadiometer (MISR), launched in 2000 – utilize multi-angle observations to geometrically determine plume heights. MISR offers a spatial resolution of approximately 275 m and a temporal resolution of around once every 7 d, making it especially useful over reflective surfaces, as it relies on geometric data rather than surface reflectance (Muller et al., 2002; Zakšek et al., 2013; Fisher et al., 2013; Val Martin et al., 2018).

Passive satellite-based ALH retrieval techniques, while offering global coverage, often simplify the aerosol vertical distribution by assuming a single homogeneous layer (Zeng et al., 2018; Xu et al., 2017, 2019). This simplification can lead to inaccurate representations of complex aerosol profiles, especially in cases of multi-layered events. In addition, these passive satellite-based methods face further constraints due to the low spatial resolution of instruments like EPIC and GOME-2, as well as low temporal resolution of sensors such as TROPOMI, GOME-2, and MISR. These constrains on resolution reduce the effectiveness of these retrievals in capturing fine-scale, rapidly evolving aerosol distribution events, such as smoke plumes from fires.

Other studies by Pashayi et al. (2023, 2024) have introduced seasonal and seasonal-independent machine learning models for AOD retrievals in multiple layers. These models seek to investigate the relationship between MODIS observations and CALIOP AOD for retrieval of multiple-layer AOD values at a spatial–temporal resolution corresponding to the MODIS AOD product. This analysis focuses specifically on the Persian Gulf region. The researchers subsequently analyze their findings using CALIOP AOD retrievals across multiple vertical layers. Although these studies have advanced the retrieval of AOD across multiple layers, the constraint of MODIS's daily temporal resolution remains a significant limitation (Wei et al., 2020).

Geostationary satellites, such as Himawari-8 (launched in 2014), which is equipped with the Advanced Himawari Imager (AHI; Da, 2015); the Geostationary Operational Environmental Satellite (GOES; launched in 2016) equipped with the Advanced Baseline Imager (ABI; Kalluri et al., 2015); and the Meteosat geostationary satellites featuring the Spinning Enhanced Visible and Infrared Imager (SEVIRI; launched in 2002; Pasternak et al., 1994), provide sub-hourly, high-resolution observations that significantly enhance global aerosol monitoring capabilities across diverse regions (Schmit et al., 2018; Zhang et al., 2019; Ge et al., 2018; Tang et al., 2019; Zawadzka-Manko et al., 2020; Witthuhn et al., 2020; Kocaman et al., 2022; Ceamanos et al., 2023). Notably, SEVIRI offers high temporal and spatial resolutions, presenting valuable opportunities to expand aerosol datasets for Europe (Stebel et al., 2021; Ajtai et al., 2021). Consequently, utilizing observations from these satellites enables the multi-layer retrieval of AOD, effectively addressing the limitations associated with temporal and spatial resolution presented in previous studies.

The retrieval of multi-layer AOD values from passive satellite observations typically entails two primary approaches: physically based (Seidel et al., 2012; Lipponen et al., 2018; Amini et al., 2021; Mehta et al., 2022) and data mining approaches (Radosavljevic et al., 2010; She et al., 2020; Chen et al., 2022). The physically based approach relies on established principles of aerosol behavior, utilizing models derived from physical laws to retrieve AOD values. This approach often involves simplifications and assumptions, such as treating the atmosphere as a single aerosol layer in most of the passive satellite-based ALH retrieval algorithms previously mentioned. While this assumption is necessary for practical implementation, it can introduce uncertainties and limit the accuracy of retrievals, particularly in complex scenes. Additionally, physically based methods are sensitive to surface reflectance. Over bright surfaces, the contribution of surface reflection to top-of-atmosphere (TOA) radiance can dominate, making it challenging to extract a clear aerosol signal, especially for low aerosol loading. This limitation underscores the need for accurate surface reflectance characterization (Xu et al., 2017, 2019; Nanda et al., 2020).

In contrast, data mining approaches offer a promising alternative by harnessing large datasets and employing learning-based algorithms to discern patterns and relationships within complex aerosol systems. Machine learning methods as a data mining approach have the potential to learn complex relationships between AOD and other atmospheric and surface variables. These methods can capture non-linear dependencies and potentially handle multi-layer scenarios more effectively than physical models. This advantage is particularly relevant for AOD retrieval in diverse and heterogeneous environments. Additionally, integrating data from multiple sources, such as meteorological, land cover, temporal, and location data, can provide more comprehensive information for a time- and location-based AOD retrieval, particularly for multi-layer scenarios (Chen et al., 2020; Lee et al., 2022; Berhane et al., 2024).

In this study, we introduce a model for sub-hourly multi-layer AOD retrieval over the European continent troposphere by integrating SEVIRI-based information with CALIOP aerosol profile products. To achieve this, two well-established machine learning models – XGBoost (XGB) and random forest (RF) – were utilized for retrieving AOD values in four distinct layers, approximately every 15 min, with a spatial resolution of 3 km × 3 km. The four tropospheric layers analyzed in this study are 0–1.5, 1.5–3, 3–5, and 5–10 km, denoted as AOD_1.5, AOD₃, AOD₅, and AOD₁₀, respectively. The selection of these layers for multi-layer AOD retrieval is based on the distinct aerosol transport mechanisms observed at these altitudes. The 0–1.5 km layer captures aerosols from local sources transported upward by updrafts from the cloud base, a process called pumping. The 1.5–3 km layer, where thermal bubbles often initiate, allows for the examination of aerosols, potentially from mid-range sources, which are lifted into the cloud with the rising bubble. The 3–5 km layer captures aerosols transported over longer distances that enter the cloud through entrainment at the cloud edges as the bubble ascends. The 5–10 km layer is designed to capture the influence of long-range-transported aerosols on cloud properties at higher altitudes. This multi-layer approach enables analysis of how local to long-range aerosol transport contributes to aerosol–cloud interactions (Zhang et al., 2021; Lebo, 2014; Marinescu et al., 2017).

To train and validate two machine learning models, we employed AOD data retrieved from the CALIOP aerosol product and EARLINET stations distributed across Europe. Model performance was qualitatively evaluated by analyzing its response to two notable aerosol events: a significant dust intrusion from 13 to 18 March 2022 and a volcanic eruption on 14 August 2023. These events offered valuable case studies to assess the model's capability in detecting and characterizing distinct aerosol signatures across these four layers. We organized the rest of the paper as follows: Sect. 2 provides a comprehensive overview of the dataset employed, while Sect. 3 details the necessary preprocessing steps and retrieval methodology. Subsequently, Sect. 4 delves into the discussion of the vertically retrieved AOD results, followed by conclusions outlined in Sect. 5.

2 Study area and data source

2.1 Study area

The study area encompasses a significant portion of the European troposphere, spanning 35 to 71° N and 7° W to 70° E, covering approximately 10.18×10⁶ km². Despite its relatively small land area, Europe exhibits a diverse geographical landscape and complex atmospheric dynamics. Urban centers in Europe face persistent air pollution issues due to industrial activities and vehicular emissions, compounded by the effects of climate change. Various aerosol types, originating from industrial processes, transportation, biomass burning, and natural events, significantly impact air quality, weather patterns, and climate dynamics across the continent. Long-range transport of aerosols, particularly from sources in Africa, such as Saharan dust storms, underscores the interconnectedness of atmospheric processes across continents and emphasizes the necessity of international cooperation in addressing air pollution and environmental challenges.

2.2 Data source

2.2.1 SEVIRI

Meteosat Second Generation (MSG) constitutes a series of four satellites managed by the Exploitation of Meteorological Satellites (EUMETSAT) and has been operational since 2004. Originally designated as MSG1 to MSG4, these satellites were subsequently rebranded as Meteosat-8 to Meteosat-11, respectively. The primary instrument on board these satellites is the Spinning Enhanced Visible and Infrared Imager (SEVIRI), a radiometer equipped with 11 spectral channels spanning the visible to the infrared spectrum. These include the visible (0.6 and 0.8 µm) channels, as well as a near-infrared (1.6 µm) channel, and provides a spatial resolution of about 3 km at the sub-satellite point and a high-resolution visible (HRV) channel offering a finer spatial resolution of 1 km at nadir. Strategically centered at various wavelengths, the thermal channels of SEVIRI include 6.2 and 7.3 µm (targeting strong water vapor absorption); 8.7, 10.8, and 12.0 µm (window channels); and 9.7 µm (for ozone absorption) and 13.4 µm (for carbon dioxide absorption). This operational system delivers full-disk Earth data, while the rapid-scan service focuses on observing the upper part of the Earth's disk, covering Europe and North Africa, with a repetition time of 15 min (Schmetz et al., 2002; Zawadzka and Markowicz, 2014). In our study, we primarily utilize SEVIRI data from Meteosat-11, the fourth and final flight unit of the MSG program, which was launched on 15 July 2015. Meteosat-11 currently operates in geostationary orbit, positioned 36 000 km above the Equator. Its coverage extends over Europe, Africa, and the Indian Ocean, spanning −81 to 81° longitude and −79 to 79° latitude. Figure 1 provides a visualization of the coverage area of SEVIRI (© EUMETSAT 2024).

Figure 1 The area covered by the SEVIRI instrument. Source: © EUMETSAT 2024 (https://data.eumetsat.int/data/map/EO:EUM:DAT:MSG:HRSEVIRI, last access: 20 January 2025).

2.2.2 CALIOP

The CALIOP instrument plays a pivotal role in the CALIPSO satellite, launched in April 2006 with the primary objective of reliably delivering high-resolution vertical profiles of global aerosol properties via an active sensing technique. Functioning as a polarization-sensitive lidar, CALIOP measures the depolarization ratio, serving as a degree of particle irregularity. CALIOP is specifically designed to observe aerosol optical properties during both the day and the night, focusing on vertical layers at wavelengths of 532 and 1064 nm. Its level 2 algorithm not only provides information on aerosol optical characteristics like the particle depolarization ratio and color ratio but also retrieves extinction coefficients. Notably, CALIOP data offer a temporal resolution of approximately 16 d, capturing insights into aerosol dynamics over time. Sampling occurs at intervals of 333 m along the orbital track, maintaining a vertical resolution of 60 m from altitudes of −0.5 to 20 km and 180 m from 20 to 30 km within the vertical profile (Winker et al., 2004, 2006, 2007). For this study, we employed CALIOP level 2 version 4.2 aerosol profile products from 2017 to 2019 to estimate multi-layer AOD values within the defined study region.

2.2.3 MODIS land cover data

In this research, we leveraged land cover (LC) data spanning 2017 to 2019, with a spatial resolution of 1 km, sourced from the global MODIS products (MCD12Q1 V6) covering Europe. These data, derived from both Terra and Aqua satellites, provide comprehensive land cover types annually from 2001. The dataset encompasses six classification schemes, elucidated in the downloadable user guide available at https://ladsweb.modaps.eosdis.nasa.gov/ (last access: 23 October 2023). Each MCD12Q1 version 6 hierarchical data format 4 (HDF4) file comprises layers for land cover types 1–5, land cover properties 1–3, land cover property assessments 1–3, land cover quality control (QC), and a land water mask (Sulla-Menashe and Friedl, 2018). Our study specifically focuses on the first classification scheme, the Annual International Geosphere-Biosphere Program (IGBP) classification.

2.2.4 Meteorological data

Meteorological data were acquired from the European Centre for Medium-Range Weather Forecasts (ECMWF) dataset, accessible at https://cds.climate.copernicus.eu/ (last access: 24 October 2023). ECMWF has been actively operational in real-time seasonal forecast systems since 1997, providing access to standard meteorological data. This dataset comprises two distinct sets of data (Copernicus Climate Change Service Climate Data Store, 2021). Firstly, version 2 of the Integrated Global Radiosonde Archive (IGRA) from 1978 integrates global radio sounding containing temperature, humidity, and wind data from various sources. The dataset is presented in the form of a global grid with a conventional grid resolution of 0.25° × 0.25°. Compared with previous-generation products, the temporal resolution has been increased from 6 to 1 h, enabling the study of diurnal variations in the troposphere. Secondly, the Radio Sounding HARMonization (RHARM) homogenized dataset offers adjusted values for temperature, relative humidity, and wind. RHARM effectively eliminates systematic effects, such as variations in measurement sensors, biases induced by solar radiation, calibration drifts, station relocations, and other factors, across 700 IGRA radiosonde stations and ship-based radio soundings. RHARM includes twice-daily (00:00 and 12:00 UTC) radiosonde data at mandatory and standard levels, featuring essential parameters like air temperature (T, K), air pressure (P, Pa), wind speed (Ws, m s⁻¹), and wind direction (Wd, degrees from north). For this study, the global grid dataset is utilized over the European continent from 2017 to 2019, as depicted in Fig. 2.

Figure 2 Map depicting the ECMWF stations for meteorological data measurements.

2.2.5 EARLINET

EARLINET, established in the year 2000 (Bösenberg et al., 2001, 2003), originated as a research project funded by the European Commission within the framework of the Fifth Framework Programme. The primary objective of EARLINET is to generate profiles of aerosol optical properties, thereby constructing an expansive, quantitative, and statistically robust database for the continental-scale distribution of aerosols. This initiative aims to enhance network operations, facilitate research on aerosol-related processes, validate satellite sensor data, advance model development and validation, integrate aerosol data into operational models, and compile a comprehensive climatology of aerosol distribution. Currently, the network comprises 30 active stations, with the majority equipped with Raman lidar featuring depolarization channels. These Raman-lidar-operating EARLINET stations typically provide profiles of aerosol extinction and backscatter coefficients without relying on significant assumptions. Figure 3 illustrates the distribution of EARLINET stations over the study area.

Figure 3 Map depicting currently active EARLINET stations. The red stars indicate the geographical distribution of EARLINET lidar stations used in this study.

3 Methodology

As noted, hyperspectral measurements in the oxygen bands enable aerosol vertical distribution retrieval by analyzing photon path length changes due to scattering at different altitudes. SEVIRI's spectral bands, however, are primarily designed for cloud and land surface observations and do not specifically cover the oxygen bands. The SEVIRI bands closest to the oxygen bands are B₁ (635 nm) and B₂ (810 nm), which are in the visible spectrum and respond to scattering by vertically distributed aerosols. The near-infrared and shortwave infrared (SWIR) bands (B₃, B₄, B₇, B₉, and B₁₀) are indirectly influenced by aerosol vertical distribution, as the accuracy of AOD retrievals using these bands can be affected by aerosol layering. While these bands may not directly provide vertical profile information, they could yield complementary data that, when combined with other wavelength retrievals, enhance the understanding of aerosol vertical distribution (Wu et al., 2017; Li et al., 2020). B₅ and B₆ offer insights into water vapor profiles, which can be incorporated into aerosol retrieval algorithms. By accounting for water vapor influence, these bands indirectly improve the accuracy of aerosol vertical distribution estimates. The ozone band (B₈) contributes to atmospheric chemistry and aerosol formation insights but does not directly reveal vertical distribution, while B₁₁'s lower scattering efficiency limits its sensitivity to vertical variations in aerosols for direct retrieval. As a result, SEVIRI's bands provide a range of potential avenues for studying aerosol vertical distribution, with both direct and indirect contributions. Although SEVIRI's bands offer valuable data for meteorological observations, such as cloud monitoring, surface temperature, and water vapor, its spectral design is not optimized for detailed monitoring of air quality or climate through atmospheric gases and aerosols, similar to TROPOMI. Thus, physical approaches for detailed multi-layer aerosol retrieval, especially for multi-layer AOD, remain challenging with SEVIRI's current spectral configuration.

Meteorological data significantly influence the vertical distribution of aerosols, with varying impacts depending on aerosol type, transport dynamics, and atmospheric conditions. Wind speed and direction drive both horizontal and vertical aerosol transport, with higher wind speeds over oceans enhancing sea salt aerosol concentrations (Kaufman et al., 1997; Yu et al., 2006; Chin et al., 2007; Tesche et al., 2009). Temperature and pressure also play critical roles; temperature inversions inhibit vertical mixing, trapping aerosols in distinct layers, while convective activity from surface heating mixes aerosols in the boundary layer, creating a more homogeneous distribution. Stable high-pressure systems promote surface accumulation by limiting mixing, whereas low-pressure systems enhance upward transport, extending aerosol atmospheric lifetimes (Tesche et al., 2009). The complexity of these interactions suggests a significant challenge for multi-layer AOD retrieval using physical approaches, as accurate modeling requires accounting for diverse meteorological influences and variations in aerosol type, transport, and vertical distribution.

Geographical location, land cover, and temporal factors significantly influence the vertical distribution of aerosols across Europe. Coastal regions tend to have elevated sea salt aerosols due to ocean surface wind activity, while continental areas, especially in winter, experience higher anthropogenic aerosol concentrations from sources like fossil fuel combustion and industrial emissions. Additionally, the latitude and prevailing wind patterns, such as easterly winds, play a role in the long-range transport of aerosols, affecting distribution both horizontally and vertically. Land cover also contributes to these dynamics: forests emit biogenic volatile organic compounds (VOCs), which can form secondary organic aerosols, while urban and agricultural areas introduce anthropogenic aerosols from activities like traffic, industrial emissions, and fertilizer use.

Temporal variations, including seasonal and diurnal changes, further complicate aerosol distribution. For example, during winter, stable high-pressure systems trap aerosols in the planetary boundary layer (PBL), while in summer, warmer temperatures enhance photochemical activity, leading to increased ozone and sulfate concentrations. Diurnal fluctuations are also evident, particularly in urban areas, where traffic and industrial activities create peaks in anthropogenic emissions during the day.

These combined effects underscore the complexity of aerosol behavior, emphasizing the necessity for an approach that integrates all relevant variables and effectively captures their interactions and influence on vertical aerosol distribution and multi-layer AOD retrieval. Machine-learning-based methodology, capable of managing large datasets and discerning intricate relationships between these variables, presents a promising solution for accurate multi-layer AOD retrieval. Our proposed model framework for estimating AOD at the mentioned four distinct layers over the European continent troposphere encompasses several sequential steps – data collection, preprocessing, partitioning, regression, and analysis of the performance of each regression model – to ascertain the most accurate one, as illustrated in Fig. 4.

Figure 4 Research framework for developing machine learning models to estimate SEVIRI multi-layer AOD values.

The process commences with data collection, detailed in the preceding section. Subsequently, preprocessing of both input and output data becomes imperative to ensure their suitability for subsequent analysis. The dataset is then partitioned into two subsets: training and testing, a pivotal step in machine learning aimed at assessing model performance and mitigating overfitting. Following data partitioning, various model structures are proposed and developed to capture the intricate relationships within the dataset. This phase entails selecting appropriate algorithms and architectures tailored to the specific task of multi-layer AOD estimation. Finally, the performance of each model is meticulously evaluated using predefined metrics to pinpoint the most accurate and reliable model for AOD estimation across the desired vertical layers. In the subsequent sections, we delve into a detailed examination of each step.

3.1 Preprocessing

To ensure a robust model for estimating AOD values at suitable 3D resolutions, this study integrates data from various sources, including satellites and ground-based observations. To address spatial–temporal sampling disparities, we employ a co-location approach where data from multiple sources, such as satellites and ground-based observations, are matched within a ±30 min time frame and within a 3 km radius of the study area (Kittaka et al., 2011; Redemann et al., 2012; Han et al., 2017; B. Liu et al., 2018). This method harmonizes disparate datasets, enhancing the reliability and comprehensiveness of our analysis. The subsequent preprocessing stages necessary for data refinement and analysis are elaborated upon in the following subsections.

3.1.1 SEVIRI

Utilizing SEVIRI data necessitates a critical preprocessing step involving co-referencing and applying geometric corrections. The Data Tailor tool, accessible at https://www.eumetsat.int/data-tailor (last access: 28 October 2024), serves as a valuable spatial resource introduced in recent years. It simplifies the definition of coordinate systems, image systems, cutting ranges, expected output types, and requisite file extensions for the output data. Estimating AOD values requires the conversion of radiance to reflectance for the SEVIRI reflective bands (VIS06, VIS08, and NIR16) and equivalent brightness temperature for the remaining eight bands. To achieve this, we computed the bidirectional reflectance factor (BRF) for the SEVIRI warm channels using Eq. (1) proposed by the European Organization for the Exploitation of Meteorological Satellites (2012):

$$\begin{array}{}\text{(1)}& r_{{\mathit{\lambda }}_i}={\displaystyle \frac{\mathit{\pi }\cdot R_{{\mathit{\lambda }}_i}\cdot d^{\mathrm{2}}\left(t\right)}{I_{{\mathit{\lambda }}_i}\cdot \cos \left(\mathit{\theta }\right(t,x\left)\right)}},\end{array}$$

where i denotes the channel number (1, VIS06; 2, VIS08; 3, NIR16; 4, HRV), $r_{{\mathit{\lambda }}_i}$ represents the bidirectional reflectance factor (BRF) for channel λ_i, $R_{{\mathit{\lambda }}_i}$ stands for the measured radiance in mW m⁻² sr⁻¹ (cm⁻¹)⁻¹, d(t) signifies the Sun–Earth distance in astronomical unit (AU) at time t, $I_{{\mathit{\lambda }}_i}$ signifies the band solar irradiance for channel λ_i at 1 AU in mW m⁻² sr⁻¹ (cm⁻¹)⁻¹, and θ(t, x) denotes the solar zenith angle in radians at time t and location x. The equivalent brightness temperature (T_b) of a satellite observation is defined as the temperature of a black body emitting the same amount of radiation. Therefore, the brightness temperature follows the form of Eq. (2).

$$\begin{array}{}\text{(2)}& T_{\mathrm{b}}={\displaystyle \frac{C_{\mathrm{2}}{v}_{\mathrm{c}}}{\mathit{\alpha }\log C_{\mathrm{1}}{v}_{\mathrm{c}}^{\mathrm{3}}/\stackrel{\mathrm{‾}}{R}+\mathrm{1}}}-{\displaystyle \frac{\mathit{\beta }}{\mathit{\alpha }}}\end{array}$$

Using the observed radiances $\stackrel{\mathrm{‾}}{R}$ (in mW m⁻² sr⁻¹ (cm⁻¹)⁻¹) and radiation constants C₁ = 2hc² and C₂ = $hc/k$, where c, h, and k represent the speed of light, Planck's constant, and the Boltzmann constant, respectively, the regression coefficients v_c, α, and β are determined through non-linear regression analysis. This analysis is conducted on a pre-calculated lookup table generated for the various SEVIRI channels, as delineated in Table 1 (Tjemkes et al., 2012).

Table 1 Values for the regression parameters.

To further enhance the preprocessing workflow, we applied the SEVIRI cloud mask product to the data. This product categorizes pixels with values 00, 11, and 22, corresponding to clear, partially cloudy, and cloudy conditions, respectively. To integrate the cloud mask into the data, we utilized the following transformation:

$$\begin{array}{}\text{(3)}& \text{Adjusted bands}=\mid \text{Cloud mask value}-\mathrm{2}\mid \cdot B_i.\end{array}$$

This ensures that data from cloudy or partially cloudy pixels are appropriately weighted or excluded, preserving the accuracy of AOD retrieval for clear-sky conditions. By incorporating the cloud mask product, the preprocessing step effectively eliminates biases introduced by cloud contamination, ensuring the reliability of subsequent AOD estimations.

3.1.2 CALIOP

In this study, to mitigate the impact of cloud contamination and retrieval errors on CALIOP AOD retrieval, our screening methods closely follow the guidelines established by Winker et al. (2013). We employ various quality filters to identify and filter aerosol pixels, including cloud and aerosol discrimination (CAD) scores, extinction QC flags, and uncertainty values. Specifically, we utilize a CAD score range outside $[-\mathrm{100},-\mathrm{20}]$ to address uncertainties in cloud–aerosol discrimination, ensuring the selection of cloud-free pixels with high confidence. Additionally, we apply extinction quality control flags with values of 0 and 1 to filter extinction retrievals with high confidence. This includes constrained retrievals utilizing transmittance measurements and unconstrained retrievals where the initial lidar ratio remains unchanged in iterations. Furthermore, we exclusively consider daytime profiles in this study. Uncertainty flags associated with extinction coefficients are employed for data screening. Range bins with an uncertainty flag value of 99.9 km⁻¹ are excluded from the analysis, following the methodology outlined by Winker et al. (2013).

3.1.3 Land cover product

Considering that the original MCD12Q1 product is stored in an HDF and utilizes the sinusoidal projection, several data preprocessing steps are required. These steps encompass format conversion, reprojection, resampling, image mosaicking, and sub-area masking. To execute these tasks, we employ the pyModis free and open-source Python-based library. This tool enables the conversion of MODIS HDF into GeoTIFF format and facilitates the conversion of data projection from SIN to WGS84/UTM. Additionally, it facilitates image mosaicking and subsetting. Moreover, to enable comparison between the MCD12Q1 and SEVIRI datasets, the spatial resolution of MCD12Q1 is resampled at 3 km using the nearest neighbor resampling method. This method preserves the gray values of the original image, unlike bilinear interpolation or cubic convolution interpolation methods, which may alter them.

3.2 Machine learning models and parameter tuning

In this study, our primary objective is to develop a machine learning model to estimate SEVIRI AOD values at various altitudes – 1.5, 3, 5, and 10 km – using CALIOP's vertical profiles across the European continent. We employ two distinct machine learning algorithms, RF and XGB, to train layering models. Both RF and XGB adopt an ensemble approach, which involves constructing and aggregating multiple decision trees (Breiman, 2001; Chen and Guestrin, 2016). In RF, each tree is built using a bootstrap sample of the data, with nodes determined by the best subset of randomly selected predictors (Breiman, 2001). These trees are then averaged to obtain a final ensemble prediction. Conversely, XGB implements the gradient boosting method, where trees are interdependent, as newly trained trees are constructed based on previous trees, incorporating their ability to predict the residuals of prior trees (Chen and Guestrin, 2016). In both RF and XGB, all trained trees are combined to make the final prediction.

We systematically explored various parameter combinations for each machine learning model. Parameters such as the number of decision trees (N_estimators), the number of variables considered for splitting at each node (max_features), and the maximum depth of each decision tree (max_depth) for RF, as well as parameters including the number of gradient boosting rounds or decision trees (n_estimators), minimum sum of instance weight (Min_sample split), maximum depth of each decision tree (max_depth), and minimum number of samples required to be at a leaf node (Min_sample leaf) for XGB, were optimized using a grid search algorithm. This algorithm exhaustively searches through a specified subset of the hyperparameter space. We set up a grid of possible values for each hyperparameter to be tuned, as illustrated in the specific search range column in Table 2. For each combination of hyperparameters in the grid, the algorithm trains the model using the training data and evaluates its performance through cross-validation. The performance of each hyperparameter combination is measured using several specified evaluation metrics. Finally, the combination of hyperparameters that results in the best performance on the validation set is selected, as shown in the optimum value column in Table 2. This optimal set of parameters is then used to train the final model on the entire training dataset. For a comprehensive overview of the optimized parameters, refer to Table 2.

Table 2 The control parameter for tuning the machine learning models.

3.3 Model training and evaluation

Data partitioning is essential for training and evaluating machine learning models, especially when working with time-series data where temporal autocorrelation might bias model performance. In this study, the dataset was partitioned by year to ensure temporal independence between training and testing data, addressing potential autocorrelation issues and enabling robust model evaluation. Specifically, the data from 2017 to 2018 were used for model training, while the 2019 data were reserved exclusively for testing. This approach was chosen following an analysis of the feature distributions of SEVIRI bands (B₁ to B₁₁); P, T, LC, Ws, and Wd; and multi-layer AOD values (AOD_1.5, AOD₃, AOD₅, and AOD₁₀) over the different years, as shown in Fig. S1 in the Supplement. The distribution represented in Fig. S1 reveals consistent patterns between 2017–2018 and 2019, with minimal variation in their shapes. This similarity confirms that the temporal separation does not introduce significant distributional shifts that might impact model generalization. In other words, the model's performance on the 2019 data would provide an unbiased evaluation of its predictive ability. This separation minimizes temporal autocorrelation, ensuring robust and unbiased model assessment.

During the training phase of our machine learning models, we leveraged datasets spanning diverse temporal periods and geographical regions where both SEVIRI and CALIOP data were accessible. However, following this training phase, the algorithms function autonomously, relying solely on SEVIRI data as their input. This advancement enables us to estimate AOD values at four specified vertical layers within each pixel of the SEVIRI dataset, based on a single SEVIRI observation along with its associated meteorological data and land cover data, covering the entire study area.

Evaluation of the multi-layer AOD-estimating models involved statistical metrics such as the coefficient of determination (R²), Pearson correlation coefficient (R), root mean square error (RMSE), and mean absolute error (MAE). The selection of the optimal model was based on higher R² and R values, along with lower RMSE and MAE scores. Additionally, we conducted a validation analysis of estimated multi-layer AOD values with EARLINET AOD profiles on a continental scale to ascertain the model's performance.

4 Results and discussion

In this paper, our primary aim is to develop a machine learning model capable of retrieving AOD across four distinct vertical layers: 1.5, 3, 5, and 10 km. To accomplish this, we utilized two well-established machine learning models, XGB and RF, previously employed in related studies. These models were trained on SEVIRI data spanning the European continent from 2017 to 2019. Our objective was to estimate sub-hourly AOD values, approximately every 15 min, at a spatial resolution of 3 km × 3 km.

To explore the relationship between AOD and potential predictor variables, we conducted a correlation analysis experiment utilizing the Pearson correlation coefficient (PCC; Benesty et al., 2009). Furthermore, we evaluated the influence of land cover and meteorological data as input variables for the machine learning models on estimating multi-layer AOD values from SEVIRI data, with a specific focus on identifying the most optimal model. Moreover, we conducted training and testing of the machine learning models across various temporal scales, including annual and seasonal analyses. Subsequently, we assessed the performance of each model using independent satellite and ground-based AOD profiles, employing evaluation metrics such as R², R, MAE, and RMSE. Finally, multi-layer AOD values for two aerosol events – Saharan dust from 13 to 18 March 2022 and a volcanic eruption on 14 August 2023 – are presented as maps to evaluate the model's ability to detect and characterize distinct aerosol signatures across these four layers. In the subsequent sections, we provide a comprehensive review of the results derived from the aforementioned assessments.

4.1 Validation of estimated AOD with satellite-retrieved AOD

4.1.1 Feature importance

According to established radiative transfer theory (Tsang et al., 1984; Zege et al., 1991), the spectral signal captured by a satellite sensor at the top of the atmosphere is intricately shaped by various factors, including the composition, size distribution, and altitude of aerosols, as well as atmospheric molecules such as water vapor. These factors have a direct impact on the retrieval of AOD values. Consequently, SEVIRI reflectance and brightness temperature across bands 1 to 11 (B₁ to B₁₁) were identified as critical features for this analysis. The relationship between AOD and all candidate features – including spatial features such as latitude (lat) and longitude (long); temporal features including year, month, and day; meteorological data like P, T, Ws, and Wd; and LC – was investigated through a correlation analysis. This analysis, illustrated in Fig. 5, utilized the PCC as the chosen filtering method. The findings underscored that the majority of selected features in this study exhibited significance levels exceeding 1 %.

Figure 5 The importance of input features in the retrieval of SEVIRI multi-layer AOD values, as determined by the PCC.

As illustrated in the PCC results (Fig. 5), aerosol dynamics in the 0–1.5 km layer of the study area are strongly influenced by geographic and temporal factors, characteristic of surface-dominated conditions. High correlations of AOD with long (23.92 %) and lat (22.41 %) highlight the impact of location-specific emissions and regional transport patterns, while significant correlations with day (11.59 %) and month (6.19 %) indicate the role of diurnal and seasonal cycles driven by emissions, meteorological changes, and boundary layer dynamics. Thermal infrared sensitivity also plays a role, as indicated by the strong correlations with B₇ (3.73 %) and B₈ (3.07 %), the highest among SEVIRI bands. This suggests a potential link between surface temperature variations and AOD in the lower atmosphere. Moreover, the correlation with B₁ (3.06 %) points to aerosols interacting with visible light, likely from urban/industrial emissions, biomass burning, and potentially dust. Meanwhile, meteorological factors such as Ws (1.31 %), Wd (1.22 %), T (1.16 %), and P (1.1 %) are relatively weak compared to geographic and temporal factors. Finally, LC (0.54 %) shows the weakest correlation, suggesting a limited direct influence on AOD, potentially masked by stronger influences from emissions, transport, and meteorological factors.

In the 1.5–3 km layer of area of interest, aerosol dynamics remain strongly correlated with long (23.15 %) and lat (21.39 %), though they decrease slightly compared to the 0–1.5 km layer. This suggests a transition toward atmospheric stability and more synoptic-scale patterns. Day (11.29 %) and month (5.86 %) still exhibit notable correlations, reflecting the ongoing impact of diurnal and seasonal cycles on AOD, albeit slightly weaker than in the lower layer. The visible band B₁ (3.86 %) shows the strongest correlation among SEVIRI bands, potentially reflecting a change in aerosol composition or properties compared to the lower layer. This could be due to an increased influence of transported aerosols, potentially with different spectral characteristics. B₇ (3.48 %) remains strongly correlated, suggesting continued sensitivity to thermal characteristics. Enhanced transport influence is evident, as Wd (1.34 %) shows a stronger correlation than Ws (1.29 %), emphasizing the role of transport pathways in AOD distribution at this altitude. Additionally, the correlation with P (1.3 %) points to an emerging link between atmospheric stability and aerosol accumulation in this layer. This aligns with other studies' discussion on PBLH and its role in controlling pollutant dispersion, as PBLH is influenced by factors like temperature and stability. LC correlation increases slightly to 0.63 %, still remaining relatively weak overall.

In the 3–5 km layer of the research area, correlations with long (19.42 %) and lat (17.32 %) decrease further, indicating a reduced influence of local sources and an increased dominance of long-range transport and large-scale atmospheric circulation patterns. Correlations with day (9.87 %) and month (6.1 %) also weaken, suggesting a shift from surface-driven cycles to broader atmospheric processes. The high correlation with B₁ (5.38 %) points to a prominent presence of aerosols interacting with visible light. Strong correlations with B₇ (4.73 %) indicate the effect of thermal characteristics on AOD at this level. Wd (1.18 %), Ws (1.1 %), T (1.09 %), and P (0.99 %) correlations remain at similar levels to those of the previous layer, suggesting continued influence of transport and stability.

In the 5–10 km layer of the target region, where free tropospheric dynamics dominate aerosol behavior, geographic correlations with long (22.92 %) and lat (21.46 %) show a slight increase compared to the 3–5 km layer, potentially indicating the influence of large-scale transport patterns in the free troposphere. Day (8.59 %) and month (4.91 %) correlations further weaken, reinforcing the diminishing impact of surface-driven cycles. B₁ (5.25) still shows the highest correlation, indicating aerosols that interact with visible light and are likely transported from lower altitudes or distant sources. B₇ (4.6 %) also maintains a strong correlation, indicating a potential sensitivity to upper-level atmospheric dynamics. Notably, P (1.41 %) now shows the strongest correlation among meteorological variables, surpassing Wd (1.13 %), T (1.1 %), and Wd (0.84 %). This emphasizes the critical role of atmospheric stability in controlling AOD distribution in the free troposphere. Despite meteorological data showing low statistical importance (below 2 %) in estimating AOD across all altitude layers, they have significant physical relevance in processes unique to each layer. LC (0.42 %) shows the weakest correlation across all layers, further suggesting its limited direct influence on AOD at this altitude.

4.1.2 Meteorological data and land cover feature selection

As previously noted, LC, T, P, Ws, and Wd are key features in AOD estimation. To further understand their impact on machine learning model performance in SEVIRI multi-layer AOD values, we conducted 16 cases of experiments with varied meteorological and LC feature combinations, validated using CALIOP AOD retrievals. Our analysis, depicted in Fig. 6 and Tables S1–S2 in the Supplement, is summarized using statistical metrics like R², R, RMSE, and MAE.

Figure 6 The impacts of input features on the retrieval of SEVIRI multi-layer AOD values, represented as R² metrics, for the RF and XGB models. Each row displays the results for the annual period and the four seasons (winter, spring, summer, and autumn). The four colors in each bar plot indicate the R² values for AOD at 1.5, 3, 5, and 10 km layers.

Our findings indicate that, for most cases across annual and seasonal datasets, combining these features with B_i has a negligible impact on the 1.5 km layer. In contrast, integrating Ws and Wd with SEVIRI bands (case 5) consistently yields the highest R² and lowest RMSE values across both annual and seasonal datasets. Incorporating these features significantly enhances R² values across models, with substantial increases ranging from 0.75 to 0.99 and 0.89 to 0.99 observed in R² and decreases ranging from 0.075 to 0.009 and 0.033 to 0.002 in RMSE, for both XGB and RF models from case 1 to case 5 in the 10 km layer. These statistical values underscore the crucial role of Ws and Wd in influencing the spatial and temporal properties of atmospheric aerosols, particularly in the 10 km layer. Physically, Ws and Wd are known to be primary drivers of aerosol transport across various spatial scales. High Ws can lift dust and other particulate matter into the atmosphere, while Wd affects the regional and long-range advection of aerosols. This mechanism is particularly impactful in the 5–10 km altitude range, where aerosols experience less drags and can be transported over long distances with minimal settling, especially in dry conditions (Pérez et al., 2006a, b; Georgoulias et al., 2016; Nicolae et al., 2019).

However, integrating T and P features, as seen in cases 3, 4, and 9, notably enhances AOD accuracy at 3, 5, and 10 km altitudes. This improvement is attributed to P reflecting changes in aerosol vertical layers, influencing aerosol diffusion capacity, while T is closely linked to atmospheric aerosol distribution by altering air movement dynamics. T affects atmospheric stability by controlling the thermal stratification of air masses; higher T can destabilize the atmosphere, promoting vertical mixing and lifting of aerosols. However, the impact of T on aerosol concentration diminishes with altitude, especially in stable layers where thermal inversion limits upward transport (Pérez et al., 2006a, b; Choobari et al., 2014; Xu et al., 2024). P, on the other hand, can reflect shifts in the PBL height, which influences the vertical distribution of aerosols. A higher PBL allows more aerosols to disperse vertically, enhancing their presence in layers such as AOD₅, while a lower PBL restricts aerosols closer to the surface. Despite these influences, T and P are generally less effective than wind-related factors (cases 6, 7, and 8), particularly in upper layers, because wind-driven advection predominantly controls the lateral movement of aerosols, especially during seasonal changes when wind dynamics vary substantially (Nicolae et al., 2019; Georgoulias et al., 2016; Ortiz-Amezcua et al., 2017; Granados-Muñoz et al., 2016). However, T and P impact PBL dynamics, but these effects may be limited in certain seasons due to more stable atmospheric conditions, where lower T and P fluctuations are less conducive to vertical mixing. For example, case 12 including T and P may show less impact on AOD₅ and AOD₁₀ layers in winter, when aerosols tend to stay close to the surface due to limited vertical convection. In contrast, warmer seasons promote vertical convection, enhancing the influence of T and P in predicting AOD values across altitude layers. Additionally, the proposed learning models can capture complex, non-linear relationships among features, but they may not always prioritize individual variables unless they strongly affect the target variable. As a result, even if T or P is expected to influence AOD, their impact may be overshadowed by more influential Ws and Wd (cases 13, 14, 15, and 16).

Conversely, LC plays a more localized role in AOD estimation, particularly in lower layers (AOD_1.5), as it impacts the sources and types of aerosols present near the surface. However, the effect of LC diminishes with altitude due to decreased influence on vertical transport; aerosols released from surface sources are progressively diluted as they disperse upward. This observation is confirmed by our findings that cases with LC (e.g., cases 2, 8, 13, 14, and 16) did not consistently outperform those with purely meteorological features, especially in higher layers (AOD₅ and AOD₁₀). Physically, this limitation arises from the fact that the vertical distribution of aerosols across different atmospheric layers over the European troposphere is more heavily influenced by continental and regional transport patterns, atmospheric stability, and meteorological conditions than localized land cover characteristics (Zhao et al., 2019).

In conclusion, the validation of our models using CALIOP AOD retrievals highlights the robustness of our findings, particularly the critical role of Ws and Wd in enhancing AOD estimation accuracy. The consistency of these findings across the RF and XGB models, evaluated at different temporal scales (annual and seasonal), highlights the critical role of Ws and Wd in AOD estimation at both the 5 and the 10 km layers. Consequently, we prioritize Ws and Wd, along with B_i, as the preferred input features for our models due to their demonstrated impact on improving AOD estimation accuracy.

4.1.3 Validation of seasonal modeling

Considering the substantial seasonal variations in aerosol distribution and meteorological conditions, we sought to evaluate whether adapting our proposed modeling approach to specific seasons could improve the accuracy of multi-layer AOD retrievals. Following the methodology outlined in Sect. 3, we partitioned the sample dataset, derived from 2017 to 2019 data, into four segments based on seasonal distinctions: winter (January, February, and March), spring (April, May, and June), summer (July, August, and September), and autumn (October, November, and December), as detailed in Table 3. Subsequently, we trained individual machine learning models on these seasonal and annual datasets spanning 2017 to 2018. For this analysis, we separately estimated SEVIRI multi-layer AOD values for the year 2019 using both the seasonal and the annual models. Detailed seasonal validation findings, including R² and RMSE metrics, are delineated as highlighted values in Tables S1 and S2.

Table 3 Number of samples used to train machine learning models in this study.

The XGB model exhibited acceptable performance across different seasons, with R² (RMSE) values for the 1.5 km layer as follows: 0.901 (0.0103) for spring, 0.889 (0.0347) for summer, 0.966 (0.0265) for autumn, and 0.881 (0.0392) for winter. In comparison, the RF model demonstrated improvement, boasting R² values of 0.97, 0.927, 0.937, and 0.999, with corresponding RMSE values of 0.0211, 0.0278, 0.052, and 0.0023 for spring, summer, autumn, and winter, respectively. Similarly, both the XGB and the RF models demonstrated satisfactory performance across other layers, with R² ranging from 0.81 to 0.95, 0.80 to 0.98, and 0.82 to 0.91 for the 3, 5, and 10 km layers, respectively. Performances of models generally tend to decrease in the upper layers compared to the lower layers due to the prevalent types and sizes of existing aerosols throughout most of the year, with aerosol distribution in Europe predominantly concentrated within the 1.5 and 3 km atmospheric layers. Consequently, R² and R metrics demonstrate higher values in these layers compared to the 5 and 10 km layers. Conversely, RMSE and MAE metrics are elevated at the 1.5 and 3 km layers but lower at the 5 and 10 km layers. This pattern arises from the typically higher aerosol concentrations occurring in the lower atmospheric layers, juxtaposed with lower AOD values observed in the 5 and 10 km layers.

However, the XGB annual model exhibited strong performance, achieving RMSE values of 0.0091, 0.0134, 0.0066, and 0.0059 and R² values of 0.993, 0.974, 0.985, and 0.981. Similarly, the RF annual model produced notable results, with R² values of 0.98, 0.962, 0.939, and 0.968, along with RMSE values of 0.015, 0.010, 0.0112, and 0.0066. In conclusion, annual models demonstrate consistently high predictive accuracy, with XGB showing slightly stronger performance in terms of both R² and RMSE. The minimal RMSE and high R² values indicate that both models effectively capture the multi-layer AOD values across the altitude layers on an annual basis. Compared to seasonal models, the annual models offer greater stability and reduced variability, as they are less affected by seasonal meteorological changes that can alter aerosol distribution, such as summer convection and winter temperature inversions. This consistency in annual performance provides a more reliable basis for long-term AOD retrieval, making annual models preferable for retrieving multi-layer AOD values for SEVIRI data. Therefore, we regarded the annual models as the desired models for retrieving multi-layer AOD values of SEVIRI.

4.1.4 Comparison of the models

Figure 7 presents scatterplots illustrating multi-layer AOD values estimated using the proposed annual XGB (Fig. 7a–d) and RF (Fig. 7e–h) models compared with CALIOP-retrieved AOD profiles over the European troposphere in 2019. Each subplot includes the number of points and mentioned metrics, i.e., R², R, RMSE, MAE, bias, and linear regression equations, to facilitate clear and thorough analysis.

Figure 7 Scatterplots comparing the estimated SEVIRI multi-layer AOD values derived from the proposed machine learning models with the CALIOP AOD profiles for the year 2019. The red line represents the linear fit between the two datasets.

Both models exhibit a strong correlation between the estimated values and retrievals. However, the XGB model demonstrates slightly superior performance, with R² (R) values of 0.993, 0.974, 0.985, and 0.981 (0.997, 0.989, 0.993, and 0.991) for the 1.5, 3, 5, and 10 km layers, respectively. In comparison, the RF model shows lower R² (R) values of 0.980, 0.962, 0.939, and 0.968 (0.993, 0.984, 0.972, and 0.985) for the same layers. The XGB model consistently outperforms in estimating AOD across all layers. However, the RF model demonstrates relatively higher accuracy in the 1.5 km layer compared to other layers, highlighting a significant performance gap when compared to XGB. Overall, the minimal variation in R², R, RMSE, and MAE across the models suggests comparable estimation capabilities. However, a detailed analysis reveals the superior accuracy of the XGB model in capturing AOD values, as evidenced by the slope values in Fig. 7a–d. In contrast, the slope values in Fig. 7e–h indicate that the RF model tends to slightly underestimate AOD values. In summary, while both models demonstrate proficiency, the XGB model outperforms the RF model, particularly in its accuracy for higher-altitude layers, thereby providing more reliable AOD estimations.

4.2 Validation of estimated AOD with ground lidar retrievals

To further validate our top-performing models, annual XGB and RF, we conducted an extensive analysis using data from eight EARLINET stations across Europe in 2020. We analyzed pixels within a 3 km radius of each station to compare and validate SEVIRI multi-layer AOD estimates against EARLINET values. Figure 8 presents the comparison using linear regression and validation metrics, with scatterplots for AOD at four layers. The XGB model demonstrates stronger agreement with EARLINET profiles, achieving R² values of 0.86, 0.79, 0.75, and 0.59 and RMSE values of 0.02, 0.01, 0.01, and 0.005. In contrast, the RF model shows weaker correlations, with R² values of 0.83, 0.26, 0.52, and 0.16 and RMSE values of 0.024, 0.022, 0.021, and 0.007. This can be attributed to XGB's ensemble nature and its ability to reduce bias through boosting, enabling it to handle complex and diverse datasets more effectively (Ahmed et al., 2023).

Figure 8 Scatterplots comparing the estimated SEVIRI multi-layer AOD values derived from the proposed machine learning models with the EARLINET AOD profiles across eight specified stations in Europe for the year 2020. The red line represents the linear fit between the two datasets. Note that the scales of the subplots vary due to the different ranges of AOD values at the various vertical layers (1.5, 3, 5, and 10 km).

When comparing the R² metrics of XGB AODs across different layers, it was found that XGB AODs exhibited lower R² values with EARLINET at the 10 km layer but showed significant improvement at the 1.5, 3, and 5 km layers, with R² values of 0.86, 0.79, and 0.75, respectively. This indicates a stronger correlation between XGB AOD estimations and EARLINET retrievals in these layers compared to the 10 km layer, which had an R² value of 0.593. This trend is consistent with other evaluation metrics. Closely scrutinizing Fig. 8, it becomes apparent that there are specific points revealing notable discrepancies between EARLINET and XGB AOD values. To determine the root cause of these outliers, the data were color-coded based on AOD values, revealing that the majority of outliers occurred when EARLINET retrieved low AOD values in each layer. At these points, the XGB model tends to overestimate. This tendency contributed to a low R² value (0.593) in the linear regression for the 10 km layer, as this layer contains small AOD values (0–0.05). Furthermore, the slope of the regression line for the XGB model ranges from 0.82 to 0.99, indicating that for every unit increase in the estimated values, the corresponding EARLINET AOD values increase by slightly less than one unit. This suggests a tendency for the XGB model to overestimate EARLINET AOD across all layers.

The statistical analysis of EARLINET AOD profiles at the eight specified stations, alongside estimated AOD values, is comprehensively presented in Table 4. The number of analyzed pairs (N) varies across stations, ranging from 12 at HPB to 387 at ATZ, providing robust validation of estimated AOD values with EARLINET AOD profiles. Metrics such as RMSE, MAE, and bias offer valuable insights into model performance at each station. A detailed analysis of the data reveals that the values of these metrics across the four layers are not consistently identical.

Table 4 Station-based statistical analysis of XGB-estimated SEVIRI vs. retrieved EARLINET AOD values. All metric values are scaled by $\times {\mathrm{10}}^{-\mathrm{2}}$.

Performance varies significantly across stations and layers, with notable discrepancies observed at ATZ (Greece), particularly at the 1.5 km layer, which exhibits the highest RMSE of $\mathrm{3.1}\times {\mathrm{10}}^{-\mathrm{2}}$. The substantial bias at this station indicates that the model tends to consistently overestimate the AOD at this altitude. Conversely, performance improves at 3 and 5 km, with RMSE values of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{2}}$ and $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{2}}$, respectively. At this station, the model demonstrates the best performance at the 10 km layer, with an RMSE of $\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{2}}$, compared to the other layers. These variations are likely due to frequent forest fires in Greece, as most smoke from these fires remains in the lower layers of the atmosphere (Nicolae et al., 2019). In contrast, the XGB model generally performs well at the SAL, HPB, and LLE stations, where both RMSE and bias are minimal. The RMSE at the 1.5 km layer, ranging from $\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{2}}$ to $\mathrm{2.2}\times {\mathrm{10}}^{-\mathrm{2}}$ at the IPR, WAW, INO, and THE stations, alongside the low RMSE and bias values across other layers, demonstrates good overall model performance at these stations. A closer examination reveals that RMSE and bias metrics are often elevated at the 1.5 and 3 km layers but lower at the 5 and 10 km layers. This pattern arises from the typically higher aerosol concentrations in the lower atmospheric layers compared with lower AOD values retrieved in the 5 and 10 km layers.

The discrepancies between the estimated and retrieved values could stem from the different measurement techniques employed by satellite and ground-based systems. EARLINET utilizes ground-based lidar systems to capture backscattered light from aerosols within the atmosphere by looking upward, whereas satellite measurements are performed from above, looking down. In this configuration, the lower atmospheric layers attenuate the lidar signal, resulting in reduced power to penetrate the upper layers. This attenuation can complicate the detection of aerosols in the upper layers (Grigas et al., 2015; Nicolae et al., 2019). Furthermore, these limitations may be attributed to the constraints associated with the utilization of CALIOP AODs, particularly their reduced precision in low-aerosol-concentration scenarios. This reduced precision arises from the low signal-to-noise ratio under clean weather conditions, which is often insufficient to accurately detect weak aerosol layers on the aerosol extinction vertical profile. Because both transmitted and scattered light must traverse this portion of the atmosphere, highly diffuse and/or tenuous scattering aerosol layers below the CALIOP detection threshold are ignored in CALIOP's estimates of column AOD. Consequently, weak aerosol layers that are not detected would not be retrieved, leading to decreased retrieved AODs under clean weather conditions (B. Liu et al., 2018; C. Liu et al., 2018).

Finally, the efficacy of the XGB model is clearly demonstrated by its ability to reliably estimate multi-layer AOD values compared to EARLINET-retrieved AOD profiles across various European regions.

4.3 Qualitative validation

For qualitative validation of the estimated multi-layer AOD values, two distinct approaches were utilized. The first approach involved a visual comparison of the spatial trends of the estimated AOD values with CALIOP AOD retrievals for 4 specific days across various seasons in 2019: 3 March (11:57 UTC), 30 April (12:42 UTC), 13 June (10:57 UTC), and 31 October (12:27 UTC). The second approach concentrated on two noteworthy aerosol events: a significant dust intrusion that occurred over Europe from 13 to 18 March 2022 and a volcanic eruption at Mount Etna in Italy on 14 August 2023. For these events, the estimation model generated multi-layer AOD values at 15 min intervals. This high temporal resolution enabled detailed analysis of aerosol behavior within each layer during these events.

Conducting the qualitative validation for an entire scene within the first approach is challenging due to the spatial and temporal resolution constraints of CALIOP. To address this limitation, SEVIRI scene pixels corresponding to CALIOP overpasses with temporal differences of less than 4 min were compared on the specified days. The results, illustrated in Fig. 9, indicate that the spatial trends of the estimations generally align well with the trends of CALIOP AOD retrievals in regions with both high and low AOD values across the four seasons. This alignment highlights the model's ability to provide reliable AOD estimates with enhanced temporal resolution, effectively complementing CALIOP AOD retrievals.

Figure 9 Spatial distributions and trends of SEVIRI-retrieved multi-layer AOD values compared to CALIOP retrievals for 4 specific days on 3 March (11:57 UTC), 30 April (12:42 UTC), 13 June (10:57 UTC), and 31 October (12:27 UTC) across various seasons in 2019.

The subsequent approach for qualitative validation of the model on a regional scale involved analyzing two aerosol events. The first event examined was a substantial Saharan dust plume that traversed western and central Europe between 13 and 18 March 2022. SEVIRI scenes captured during this event, illustrated in Fig. 10, were visualized at specific hours to assess how well the aerosol vertical dispersion behavior aligned with the dynamic characteristics of the event.

Figure 10 Spatial distributions of SEVIRI-estimated multi-layer AOD values during the Saharan dust event from 13 to 18 March 2022, at selected times.

On 13 March, at 13:42 UTC, the event is characterized by concentrated aerosol presence at 1.5 km over southern Europe, particularly Spain, Italy, and the Mediterranean, with progressive dispersion at higher altitudes reaching central Europe. By 14 March at 14:42 UTC, the dust plume expands further into central and northern Europe at 1.5 km, while mid-altitudes (AOD₃ and AOD₅) show notable aerosol presence over the Iberian Peninsula and central Europe, reflecting the dust's horizontal transport. At 10 km, concentrations remain low, indicating limited vertical penetration. Comparing these 2 d reveals a clear intensification and northward progression of the plume. On 15 March, at 10:12 UTC, AOD_1.5 shows sustained high concentrations over southern Europe, with increased impact in southeastern Europe, including the Balkans. At mid-levels, the plume spreads further into central and eastern Europe, reaching countries like Germany and Poland, while at 10 km, faint concentrations persist over parts of southern and central Europe. By 16 March at 09:57 UTC, the plume shows signs of dissipation, with reduced intensity at 1.5 km across Spain, Greece, and the Mediterranean and weaker signals at mid-altitudes over central Europe. On 17 March at 08:27 UTC, the dust plume is confined mainly to southern Europe at lower altitudes, while mid-altitudes exhibit limited spread and intensity, and 10 km shows negligible influence.

Finally, the spatial distributions on 18 March at 14:42 UTC, as shown in the figure, reveal a significant reduction in aerosol concentrations at all altitudes. At 1.5 km, the AOD values remain notable over southern Spain and the western Mediterranean but are much weaker compared to earlier days. At 3 km, the dust plume is localized over the Iberian Peninsula, with faint traces extending toward southeastern Europe. Higher altitudes, represented by AOD₅ and AOD₁₀, exhibit very low concentrations, indicating minimal vertical transport of the dust on the final day of the event. This vertical behavior, compared to earlier days, reflects the dissipation phase of the Saharan dust plume, as atmospheric processes like mixing, dilution, and deposition progressively weaken its intensity. Additionally, during the Saharan dust event, sub-hourly estimated multi-layer AOD values for the European continent (from 07:12 to 16:12 UTC daily on 18 March 2022) were compiled into an animation that includes 37 SEVIRI scenes over Europe, as shown in Pashayi et al. (2025a) (https://doi.org/10.5446/69730).

The other event in this approach is the Mount Etna eruption, located on the eastern coast of Sicily, Italy, one of the most active volcanoes in the world and a prominent feature in the Mediterranean region. On 14 August 2023, Mount Etna erupted in a significant volcanic event that released vast quantities of ash and aerosols into the atmosphere. This eruption produced an ash plume reaching up to 8200 m above the crater and spreading southward over the Mediterranean. The event provided another opportunity to visualize the aerosol vertical dispersion behavior aligned with the dynamic characteristics of the event. The AOD spatial distribution maps in Fig. 11 revealed distinct layers of aerosol, with a clear upward transport of particles, especially in the lower- to mid-level layers (1.5 to 3 km) at 07:27 UTC, suggesting strong vertical convective activity driven by the eruption's intensity. As time progressed, the AOD values in the upper layers increased, particularly at the 5 km altitude around 08:27 UTC, signalling a greater vertical transport of aerosols. This upward transport continued throughout the day, with the highest values observed at 10 km by 09:57 UTC, suggesting the plume had reached greater altitudes. By 13:42 UTC, the AOD values in the 10 km layer were at their peak, reflecting the maximum extent of vertical transport before the aerosols began to disperse more horizontally. After this point, the plume's vertical extent started to decrease (15:57 UTC), likely due to the long-distance transport of the aerosols across the Mediterranean region. This multi-layer analysis highlights the dynamic behavior of volcanic aerosols, providing valuable insight into their dispersal patterns. An animation of sub-hourly SEVIRI multi-layer AOD estimations during the eruption (from 07:27 to 15:57 UTC) was also generated, represented as a video in Pashayi et al. (2025b) (https://doi.org/10.5446/69731), providing a detailed view of the plume's evolution over time.

Figure 11 Spatial distributions of estimated multi-layer AOD values during the volcanic eruption on 14 August 2023, at selected times.

5 Conclusion

This study proposes a model that integrates satellite TOA reflectance data from the SEVIRI sensor, meteorological data, and land cover to estimate AOD across distinct altitude layers at 1.5, 3, 5, and 10 km. Using CALIOP AOD profiles as reference data, RF and XGB models were trained on a dataset spanning 2017 to 2018. The trained models were then used to estimate SEVIRI-based multi-layer AOD values over Europe for 2019 and 2020. The 2019 estimates were compared with CALIOP AOD retrievals, while the 2020 estimates were evaluated against EARLINET AOD retrievals, yielding the following insights.

Both RF and XGB models demonstrate high accuracy in estimating sub-hourly SEVIRI multi-layer AOD values at approximately 15 min intervals. The XGB model shows slightly superior performance, achieving R² values ranging from 0.97 to 0.99 across different layers when compared to CALIOP and from 0.59 to 0.87 when evaluated against EARLINET retrievals. Incorporating meteorological data such as T, P, Ws, and Wd, along with LC data, during model training significantly enhances the performance of the proposed frameworks. These additional features, often excluded in traditional physical AOD retrieval methods that rely solely on atmospheric radiative transfer models, greatly improve the accuracy of SEVIRI multi-layer AOD estimations. Among the meteorological variables, Ws and Wd are the most influential, resulting in higher R² values and lower RMSE across all estimated layers.

A qualitative validation was conducted by comparing the spatial trends of the estimated AOD values with CALIOP AOD retrievals. The analysis focused on SEVIRI pixels corresponding to CALIOP overpasses on 3 March, 30 April, 13 June, and 31 October 2019, with temporal differences of less than 4 min. The results demonstrate strong agreement between SEVIRI estimates and CALIOP retrievals across varying AOD levels, highlighting the model's capability to provide reliable high-resolution AOD estimates that effectively complement CALIOP data. Additionally, the study successfully estimated multi-layer AOD at 15 min intervals for two real events – a Saharan dust plume that swept across western and central Europe between 13 and 18 March 2022 and the Mount Etna eruption on 14 August 2023. The results are consistent with the physical characteristics of these phenomena, such as Saharan dust long-range transport in the upper layers of the atmosphere and a gradual increase in AOD values over time from the lower to higher tropospheric layers during volcanic events. This approach enables detailed monitoring of aerosol behavior across vertical layers of the troposphere, providing valuable insights into the dynamics of such events. In conclusion, the XGB model can estimate detailed sub-hourly 3 × 3 km² multi-layer AOD values, providing valuable insights into aerosol properties.

Our research, centered on the troposphere over Europe and validated with ground- and satellite-lidar-based AOD retrievals, provides a foundation for future studies to develop a more comprehensive approach for multi-layer AOD retrievals by incorporating an ensemble of geostationary meteorological satellites. Moreover, the current framework utilizes a limited range of input features, omitting important variables such as precipitation, NDVI, and land use, which significantly affect AOD dynamics. Future efforts will focus on improving model accuracy by including these additional factors and exploring the potential applications of the model's outputs in areas like aerosol transport analysis, air quality assessment, and climate studies to enhance its practical relevance.