Consistency of Aerosol Optical Properties between MODIS Satellite Retrievals and AERONET over a 14-Year Period in Central–East Europe

Abstract

Aerosols influence Earth’s climate by interacting with radiation and clouds. Remote sensing techniques aim to enhance our understanding of aerosol forcing using ground-based and satellite retrievals. Despite technological advancements, challenges persist in reducing uncertainties in satellite remote sensing. Our study examines retrieval biases in MODIS sensors on Terra and Aqua satellites compared to AERONET ground-based measurements. We assess their performance and the correlation with the AERONET aerosol optical depth (AOD) using 14 years of data (2010–2023) from 29 AERONET stations across 10 Central–East European countries. The results indicate discrepancies between MODIS Terra and Aqua retrievals: Terra overestimates the AOD at 16 AERONET stations, while Aqua underestimates the AOD at 21 stations. The examination of temporal biases in the AOD using the calculated estimated error (ER) between AERONET and MODIS retrievals reveals a notable seasonality in coincident retrievals. Both sensors show higher positive AOD biases against AERONET in spring and summer compared to fall and winter, with few ER values for Aqua indicating poor agreement with AERONET. Seasonal variations in correlation strength were noted, with significant improvements from winter to summer (from R² of 0.58 in winter to R² of 0.76 in summer for MODIS Terra and from R² of 0.53 in winter to R² of 0.74 in summer for MODIS Aqua). Over the fourteen-year period, monthly mean aerosol AOD trends indicate a decrease of −0.00027 from AERONET retrievals and negative monthly mean trends of the AOD from collocated MODIS Terra and Aqua retrievals of −0.00023 and −0.00025, respectively. An aerosol classification analysis showed that mixed aerosols comprised over 30% of the total aerosol composition, while polluted aerosols accounted for more than 22%, and continental aerosols contributed between 22% and 24%. The remaining 20% consists of biomass-burning, dust, and marine aerosols. Based on the aerosol classification method, we computed the bias between the AERONET AE and MODIS AE, which showed higher AE values for AERONET retrievals for a mixture of aerosols and biomass burning, while for marine aerosols, the MODIS AE was larger and for dust the results were inconclusive.

1. Introduction

Aerosols, consisting of a heterogeneous mixture of solid and liquid particles suspended in the atmosphere, exert significant influence on the Earth’s climate system. Their intricate interactions with clouds [1,2,3,4] and radiation [5,6] represent a complex component of the climate system, resulting in profound impacts on atmospheric dynamics, energy distribution, and climate variability [6,7,8]. They may significantly affect hydrological cycles, air quality, ecosystems, visibility, and precipitation and exert a great influence over human health and economic activity [9,10]. Consequently, it is very important to trace the aerosol distribution and variability in time and space to evaluate the aerosol radiative forcing affecting the climate [11,12]. Disentangling the complex interrelationships among aerosols, clouds, radiation, and climate has emerged as a paramount research focus in contemporary atmospheric science. These require a multidisciplinary approach, integrating observations from ground-based, airborne, and satellite platforms with advanced numerical modelling techniques. Observational studies, such as ground-based measurements using sun photometer and lidar networks, e.g., AERONET—Aerosol Robotic Network [13], EARLINET—the European Aerosol Research Lidar Network [14], or MPLNET—the NASA Micro-Pulse Lidar Network [15], provide crucial insights into the spatial and temporal distribution of aerosols, clouds, and associated radiative properties across different atmospheric settings. The ground-based techniques allow for the continuous acquisition of data over long periods of time (years and decades) but only at specific observational locations. These sites are often used as ground-based references for satellite data validation [16,17]. Spaceborne remote sensing with large spatial coverage enables us to derive the aerosol distribution over the globe, especially in remote regions, where ground-based observations are scarce or unavailable. A limitation of the remote sensing datasets is their rather sparse temporal and/or spatial retrieval intervals, particularly during the night (except for lidars) and often constrained by the presence of clouds. There are several datasets of global aerosol properties over land and oceans available from various satellite sensors, e.g., MODIS [18], CALIOP [19], MISR [20], POLDER [21], and OMI [22], that cover rather long timespans. The main source of errors when retrieving aerosol properties from satellite data is attributed to the extreme complexity and variability of Earth’s atmosphere–surface system. Advancements in remote sensing technologies and instrumentation have facilitated unprecedented capabilities for characterising aerosol and cloud properties with enhanced spatial and temporal resolutions even in complex atmospheric settings. These observational datasets, coupled with sophisticated retrieval algorithms and data assimilation techniques, offer invaluable opportunities for refining our understanding of aerosol optical and microphysical properties, aerosol–cloud–radiation interactions, and the aerosol impact on climate.

Significant reductions in near-surface aerosol concentrations and aerosol optical depth (AOD) have been observed in Europe during the last few decades from long-term station measurements and satellite retrievals [23,24,25]. The decrease in aerosols over Europe was mainly attributed to continuous reductions in local European anthropogenic emissions of aerosols and precursor gases since the 1980s [26], which are a result of policies for improving air quality. The decrease in aerosol concentrations is considered a cause of the increase in surface solar radiation over Europe since the 1980s [27] as well as a contributor to the Eastern European warming [28], arctic amplification [29], and increased atmospheric visibility over Europe [9] during the past three decades. Eastern Europe, which includes countries with different levels of economic and social development, is one of the least studied regions of the world in terms of atmospheric composition. Eastern Europe covers various kinds of landscapes, from plains and mountains to coastal and interior regions, leading to a great variation in aerosol characteristics over different regions. It also experiences a steady growth of industrialization and urbanization, resulting in anthropogenic activity linked to various sources of aerosol formation. A few studies have investigated certain areas within the region. For example, Sitnov et al. (2013) conducted smoke aerosol monitoring over the European part of Russia in the period of massive forest and peat fires [30]. Barnaba et al. (2011) discussed the impact of fire-related aerosols on the yearly AOD cycle over Europe [31]. Ettehadi et al. (2022) compared the satellite and AERONET AOD over southeast Europe [32]. Logothetis et al. (2020) used AERONET data to classify the aerosol types in several regions, including Central and Eastern Europe [33]. Nicolae et al. (2019) used collocated AERONET and EARLINET data to characterise the types of long-range transported aerosols over Europe [34]. Aerosol climatology and trends were investigated by Carstea et al. (2019) using AERONET observations over continental southeast Europe [35]. Studies investigating the atmosphere over large cities were carried out: Evgenieva et al. (2009) studied atmospheric aerosol optical characteristics over Sofia [36]; Stefanie et al. (2023) described the temporal variation of aerosol optical properties over Cluj-Napoca [37]; Posyniak et al. (2016) examined the long-term variability of aerosol optical thickness over Belsk [38]; Elansky et al. (2018) studied air quality and pollutant emissions [39], while Chubarova et al. (2016) investigated the temporal variability of the AOD in the Moscow region [40]; and Rupakheti et al. (2023) studied aerosol optical properties over Chisinau, Moldova, over two decades [41]. Aerosol–climate interactions were discussed by Markowicz et al. (2021), while Filonchyk et al. (2021) showed the impact of the COVID-19 lockdown on air quality over Poland [42,43]. Nevertheless, studies of temporal and spatial aerosol distribution in the whole region are very limited [44,45,46,47], demonstrating the need of further integrative studies.

This study aims to explore the spatiotemporal dynamics of aerosols across ten Central–East European countries, namely Austria, Belarus, Bulgaria, Estonia, Lithuania, Moldova, Poland, Romania, Slovakia, and Ukraine, over a 14-year period from 2010 to 2023. It examines their optical characteristics and seasonal variations, utilizing both ground-based methodologies and satellite retrievals. The dataset encompasses retrievals of aerosol properties, including the aerosol optical depth (AOD) and the Ångström exponent (AE), sourced from the moderate resolution imaging spectroradiometer (MODIS) and the aerosol robotic network (AERONET). The analysis in Section 3 initially evaluates the variability of the aerosol AOD derived from the MODIS and AERONET across four seasons. Additionally, the study presents a seasonal comparison and statistical correlation between the MODIS and AERONET in terms of AOD and AE measurements. Furthermore, the research conducts an analysis of AOD biases between the MODIS and AERONET, along with an examination of AOD trends. Finally, Section 4 explores the classification of aerosols, computed using AERONET observations of the AOD and AE considering a location-specific analysis (urban vs rural). Conclusions are presented in Section 5.

2. Methodology and Data Selection

2.1. MODIS Aerosol Data

The MODIS sensors onboard the Terra and Aqua satellites have been pivotal in acquiring columnar aerosol properties since 2000 and 2002, respectively. As constituents of the A-Train constellation, the Aqua MODIS satellite passes the equator at 13:30 local time, while Terra does so at 10:30 local time. Equipped with 36 spectral bands featuring varying spatial resolutions (250, 500, and 1000 m) and a wide swath of 2330 km in the ±55° scanning mode, the MODIS instruments ensure nearly daily global coverage by each satellite (Terra/Aqua) [48]. To retrieve aerosol optical depth, the MODIS employs three operational algorithms: Dark Target (DT), Deep Blue (DB), and Multi-Angle Implementation of Atmospheric Correction (MAIAC), which are regularly updated to enhance their accuracy [48,49,50]. The obtained AOD data are invaluable for validating, evaluating, or constraining aerosol models, conducting dynamic air pollution analyses, and monitoring air quality. The study was conducted over the Central–East European region, characterised by diverse surface types, which require the use of a combined aerosol product. The MODIS-Aqua Dark Target Land and Ocean Collection 6.1 Level 2 aerosol operational product with a spatial resolution of 10 by 10 km was selected for this study for both the Terra MODIS (MOD04_L2) [51] and Aqua MODIS (MYD04_L2) [52]. The products used are ‘Optical_Depth_Land_And_Ocean’ for the AOD at 550 nm and ‘Deep_Blue_Angstrom_Exponent_Land’ for the Ångström exponent at 440/675 nm. This product integrates data from two independent algorithms to retrieve aerosol parameters over both land and ocean surfaces. The reported uncertainty range of the DT AOD over land is ±0.05 + 15% of the AOD, whereas over the ocean, it is ±0.03 + 5% of the AOD [48]. These datasets are accessible via NASA’s Distributed Active Archive Centers (DAACs) at https://www.earthdata.nasa.gov, accessed on 1 December 2023.

2.2. AERONET Data

The AERONET initiative, a collaborative effort between the National Aeronautics and Space Administration (NASA) and LOA-PHOTONS, utilises the CIMEL Electronique CE-318 multiband sun photometer for systematic observations of global spectral aerosol optical depth [13,53]. These observational data, accessible through the NASA Goddard Space Flight Center (GSFC) website (https://aeronet.gsfc.nasa.gov/, accessed on 1 December 2023), are presented in three algorithmic iterations: Versions 1.0, 2.0, and 3.0 (V3) [54]. Each version introduces distinct quality levels—Levels 1.0 (unscreened), 1.5 (cloud screened), and 2.0 (cloud screened and quality assured). Spectral AOD, measured within the 0.34–1.64 μm range, exhibits low uncertainty within ±0.01–0.02 with a high temporal resolution (every 3–15 min) [55]. AERONET serves as a benchmark for satellite inversion correction, validation, and aerosol characterisation.

The Version 3.0 Level 2.0 (L2) AOD data from 29 AERONET stations across 10 countries (spanning January 2010 to December 2023) were employed, providing AOD data with fully automatic near-real-time quality controls and cloud screening. In this study, we only selected sites that offer at least three years of available level 2.0 data (see

Table 1. Location of 29 AERONET stations in 10 European countries and the attributed location type, urban or rural sites, mean values of Ångström exponent at 440/675 nm, mean AOD at 550 nm, and the available number of daily measurements between 2010 and 2023.

Country	Station	Location	Latitude	Longitude	Mean AE 440/675 nm	Mean AOD 550 nm	Nr. AOD
AU	Kanzelhohe_Obs	Rural	46.677	13.901	1.32	0.087	842
AU	Vienna_BOKU	Urban	48.237	16.331	1.49	0.149	957
AU	Vienna_UNIVIE	Urban	48.221	16.355	1.48	0.165	493
BG	Galata_Platform	Rural	43.044	28.193	1.48	0.156	1976
BG	Sofia_IEBAS	Urban	42.653	23.386	1.50	0.157	462
BY	Minsk	Urban	53.920	27.601	1.47	0.152	1435
EE	Toravere	Rural	58.264	26.466	1.35	0.112	1723
LT	Irbe_lighthouse	Rural	57.750	21.722	1.41	0.117	365
MD	Moldova	Urban	47.000	28.815	1.50	0.168	1465
PL	Belsk	Rural	51.836	20.791	1.45	0.185	1292
PL	Debrzyna_PULS	Rural	53.782	16.592	1.38	0.135	207
PL	POLWET_Rzecin	Rural	52.762	16.309	1.39	0.154	181
PL	Raciborz	Rural	50.083	18.191	1.43	0.178	1032
PL	Strzyzow	Rural	49.878	21.861	1.50	0.152	1061
PL	Warsaw UW	Urban	52.210	20.982	1.39	0.167	593
RO	Bucharest_Inoe	Urban	44.348	26.028	1.52	0.216	759
RO	CLUJ_UBB	Urban	46.768	23.551	1.53	0.184	1231
RO	Eforie	Rural	44.075	28.632	1.48	0.176	1106
RO	Gloria	Rural	44.599	29.359	1.50	0.165	1719
RO	Iasi_Loasl	Urban	47.193	27.555	1.52	0.178	1408
RO	Magurele_Inoe	Urban	44.348	26.030	1.53	0.175	1643
RO	Section-7_Platform	Rural	44.545	29.446	1.42	0.144	620
RO	Timisoara	Urban	45.746	21.227	1.46	0.183	970
SK	Poprad-Ganovce	Rural	49.035	20.322	1.53	0.128	1279
UA	Kyiv	Urban	50.363	30.496	1.47	0.173	1674
UA	Kyiv_AO	Urban	50.452	30.498	1.39	0.168	247
UA	Lugansk	Urban	48.570	39.364	1.39	0.204	60
UA	Martova	Rural	49.936	36.953	1.41	0.148	66
UA	Sevastopol	Rural	44.615	33.517	1.51	0.166	659

).

For the present investigation, the AOD at 550 nm was calculated through the linear interpolation of AOD values at 440 nm and 675 nm. This computation was conducted for the purpose of collocating with the AOD retrieved by MODIS, which was also obtained at 550 nm. The Ångström exponent (AE) from AERONET was provided based on measurements at 440/675 nm, while we additionally derived the AE at 440/870 nm to facilitate the aerosol classification described in Section 4.

2.3. Collocation and Intercomparison

For this study we used daily AERONET Level 2, Version 3 aerosol data from 29 monitoring stations located in Central–East European (CEE) countries, as outlined in Table 1. These datasets encompass AOD measurements at 550 nm and the Ångström exponent at 440/675 nm spanning over a 14-year period from 2010 to 2023. Given that AERONET stations record measurements at intervals ranging from 3 to 15 min, a minimum threshold of 8 individual measurements was established as a criterion for extracting daily means. Table 1 also details the number of daily means selected from each individual station. Additionally, the MOD04_L2 and MYD04_L2 aerosol products from the Terra and Aqua MODIS instruments, part of the Collection 6.1 data release, were used to acquire AOD data at 550 nm and the Ångström exponent at 440/675 nm for the same time frame as the AERONET dataset. For this study, a spatial resolution of 10 × 10 km was adopted, as it was deemed appropriate for regional-scale analyses and maintained consistency with previous investigations carried out in the same geographical context [44,45]. Regarding spatial collocation, the study adopted a 30 km radius as a compromise between ensuring spatial representation and achieving a statistically significant number of corresponding measurements. Matchups were excluded if less than 25% of pixels fell within the selected 30 km radius. After applying the collocation criteria, the study retained 51.27% matches with Terra MODIS and 43.45% matches with Aqua MODIS from the total AERONET daily measurements available between 2010 and 2023. Supplementary Table S1 shows the matching percentages for individual AERONET stations. The lower percentage of Aqua MODIS collocations compared to Terra MODIS is attributed to increased cloud coverage during the afternoon hours as well as the Aqua MODIS overpass coinciding with low-light or night-time conditions during winter months. These factors contribute to fewer successful AOD retrievals and, consequently, fewer collocations with AERONET measurements. Spatial distribution and temporal variations were analysed based on daily values, which were aggregated on seasonal and annual bases.

3. Analysis of Aerosol Properties Derived from MODIS and AERONET

This section presents an in-depth analysis and interpretation of the comparison between AOD data derived from MODIS satellite observations and ground-based measurements obtained from AERONET stations. This comparison is essential for evaluating the accuracy and reliability of satellite-derived AOD estimates, which play a crucial role in understanding aerosol distribution and its impact on climate. By comparing these datasets, we seek to assess the performance of the MODIS in capturing aerosol characteristics and variability and to identify potential biases or discrepancies in the region. We present the statistical results for both MODIS sensors to enhance the assessment of daily representativeness at specific locations relative to each individual AERONET station.

3.1. Seasonal Analysis of AOD in Central–East Europe for the Period 2010–2023

Aerosol concentrations exhibit significant seasonal variations influenced by factors such as meteorological conditions, seasonal changes in emission sources, and atmospheric circulation patterns. Comparing the MODIS AOD with the AERONET AOD allows for the assessment of the consistency and accuracy of satellite-derived aerosol retrievals. Evaluating the agreement between these two datasets across different seasons provides valuable information about the performance of remote sensing techniques under varying atmospheric conditions. Discrepancies between MODIS and AERONET AOD measurements may also reveal biases or limitations in satellite retrieval algorithms.

3.1.1. AOD Maps from MODIS Retrievals and Measurements at AERONET Stations

The MODIS AOD maps and AERONET AOD provide useful information regarding the spatial variations in the AOD over CEE, illustrating the agreement and discrepancies between satellite and ground-based measurements during different seasons.

In Figure 1 four maps depicting the seasonal variation in the aerosol optical depth (AOD), derived from MODIS Terra (MOD04_L2) datasets, for winter, spring, summer, and autumn seasons are presented for the CEE region covering the period from 2010 to 2023. Simultaneously, mean AOD values obtained from 29 AERONET stations distributed across the region are depicted on these maps. Seasonal variations in the AOD are discernible, with the lowest values retrieved during winter and the highest during summer. Persistent increases of AOD levels throughout the seasons is observed at specific AERONET stations, such as the Bucharest station, in southern Romania. The disparity in the number of AOD measurements between winter and summer is attributed to meteorological factors, including heightened cloud cover and precipitation during winter, which impedes both AERONET photometer measurements and MODIS retrievals. Additionally, the presence of bright surfaces, particularly snow, during winter months impacts MODIS retrievals. Winter maps indicate significant data gaps (approximately 30%) primarily over mountainous regions, such as the Carpathian Mountains and the Alps, as well as in the northern and eastern parts of the studied region. Elevated AOD values are predominantly observed in Poland and Estonia during winter, mean values occasionally exceeding 0.25, while the overall average AOD across the region remains below 0.1. However, caution is warranted in interpreting these high mean values, as they may be influenced by limited measurements and subpixel-level cloud or snow contamination [56]. As spring progresses, MODIS Terra coverage improves, with only sporadic missing data observed over mountainous areas and in the northern regions. The southeastern area of Eastern Europe, notably Ukraine, exhibits the lowest AOD levels during this season, with mean values below 0.1. Conversely, elevated AOD values are recorded in Poland and the Czech Republic, exceeding 0.2. Stations in Chisinau (Moldova), Warsaw (Poland), and Bucharest (Romania) exhibit the highest AOD values, with mean values surpassing 0.25. Moreover, MODIS data indicate elevated AOD levels around 0.2 over the Black Sea, while AERONET stations report an average AOD of approximately 0.15. The discernable difference between the Black Sea area and Eastern Ukraine may be attributed to changes in surface reflectance estimations between the land and ocean DT algorithm. Locations with higher AOD are also observed over wider stretches of the Danube and Dnipro rivers, highlighting potential limitations in the retrieval algorithm’s sediment masking capability. In summer, MODIS Terra retrieves an average AOD exceeding 0.2 across Eastern Europe, with exceptions noted in areas with cleaner air, such as the Carpathian Mountains and the Alps. Lower AOD values are observed in Crimea and southern Ukraine, ranging between 0.10 and 0.15. During autumn, AOD levels decrease, reaching a maximum of 0.2 in Poland and Estonia. Conversely, lower AOD values are observed in rural areas of eastern Ukraine (approximately 0.05) as well as in Moldova, Romania, Austria, and Bulgaria, with values around 0.1. AERONET measurements generally align with MODIS observations, indicating AOD values between 0.15 and 0.2, except for stations in Bucharest and Timisoara, which exhibit higher mean values. Additionally, stations located along the Black Sea coast register mean AOD values around 0.1.

The maps in Figure 2 show the seasonal variability of the AOD derived from MODIS Aqua. Seasonal patterns in the AOD demonstrate similarities between MODIS Terra and MODIS Aqua, with the lowest AOD values typically occurring during winter and autumn and the highest during summer. However, notable differences arise, particularly during winter, where MODIS Aqua exhibits more than 50% missing data across the region. This discrepancy can be attributed to several factors. First, the presence of bright surfaces covered by snow impacts the accuracy of retrievals from both MODIS platforms. Additionally, the temporal offset between MODIS Terra and MODIS Aqua overpasses further complicates AOD retrievals during winter months. Specifically, MODIS Aqua passes over the region approximately 3–4 h later than MODIS Terra, which, given the local time in Eastern Europe, could coincide with or occur after sunset. This temporal misalignment inhibits AOD retrievals from MODIS Aqua, contributing to the observed data gaps during winter. Thus, while both MODIS sensors exhibit similar seasonal trends in AOD, the delayed overpass of MODIS Aqua and the associated challenges with low-light retrievals result in larger amounts of missing data during winter months in the region.

3.1.2. Seasonal Comparison between MODIS and AERONET AOD Retrievals

In this section, the study investigates a comparative analysis of AOD data derived from MODIS satellite observations and ground-based measurements. The statistical analysis makes use of the linear regression function to calculate the R-value, the R-squared (R²) value, and the root mean square error (RMSE) to assess the correlation between the MODIS and AERONET AOD. This examination offers a comprehensive exploration of the spatiotemporal variability in the AOD over Central–East Europe across different seasons, providing valuable insights into the agreement and disparities between satellite and ground-based measurements.

Figure 3 illustrates the comparative analysis between AOD measurements obtained from the AERONET at 550 nm and the MODIS Terra AOD at 550 nm, with collocation as detailed in Section 2.3, encompassing four distinct seasons, spanning from 2010 to 2023. The plots depict the statistical correlation between the two datasets for each season across all AERONET stations. Additionally, marginal histograms display the distribution of data points per bins of an AOD of 0.025 for both AERONET and MODIS datasets. Notably, both datasets exhibit a lognormal distribution for the AOD. The scatter plot is colour mapped using the corresponding AERONET (AE), calculated with AOD values at 440 and 870 nm. The number of data points included in the analysis varies across seasons, with the highest count observed during summer (5904 points) and the lowest during winter (888 points). In spring and autumn, the available data points amount to 3721 and 3600, respectively.

Figure 3a specifically presents the comparison during winter, where the statistical results indicate an R² value of 0.58 and an RMSE of 0.048 (p-value of 0.015). The AE values predominantly range between 1.2 and 2.1, suggesting the presence of mixed or polluted aerosols. Additionally, some smaller AE values (around 0.3) are observed at larger AOD values, potentially indicating dust events originating from Northern Africa. As spring progresses (Figure 3b), the correlation between the two datasets strengthens, with an R² value of 0.62 (p-value of 1.106 × 10⁻⁶). The RMSE of 0.069 aligns with the spread of data at larger AOD values, while a bias is noticeable at small AOD values, where AERONET AOD values exceed 0.01 and the MODIS exhibits negative AOD values up to −0.05. Negative values in the MODIS AOD product can occur due to calibration uncertainties, surface reflectance estimation errors, cloud contamination, atmospheric correction uncertainties, and low aerosol loading, all of which contribute to the inherent uncertainties in the retrieval process. These negative AOD values are typically close to zero and fall within the uncertainty range of the MODIS AOD product. Most AE values fall within the range of 0.6 to 2.1, with occasional points below 0.3 potentially linked to distinct dust events, as discussed in Ajtai et al. (2020) and references therein [57]. During summer (Figure 3c), the correlation between the datasets further improves, reflected in an R² value of 0.76 with an RMSE of 0.063. (p-value = 0). Although a bias is still discernible, the MODIS exhibits fewer points below 0 compared to spring. The AE predominantly ranges above 1.5 to 2.1 at larger AOD values, with lower AE values observed at lower AODs, alongside scattered points of low AE (<0.3) across the dataset. In autumn (Figure 3d), the correlation between the two datasets remains robust, characterised by an R² of 0.75 and an RMSE of 0.05 (p-value = 0.0976). AE values predominantly range between 0.6 to 2.1, with few exceptions noted.

Figure 4 presents the comparison between the AERONET AOD at 550 nm and collocated MODIS Aqua AOD at 550 nm across the four seasons spanning from 2010 to 2023. During winter, the total number of collocated points is 662, which is 226 points fewer than those observed for the AERONET–MODIS Terra comparison for the same time period. This discrepancy aligns with the observations made in Figure 1a and Figure 2a, where MODIS Aqua exhibited approximately 20% more missing data than MODIS Terra over the entire Eastern Europe region during winter. Similarly, for the other seasons, a lower number of collocated data points is observed compared to the AERONET–MODIS Terra analysis presented in Figure 3: 3127 in spring, 5051 in summer, and 3122 in autumn. The observed difference in the number of retrievals between MODIS Aqua and MODIS Terra can be explained by the satellite overpass time relative to the local solar time. MODIS Aqua has an equatorial crossing time of approximately 13:30 p.m., which is later in the day compared to MODIS Terra’s equatorial crossing time of around 10:30 a.m. [58]. This later overpass time for MODIS Aqua results in fewer successful AOD retrievals due to higher cloud coverage. The increase in atmospheric turbulence and cloud coverage in the afternoon hours can hinder the retrieval of the AOD and collocation with AERONET data. During spring in Eastern Europe, mornings typically start with cooler temperatures, which may lead to more stable atmospheric conditions compared to the warmer afternoons. As the day progresses and temperatures rise, the atmosphere becomes more prone to convective processes, which can lead to increased atmospheric turbulence and cloud development. During summer and autumn, atmospheric turbulence and cloud cover may exhibit similar trends to those observed in spring. Mornings in Eastern Europe during summer can start relatively cool, with temperatures rising quickly throughout the day. This rapid warming leads to pronounced convective activity in the atmosphere, resulting in increased turbulence and cloud development, particularly in the afternoon and early evening. In general, in autumn, as temperatures tend to be cooler overall compared to summer, convective activity may be somewhat less pronounced. However, there can still be significant variability, and turbulence and cloud coverage may still be higher in the afternoon under certain weather conditions [59,60,61]. Despite the seasonal variations, the AERONET typically measures the atmospheric column in the morning when cloud coverage is low. However, MODIS Aqua, passing later in the day, may encounter cloudy skies, inhibiting the retrieval of the AOD and collocation with AERONET data.

Furthermore, the statistical analysis reveals lower values for R² in the MODIS Aqua–AERONET comparison compared to the MODIS Terra–AERONET analysis. In winter, the R² is 0.53 with an RMSE of 0.055 (p-value of 4.616 × 10⁻⁶). In spring, the R² is 0.60 with an RMSE of 0.069 (p-value of 0). During summer, the R² is 0.74 with an RMSE of 0.062 (p-value of 6.986 × 10⁻⁵), while in autumn, the R² is 0.70 with an RMSE of 0.05 (p-value of 0). Throughout all seasons, the AE values predominantly range between 0.6 and 2.1, with exceptions of lower AE values (<0.3) potentially related to dust events.

The analysis reflects the fact that correlations improve with the increase in the number of observations, while the RMSE improves with a decrease in AOD values. In winter the lower R² values and low RMSE can be attributed to several factors. The lower number of observations may contribute to a decreased statistical robustness, leading to lower correlations. Also, the lower variability in AOD values during winter and autumn, with generally lower overall AOD values compared to spring and summer, can lead to a lower RMSE, as seen in this study. Since the observation count between spring and autumn is very similar, the small decrease in R² and RMSE in spring can be attributed to the increase in dust events and biomass burnings. Based on the p-values, all the correlations analysed are statistically significant.

The MODIS AOD uncertainty, expressed as EE = ±(0.05 + 0.15 × AOD) AOD units, accounts for systematic biases with the 0.05 component, while the 0.15 × AOD component represents AOD-dependent uncertainties that become more significant at higher AOD values. During seasons with generally higher AOD values, such as spring and summer, the AOD-dependent uncertainty component becomes more prominent. These higher uncertainties associated with larger AOD values can contribute to larger residuals, resulting in higher RMSE values. During seasons with generally lower AOD values, such as winter and autumn, the AOD-dependent uncertainty component becomes less significant, but the systematic biases can still impact the agreement between datasets.

Between Figure 3 and Figure 4, there is a slight deviation in slope values, indicating either a positive or negative bias of the MODIS AOD. Thus, in winter, both MODIS sensors (Figure 3a and Figure 4a) underestimate the AOD with Aqua, showing a larger deviation then Terra indicated by slope values of 0.82 and 0.88, respectively. In spring and autumn, both sensors seem to agree well with slope values close to 1. In summer, MODIS Terra slightly overestimates the AOD more than MODIS Aqua, for a slope of 1.11 vs. 1.04. Section 3.2 provides an analysis of the mean temporal bias and bias distributions at each AERONET location.

3.2. Analysis of the AOD Biases Retrieved from the MODIS at AERONET Stations

In this section, we first analysed the temporal mean bias between the MODIS (Terra and Aqua) and AERONET AOD. For this, we calculated the error ratio between both the MODIS Terra and Aqua AOD with respect to the AERONET AOD using the following equation:

ER = (MODIS AOD − AERONET AOD)/EE

(1)

where EE is the estimated error:

EE = ±(0.05 + 0.15 × AERONET AOD)

(2)

The ER compares the actual error (bias) to the EE. Values of −1  ≤  ER  ≤  1 mean that the actual errors are smaller than the EE, whereas |EE|  >  1 indicates a poor match. Figure 5, which illustrates the temporal variation of the ER for both MODIS Terra and Aqua (a) and the corresponding count of collocated daily means for these two instruments (b). We can observe a strong seasonality in the number of coincident retrievals, with a very low number of days available during winter. Correspondingly, the MODIS Aqua ER values during winter are the most negative, with a few values less than −1, indicating a poor match with AERONET measurements. On the other hand, the ER values are visibly larger for MODIS Terra but within the bounds of uncertainty (exception in March 2012 when ER = 1.02).

The temporal mean biases for the entire time series are 7.4 × 10⁻⁴ for MODIS Terra with an SD of 0.027 and −0.023 for MODIS Aqua, with an SD of 0.034. The corresponding temporal mean ER of 0.008 and −0.325 for Terra and Aqua, respectively.

Both sensors have higher positive biases against the AERONET in spring and summer than in fall and winter. The consistent higher biases observed in Terra in comparison to Aqua are consistent with the analysis outlined earlier. This bias difference between the two MODIS sensors appears to persist throughout the analysed time span. The greater number of collocations in the Terra dataset compared to Aqua’s could be due to diurnal cloud patterns, resulting in cloudier conditions during Aqua’s afternoon overpass than Terra’s morning pass. This increased cloud cover in the afternoon might limit the number of feasible collocations. Alternatively, discrepancies in instruments affecting retrievals could also contribute to these differences.

In Figure 6, we represent the biases in AOD values (ΔAOD) at 550 nm between ground-based AERONET measurements and satellite-retrieved MODIS data for 29 locations. These station-specific ΔAODs are plotted as the probability density function (PDF) for both the Terra (Figure 6a) and Aqua (Figure 6b) MODIS. Since AERONET measurements of the AOD are widely considered a reference, this analysis indicates the site-specific accuracy of the MODIS AOD product. The sequence of the PDFs, arranged from top to bottom, is determined by the mean MODIS bias displayed in progression from the most negative to the most positive values. Thus, each site-specific mean bias is indicated along the colour bar. To better represent the transition between negative and positive bias, we selected a divergent colour scheme that progresses from blue, indicating negative values, to red, signifying positive values, with a lighter shade indicating values closer to 0. The vertical black line was plotted to represent a consistent visual separation between positive and negative regions of the PDFs. Since we consider AERONET measurements as ‘ground-truth’, positive mean bias values indicate an underestimation of the MODIS AOD compared to the AERONET, while negative values signify that MODIS retrievals overestimate those of the AERONET. For a more comprehensive analysis on site-specific biases,

Table 2. Number of daily collocated MODIS Terra and MODIS Aqua AOD with AERONET AOD, and their mean biases and confidence intervals for individual stations.

Station	MODIS Terra (MOD)			MODIS Aqua (MYD)
	Nr. AOD	Mean Bias	CI 95%	Nr. AOD	Mean Bias	CI 95%
Belsk	628	0.007	−0.004, 0.020	427	0.040	0.025, 0.055
Bucharest_Inoe	459	0.037	0.023, 0.051	366	0.066	0.049, 0.082
CLUJ_UBB	748	0.019	0.008, 0.029	575	0.050	0.037, 0.061
Debrzyna_PULS	77	−0.038	−0.089, 0.013	62	0.021	−0.017, 0.058
Eforie	672	0.003	−0.007, 0.013	656	0.023	0.010, 0.035
Galata_Platform	1174	−0.028	−0.035, −0.019	1180	0.007	−0.00, 0.0153
Gloria	1053	−0.024	−0.032, −0.014	1036	−0.008	−0.017, 0.000
Iasi_Loasl	905	0.014	0.004, 0.024	734	0.053	0.041, 0.064
Irbe_lighthouse	226	−0.085	−0.106, −0.063	246	−0.035	−0.051, −0.018
Kanzelhohe_Obs	232	−0.086	−0.104, −0.068	188	−0.063	−0.081, −0.045
Kyiv	754	−0.012	−0.024, 0.001	591	0.019	0.006, 0.031
Kyiv_AO	106	0.015	−0.010, 0.040	71	0.007	−0.039, 0.052
Lugansk	21	−0.021	−0.071, 0.030	17	0.155	0.081, 0.229
Magurele_Inoe	935	0.001	−0.008, 0.011	748	0.026	0.013, 0.038
Martova	41	−0.010	−0.051, 0.031	30	−0.049	−0.107, 0.009
Minsk	516	0.014	−0.000, 0.028	395	0.014	−0.003, 0.030
Moldova	847	−0.001	−0.012, 0.010	656	0.034	0.020, 0.047
POLWET_Rzecin	79	0.011	−0.017, 0.038	67	0.010	−0.021, 0.041
Poprad-Ganovce	406	−0.051	−0.066, −0.036	313	−0.035	−0.050, −0.018
Raciborz	396	0.003	−0.013, 0.019	338	0.041	0.026, 0.055
Section-7_~Platform	362	−0.032	−0.046, −0.018	378	−0.006	−0.021, 0.008
Sevastopol	420	−0.018	−0.031, −0.004	377	−0.002	−0.017, 0.013
Sofia_IEBAS	242	−0.039	−0.059, −0.017	187	0.014	−0.006, 0.034
Strzyzow	526	0.011	−0.001, 0.023	447	0.019	0.005, 0.031
Timisoara	650	0.007	−0.005, 0.019	523	0.035	0.021, 0.049
Toravere	636	−0.041	−0.053, −0.029	527	−0.007	−0.020, 0.007
Vienna_BOKU	444	−0.026	−0.040, −0.012	381	0.008	−0.007, 0.023
Vienna_UNIVIE	268	0.006	−0.012, 0.025	235	0.007	−0.016, 0.030
Warsaw UW	290	−0.015	−0.032, 0.001	211	0.047	0.025, 0.068

additionally provides mean biases for both MODIS sensors. The biases presented in Table 2 facilitate a straightforward analysis, highlighting the sensor-specific differences more effectively.

A closer examination of Figure 6 reveals variations in biases across the AERONET sites. MODIS Terra exhibits a bias range spanning from −0.086 (Kanzelhohe_Obs) to 0.037 (Bucharest_Inoe), while MODIS Aqua shows values ranging from −0.063 (Kanzelhohe_Obs) to 0.155 (Lugansk). Between the two sensors, MODIS Terra shows mean negative biases at 16 out of 29 sites, while MODIS Aqua exhibits this behaviour in 8 locations. This suggests that MODIS Terra AOD retrievals are mostly overestimated, while Aqua retrievals are mostly underestimated compared to AERONET observations. This tendency was also observed in Section 3.1.2. Even so, the differences in site-specific biases between the two sensors is typically below 0.05. Larger deviations, such as those observed at Debrzyna_PULS, Martova, and Lugansk, can be attributed to a limited number of matching retrievals. In cases with such a low sample size, any results should be considered statistically irrelevant. The MODIS’s largest negative bias at Kanzelhohe_Obs can be attributed to site-specific factors since the station is located in the Alpine region, registering the lowest mean AOD of 0.115. Other locations with an AOD < 0.15 register a similar negative bias due to an overestimation of surface reflectance in low AOD conditions [62,63]. Any difference between the two sensors based on a location-specific criterion (urban vs. rural environment) is for the most part inconclusive. Despite these facts, the Terra MODIS consistently provides a negative mean bias in five out of six sites located in coastal or oceanic environments (Galata_Platform, Gloria, Section-7_Platform, Irbe_lighthouse, Sevastopol), while at the sixth site (Eforie), the bias is close to 0.

Regarding ΔAOD distributions, most PDFs exhibit a dominant peak situated close to 0, which indicates a good agreement between the AERONET and MODIS AOD. Locations that show a second large peak lack symmetry in the ΔAOD and are mostly skewed in one direction. This observation holds true for Debrzyna_PULS, Martova, Lugansk, and POLWET_Rzecin, which exhibit the lowest number of matching AOD pairs. Other examples with two peaks include Poprad-Ganovce and Irbe_lighthouse, both negatively skewed and both averaging low AOD values, which the MODIS tends to overestimate [62,63]. Locations that exhibit a central dominant peak and two secondary symmetrical peaks, such as Strzyzow and Vienna_BOKU, are only present in either the Terra or Aqua dataset, demonstrating some inconsistency between the two retrievals. Most secondary peaks, either skewed or symmetrically distributed, typically fall between ±0.1 and ±0.2, denoting inconsistencies in conditions of low aerosol loading, observed in Figure 4. PDFs dominated by more elongated distributions include urban sites with a higher count of matching AOD pairs, such as Magurele_Inoe, Iasi_Loasl, Moldova, Kyiv, CLUJ_UBB, and Timisoara. In this case, a more spread-out PDF can indicate multiple aerosol types, resulting in inconsistencies between the AERONET and MODIS retrievals. Tail distributions are mostly inconsistent between the two MODIS datasets and do not follow any concrete pattern. Apart from sites exhibiting a low number of matching AOD pairs, there are no discernible outliers within the distributions of the ΔAOD. The overall range of biases observed in the present study is consistent with values reported in prior investigations. Specifically, Filonchyk et al. (2020a, 2020b) documented bias estimates between ±0.08 and ±0.04 at Minsk, Kyiv, Moldova, Belsk, and Eforie when evaluating both MODIS datasets between 2002 to 2019 [44,45].

By assessing the confidence intervals (CI) in Figure 7, we can discern the significance of the MODIS AOD biases.

Of the 29 locations, 16 for Terra and 13 for Aqua have CIs overlapping zero suggesting that at these locations, the mean biases may not be statistically significant. The largest CIs are present in locations with the lowest number of observations as opposed to narrow CIs at locations with a high number of observations. In both MODIS datasets, CIs that do not overlap are present in urban vs seaside locations. This could be attributed to differences in the retrieval scheme between the “land” and “ocean” algorithms. Some locations exhibit large CIs, and we consider the biases statistically insignificant due to the low number of observations and large variability. The CI values for each location are present in Table 2.

Based on this analysis, it is difficult to discern whether the MODIS bias can be attributed as purely site-specific (meteorological and surface reflectance), sensor-specific (viewing angle and calibration dependencies), or algorithmic dependent. As such, a more in-depth analysis is needed to fully address this issue.

3.3. Trend in AOD over Central–East Europe Derived from AERONET Network and Collocated MODIS Terra and Aqua between 2010 and 2023

Figure 8a depicts the monthly mean AERONET AOD at 550 nm across the CEE region from January 2010 to June 2023 alongside its associated standard deviation. The trends in AOD depicted in this plot were derived from daily measurements, collected from 29 AERONET stations. Notably, a decreasing trend is discernible over this timespan, the trend supported by the p-value of monthly means of 2.32 × 10⁻⁵. For instance, the mean annual AOD value in 2010 stood at 0.174, which notably decreased to 0.116 by the initial six months of 2023 (see Figure S3 in the Supplementary Materials). The sharpest decrease in the AOD is observed between 2020 to 2021, coinciding with the COVID-19 lockdown across Europe. This decrease in the AOD was also identified by Filonchyk et al. (2021) over Poland within the same study period [43]. Our findings support the observed trends from prior research, showcasing a monthly mean AOD trend of −0.00027 and an annual trend of −0.0033 across the 10 countries where the AERONET measurements were collected. Other studies, such as Nicolae et al. (2019), also support a descending AOD trend between 2007 and 2017, indicating a decrease in aerosol loading across CEE [34].

Figure 8b,c present the monthly mean AOD trends computed using MODIS-collocated data for Terra and Aqua. While these two plots show a smaller decreasing trend than for the AERONET AOD, they are supported by the p-values of 1.18 × 10⁻³ and 4.32 × 10⁻⁴, respectively. For MODIS Terra, the mean AOD decreases from 0.142 to 0.1, with a monthly mean AOD trend of −0.00023, and for MODIS Aqua, the mean AOD decreases from 0.113 to 0.075, with a monthly mean AOD trend of −0.00025.

The observed decrease in the AOD can be attributed to global and regional efforts aimed at enhancing air quality through policy interventions and technological advancements. Additionally, broader economic and societal changes within the region may also contribute to this trend.

4. Aerosol Types and Analysis of AE Biases Retrieved from MODIS at AERONET Stations

4.1. Classification of Aerosol Types and Individual Contribution to Atmospheric Aerosols

The commonly used AERONET aerosol classification, as described by Dubovik et al. (2002) and employed by Raptis et al. (2020) for the study of aerosol properties above Athens, relies on inversion products and data from stations with well-known dominant aerosol types to establish threshold values for each class [64,65]. However, the accuracy of these thresholds may be compromised in cases of complex aerosol mixtures. It is important to note that urban atmospheric mixtures rarely consist of a single aerosol type, requiring the classification to be constrained to characterising the prevailing type(s) closest to the recorded optical properties. This classification methodology combines direct measurements and inversion products, aiming to identify predominant aerosol types and their proportions. This is particularly relevant in urban environments, where a single dominant aerosol type is unlikely. Nevertheless, uncertainties in this classification, as highlighted by Dubovik et al. (2002), stem from the precision of individual retrievals, influenced by factors such as measurement errors, systematic instrumental biases, and assumptions related to aerosol properties [64].

Figure 9 and Figure 10 similarly demonstrate the aerosol classification throughout the CEE region, employing the methodology devised by Dubovik et al. (2002) and Raptis et al. (2020) [64,65], aimed at identifying the prevailing aerosol type among the 29 AERONET stations. In Figure 8, six distinct aerosol types (biomass or smoke, continental, dust, marine, mixed, and polluted) are discerned utilizing the calculated AE at 440/870 nm and the aerosol optical depth (AOD) at 440 nm. Furthermore, we conducted an analysis of the aerosol classification, examining monthly variations to observe the seasonal dynamics of different aerosol types across the CEE region as well as across urban and rural sites.

Figure 10 is comprised by three pie charts, which represent the aerosol composition in total (a), urban areas (b), and rural areas (c). Below these charts, bar graphs display the distribution of aerosols for each month corresponding to the total, urban, and rural analyses. Urban sites were selected based on proximity to cities with populations exceeding 50,000 within a 5 km radius, while all other stations were categorized as rural.

The statistical distribution of aerosol classes reveals a predominance of mixed aerosols, constituting 33% of the total, 34.9% in urban areas, and 31.4% in rural sites. During winter months, this percentage decreases in favor of the continental category. Generally, continental and polluted aerosols contribute to the total aerosol load with 22.9% and 22.2%, respectively; 20.3% and 24.6% in urban locations; and 24.8% and 20.7% in rural locations. Bar graphs indicate a significant contribution of polluted aerosols during summer months (June, July, August). Dust contributions primarily result from Saharan dust intrusions, peaking in spring (April, May) at nearly 20% of the monthly aerosol loading. On average, the annual dust contribution across all stations is around 8.5%.

Biomass-burning aerosols are predominantly visible during spring and summer, contributing approximately 5% to 10% of the total monthly aerosol loading and about 5.3% of the total aerosol loading, with 6.6% in urban locations and 4.2% in rural sites. The spring and summer peaks in biomass-burning aerosols were also observed by Barnaba et al. (2011) in Central and Eastern Europe between 2002 and 2007, highlighting seasonal consistency in multidecadal observations [31]. These aerosols originate mainly from forest fires and biomass burning, as indicated in Figure 7 in Folonchyk et al. (2020a) [44]. Marine aerosols exhibit a greater presence in rural areas, which include AERONET stations along the coast, accounting for 10.3%, compared to only 5.2% in urban locations. Monthly percentages decline during summer months from almost 20% to 5% in rural locations and from 10% to 2% in urban locations.

4.2. MODIS AE Biases for Different Types of Aerosols

Based on the six aerosol classes described by the aerosol typing methodology conducted in Section 4.1, we analysed the biases between the AERONET AE and MODIS Terra AE, see Figure 11. The MODIS AE 440/675 nm data were compared with the coincident collocated AERONET at 440/675 nm. The difference between the AERONET and MODIS AE shows a positive bias for biomass-burning, continental, mixed, and polluted types of aerosols, with the AERONET AE appreciably larger values than the MODIS AE (mean AE difference values range between 0.169 for continental aerosols to 0.671 for biomass-burning aerosols). The spread in bias is quite large, with a standard deviation (SD) varying between 0.282 to 0.57, but the confidence intervals (CI95%) are quite narrow around the mean. In other words, based on our type of classification, the MODIS AE underestimates these types of aerosols. In the case of marine aerosols, the mean bias is negative (mean = −0.262 and SD = 0.296), suggesting that the MODIS AE is larger than the AE retrieved by the AERONET, with a narrow CI95% between −0.282 and −0.241. In the case of dust aerosols, the mean AE difference is almost 0, with a large SD of ±0.51 and CI 95% between −0.031 and 0.031, being statistically inconclusive. These values suggest that based on our aerosol typing, the coincident MODIS AE inadequately classifies large particles. The MODIS AE parameter is obtained as a “derived” rather than a “retrieved” product of the retrieval algorithm and is typically considered unreliable in low AOD conditions, <0.2 [48]. Similar results were obtained for the biases calculated between the AERONET AE and the MODIS Aqua AE and are presented in the Supplementary Materials, Figure S4.

There are significant differences between the MODIS and AERONET AE that can be attributed to spatial and temporal sampling procedures, aerosol type classification, and model assumptions and retrieval algorithms of the satellite data.

5. Conclusions

Our study provides a comprehensive analysis of the aerosol optical depth (AOD) across Central–East Europe (CEE) for the period 2010–2023, focusing on seasonal variations and comparing data derived from MODIS satellite observations with ground-based measurements from 29 AERONET stations. Through this analysis, we aimed to assess the MODIS’s performance in capturing aerosol properties and variability and to identify potential biases or inconsistencies.

Comparing the MODIS AOD with the AERONET AOD allowed us to assess the consistency and accuracy of satellite-derived aerosol retrievals. Evaluating the agreement between these two datasets across different seasons provided valuable insights into the performance of the two types of remote sensing techniques. Discrepancies identified between MODIS and AERONET AOD measurements may point to biases or limitations in satellite retrieval algorithms, particularly under specific seasonal or environmental conditions, thus facilitating potential improvements in satellite retrieval techniques.

The analysis of AOD maps from MODIS retrievals alongside AERONET data provided a comprehensive overview of spatiotemporal variations in the AOD over CEE across different seasons. Seasonal trends in the AOD were evident, with the lowest values during winter and the highest during summer. We observed a persistent elevation of AOD levels throughout the year at specific AERONET stations, indicating potential local sources or specific regional transport patterns. Data gaps, particularly during winter, were attributed to meteorological factors and surface conditions impacting both ground-based and satellite measurements. Furthermore, our comparison between MODIS Terra and MODIS Aqua revealed similarities in seasonal patterns of AOD, although notable differences emerged, particularly during winter. The presence of bright surfaces covered by snow and the temporal offset between MODIS Terra and MODIS Aqua overpasses contributed to pronounced missing data during winter months in the region.

The intercomparison analysis of AOD data derived from MODIS satellite retrievals and AERONET stations provided a thorough exploration of spatiotemporal variability in the AOD over the region, offering valuable insights into the agreement and disparities between satellite and ground-based measurements. The evaluation of the AERONET AOD at 550 nm and the MODIS Terra AOD at 550 nm across four distinct seasons from 2010 to 2023 revealed seasonal variations in the correlation strength between the datasets, with notable improvements observed from winter to summer (from R² of 0.58 in winter to R² of 0.76 in summer). The similar analysis with MODIS Aqua AOD retrievals revealed lower correlation coefficients (from R² of 0.53 in winter to R² of 0.74 in summer) compared to the MODIS Terra–AERONET comparison, attributed partly to more missing data points during winter due to the later overpass time of MODIS Aqua. Seasonal variations in atmospheric turbulence and cloud coverage can influence the retrieval of the AOD and collocation with AERONET data. The statistical analysis revealed slight deviations in slope values between MODIS Terra and MODIS Aqua, indicating potential biases in AOD retrievals. Specifically, both MODIS sensors underestimated the AOD during winter, with MODIS Aqua exhibiting a larger deviation. However, in spring and autumn, both sensors showed good agreement with AERONET measurements, while in summer, MODIS Terra slightly overestimated the AOD compared to MODIS Aqua.

The examination of temporal biases in the AOD using the calculated estimated error (ER) between AERONET and MODIS retrievals reveals a notable seasonality in coincident retrievals, particularly with fewer available days during winter. In winter, MODIS Aqua’s ER values indicate poor agreement with AERONET measurements (when ER < −1), while for the MODIS Terra AOD bias determination, the ER values are within uncertainty bounds (−1 < ER < 1). Both sensors show higher positive biases against the AERONET in spring and summer compared to fall and winter. The consistently higher biases in Terra compared to Aqua are in line with previous analyses and persist throughout the analysed period. Furthermore, the analysis of AOD biases between AERONET and MODIS retrievals at individual AERONET stations revealed variations in biases across sites and sensors. MODIS Terra exhibits a bias range spanning from −0.086 to 0.037, while MODIS Aqua shows values ranging from −0.063 to 0.155. The results show that MODIS Terra mostly overestimates the AOD, while MODIS Aqua mostly underestimates the AOD compared to AERONET observations. Larger deviations in biases are attributed to factors such as limited matching retrievals, site-specific conditions, or inconsistencies in aerosol types. While most locations show good agreement between the AERONET and MODIS AOD, some exhibit asymmetrical distributions in the ΔAOD, indicating potential inconsistencies in retrievals. However, discerning the exact sources of the MODIS bias requires further investigation. These findings underscore the importance of considering both spatial and temporal factors when interpreting satellite-derived aerosol measurements.

In a series of recent investigations conducted by Folonchyk et al. (2020a and 2020b), an examination of the seasonality and trends of aerosol concentrations across ten Eastern European countries were investigated spanning the period from 2002 to 2019 [44,45]. This investigation employed data acquired from the MODIS Terra and Aqua satellites. The analyses show a discernible downward trend in the aerosol optical depth throughout the studied timeframe, showcasing annual trends ranging from −0.0050 in Belarus to −0.0029 in Russia. In our study, the AERONET aerosol trend analysis over the CEE region shows a monthly trend of −0.00027 and an annual trend of −0.0033 across the 10 countries where AERONET observations were collected. Similarly, the monthly MODIS Terra and Aqua AOD trends follow a decrease, with slopes of −0.00023 and −0.00025, respectively. The trends are supported by the associated p-values smaller than 0.05. These findings are in line with previous research studies. The observed decrease in the AOD is likely attributed to global and regional efforts in the past few decades aimed at enhancing air quality through policy interventions and technological advancements alongside broader economic and societal changes within the region.

Lastly, an aerosol classification method, which identifies the predominant aerosol types and their atmospheric contribution, was explored, acknowledging that urban atmospheric mixtures rarely consist of a single aerosol type. We characterised aerosol types across the CEE region using data from 29 AERONET stations. Six distinct aerosol types—biomass-burning or smoke, continental, dust, marine, mixed, and polluted—are identified based on the AE and AOD. Monthly variations reveal the seasonal dynamics of these aerosol types across urban and rural sites. The statistical analysis indicates a predominance of mixed aerosols, followed by continental and polluted aerosols. Polluted aerosols show significant contributions during summer months, while dust contributions peak in spring due to Saharan dust intrusions. Biomass-burning aerosols are more prevalent in spring and summer, originating from forest fires and biomass burning. Marine aerosols exhibit a greater presence in rural areas, with the lowest occurrence in summer months. Overall, these findings underscore the complexity of aerosol composition across the CEE region and highlight the importance of considering seasonal dynamics and urban–rural distinctions in aerosol studies.

Based on our aerosol classification method, we also compared the biases between the AERONET AE for the six different aerosol types with the corresponding collocated MODIS AE retrievals. We noticed a positive bias for biomass-burning, continental, mixed, and polluted aerosols, with the AERONET AE larger than the MODIS AE. The dust results were inconclusive, and the biases for marine aerosols presented a negative bias, with MODIS AE values larger than AERONET AE values. The confidence intervals calculated for the mean biases were narrow and proved the analysis significant.

Our study contributes to a better understanding of aerosol dynamics in the CEE region and highlights the challenges and opportunities associated with remote sensing techniques for monitoring aerosol variability across different seasons and environmental conditions. Further research and advancements in satellite retrieval algorithms are warranted to enhance the accuracy and reliability of aerosol retrievals, particularly in regions prone to seasonal variability and challenging surface conditions.