Two decades of high-resolution aerosol product over the Sierra Nevada Mountain region (SE Spain): Spatio-temporal distribution and impact on ecosystems

Atmospheric Research 308 (2024) 107515 Available online 5 June 2024 0169-8095/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Two decades of high-resolution aerosol product over the Sierra Nevada Mountain region (SE Spain): Spatio-temporal distribution and impact on ecosystems Ana del ´ Aguila a , b , * , Domingo Alcaraz-Segura c , d , e , Javier Martínez-L´ opez a , d , e , Thedmer Postma a , e , Lucas Alados-Arboledas a , b , Regino Zamora a , e , Francisco Navas-Guzm´ an a , b a Andalusian Institute for Earth System Research (IISTA-CEAMA), University of Granada, 18006 Granada, Spain b Department of Applied Physics, University of Granada, 18071 Granada, Spain c Department of Botany, Faculty of Science, University of Granada, 18071 Granada, Spain d Centro Andaluz para el Cambio Global - Hermelindo Castro (ENGLOBA), University of Almería, 04120 Almería, Spain e Department of Ecology, Faculty of Science, University of Granada, 18071 Granada, Spain ARTICLE INFO Keywords: Aerosol AOD Ecosystem Mountain MAIAC AERONET ABSTRACT Atmospheric aerosols play a pivotal role in shaping our environment, impacting climate, human health, and ecosystems. Characterizing the influence of aerosols on ecosystems, especially in mountain environments, is a challenging task due to their complex-orography and scarcity of aerosol ground stations. Satellite-based aerosol data can improve our knowledge over such complex-orography areas. Thus, we have analyzed the Aerosol Optical Depth (AOD) product from the MODerate resolution Imaging Spectrometer (MODIS) sensor produced by the inversion algorithm MultiAngle Implementation of Atmospheric Correction (MAIAC) over the last two decades for the period 2001–2022 with a spatial resolution of 1 ×1 km. Our study focuses on the Sierra Nevada Mountain region and National Park in Southeastern Spain. As a first step, we have validated the AOD from MODIS+MAIAC against three AERONET stations at different altitudes (680 m, 1800 m, and 2500 m above sea level (a.s.l.)). MODIS+MAIAC AOD showed good agreement with the ground-based AOD observations, with R values ranging from 0.75 to 0.82, RMSE values ranging from 0.047 to 0.066 and having 80% of the samples within the expected error (EE) of the product. The MODIS+MAIAC AOD product is able to characterize the finescale features of such complex-orography area and hence evaluate the spatio-temporal distribution of the AOD over the mountainous region. We have generated the most extended AOD dataset for a mountainous region, spanning the past two decades. We have deepened into the spatial and seasonal AOD patterns from 2001 to 2022, unveiling elevated AOD values near valleys and urban areas. In general, the AOD values decrease with increasing altitude with the exception of snow-covered areas at high altitudes (>2800 m a.s.l.), which might affect aerosol retrieval and provide bias due to higher-reflecting surfaces and pixel removal. For the first time, the relationship of aerosol loading with ecosystem type has been assessed in the protected environment of Sierra Nevada Natural Park. Monthly AOD trends across different ecosystem types and altitudinal ranges are analyzed in detail over the last two decades. In addition, Generalized Linear Models (GLM) are applied to reveal significant correlations between ecosystems and AOD, irrespective of altitude, latitude or longitude. Based on the interannual variation of AOD over the last two decades, we have analyzed the relationship of AOD with the different ecosystems of Sierra Nevada at 500 m elevation ranges. The patterns of the ecosystem's types are maintained over the elevation ranges 1200–1700 m and 1700–2200 m a.s.l., which demonstrates that land-type has an impact on the AOD product. Furthermore, it is observed that forest-like ecosystems tend to present lower AOD compared with baresoil or low-growth vegetation ecosystems. In addition, the areas of the mountain closer to Granada city present generalized higher AOD values on the western part of the mountain, regardless of the ecosystem, showing the significant influence of the proximity of urban sites over the ecosystems and the potential impact on the environment. * Corresponding author at: Andalusian Institute for Earth System Research (IISTA-CEAMA), University of Granada, 18006 Granada, Spain. E-mail address: [email protected] (A. del ´ Aguila). Contents lists available at ScienceDirect Atmospheric Research journal homepage: www.elsevier.com/locate/atmosres https://doi.org/10.1016/j.atmosres.2024.107515 Received 11 March 2024; Received in revised form 20 May 2024; Accepted 3 June 2024 Atmospheric Research 308 (2024) 107515 2 1. Introduction Atmospheric aerosols have a strong impact on the Earth radiative balance, climate change and the physicochemical processes (IPCC, 2023) that affect vegetation and ecosystems. Specifically, highmountain protected environments are of great interest in terms of biodiversity richness and atmospheric dynamics. Aerosols influence ecosystem photosynthesis and evapotranspiration by means of changes in climatic conditions (Zhou et al., 2021). Furthermore, aerosols can impact mountain environments in different ways: (1) by attenuating solar radiation reaching the surface (Li et al., 2020; Li et al., 2022; Chen et al., 2022), which might be a source of stress for plants and animals in high-mountains; (2) fertilization of terrestrial vegetation, which might induce an increase in aerosol depositions on vegetated-areas (JuanOvejero et al., 2022) or water ecosystems fertilization, particularly small-lake lentic systems; and (3) accelerating snow melting, affecting water provisions and regulation services. However, the impact of atmospheric aerosols on complex-orography areas like high-mountain ecosystems is scarcely investigated and presents a critical knowledge gap. In fact, several studies agree that the effects of aerosols on plants require research effort (e.g. Burkhardt and Grantz, 2016) as they can affect ecosystems by means of changes in climatic conditions. The Sierra Nevada Natural Park (SNNP) at Southeastern Spain is a high-mountain ecosystem and complex-orography area and hence, a biodiversity hotspot in the Mediterranean basin, which makes it an exceptional natural laboratory for ecological research within the longterm biosphere reserve (Cazorla et al., 2023). Furthermore, Sierra Nevada is close to Granada urban environment and numerous studies have investigated the impact of aerosols on the Granada city area (del ´ Aguila et al., 2018; Lyamani et al., 2008, 2010; Navas-Guzm´ an et al., 2013; Titos et al., 2017). However, aerosols over the mountain of Sierra Nevada on a regional scale have not yet been explored. Thus, we have investigated the Aerosol Optical Depth (AOD) at 550 nm as the variable that accounts for the extinction of solar radiation in the sampled atmospheric column attributed to aerosol particles. AOD is a columnar property that includes both anthropogenic and natural aerosols and whose temporal and spatial variability is governed by processes such as emissions or transport. Aerosols can be transported at highor low-level from local or remote areas. Dust aerosols play a major role in AOD quantification since it is the most abundant component in terms of mass (Zender et al., 2004). Dust intrusions are frequent in the region of study due to its proximity to North Africa, which is one of the major dust sources of the planet and thus, showing high dust optical depths close to the Mediterranean region next to North Africa (Gkikas et al., 2021, 2022). Ground-based aerosol measurements can provide accurate information; however, their punctual location prevents them from being representative of larger complex areas. Satellites can provide global coverage, allowing accurate monitoring of larger areas. The MODerate resolution Imaging Spectrometer (MODIS), along with the Multi-angle Implementation of Atmospheric Correction (MAIAC), provides a high spatial resolution AOD of 1 km ×1 km. The MAIAC algorithm is found to perform well in Europe (Di Antonio et al., 2023), in Central Asia (Chen et al., 2021) and globally (Qin et al., 2021). Also, it has been tested in specific regions such as the Po Valley, in northern Italy (Arvani et al., 2016) or in the C´ ordoba region in Argentina (Della Ceca et al., 2018). However, few studies have validated and evaluated this algorithm over mountainous regions. A previous study characterized the AOD in the mountains of the Alps at high spatial resolution with MAIAC v.5 (Emili et al., 2011). MAIAC Collection v.6 (Lyapustin et al., 2018) accuracy depended on the surface properties and comparisons revealed high confidence retrievals over cropland, forest, savanna, and grassland covers (Martins et al., 2017). However, the problem of residual cloud/ snow contamination within Collection v.6 was present and could provide high values over bright surfaces compared to over vegetated areas. Therefore, the AOD retrieval over high-altitude mountains could lead to high AOD values, which can be noted in the Alps region in a recent study (Di Antonio et al., 2023). The effect of high AOD values has been significantly improved for the AOD retrieval in the latest MAIAC Collection v6.1 (Lyapustin, 2022) and has been used in this research study. Furthermore, the latest collection of MODIS+MAIAC has allowed us to perform a long-term study over the last two decades from 2001 until 2022 in the complex-orography area of Sierra Nevada. One of the key aspects of this study is the evaluation of the uncertainties of satellite/model products in complex terrain using different ground-based remote sensing instruments (sun-photometers). Factors such as complex topography, variability in surface Bidirectional Reflectance Distribution Function (BRDF), presence of shadows, cloud screening and snow cover make the accurate aerosol retrieval from satellite observations challenging (Rao et al., 2022). These difficult conditions in capturing aerosol properties in such complex environments often lead to enhanced uncertainty in satellite-derived aerosol products over high mountainous areas, as reported in Lyapustin (2022). Thus, we have validated the AOD from MODIS+MAIAC against three AERONET locations at different altitudes: Granada (680 m above sea level (a.s.l.)), Cerro_Poyos (1809 m a. s.l.) and Albergue_UGR (2500 m a.s.l.). To the best of our knowledge, this is the first study that characterizes atmospheric aerosols in a complex-orography of a high-altitude mountain using MAIAC v6.1. Henderson et al. (2022) pointed out the importance of accurately representing land surface properties for characterizing the radiative characteristics in atmospheric convection. Thus, the land-cover type of Sierra Nevada can be relevant for aerosol retrieval. Aerosol retrieval from the MAIAC algorithm can be influenced by surface reflectance or by the convective effects of the atmosphere over different type of ecosystems. In this line, forests might present lower aerosol loading than other types of vegetation like shrublands or croplands. Long-term information on such complex-orography as the high-altitude mountains of Sierra Nevada constitutes an opportunity to understand the aerosol behavior in other complex-orography areas in the Earth. Thus, ecosystem types, aerosol loading, and interannual variability are studied for the long-term database over SNNP, including all types of aerosols. The aim of the study is to analyze the impact of ecosystems or land cover types on aerosols. The novelties of this study are as follows: (1) The AOD from MAIAC v6.1 has been used for characterizing the AOD at high spatial resolution over a complex-orography area; (2) The period of the analysis covers more than 20 years (2001−2022), being the longest climatology of AOD with satellite data for a mountain site to date; (3) the spatio-temporal distribution of AOD is assessed through the interannual mean maps and monthly interannual mean maps, summarizing the trends of AOD for the latest 20 years at high-spatial resolution; (4) the relationship between the ecosystem types of Sierra Nevada, elevation and AOD is examined. In this regard, we have used MODIS+MAIAC AOD at 550 nm product for the latest collection V6.1 to evaluate the spatial and temporal variability of aerosols in Granada and Sierra Nevada for the period 2001–2022. This multidisciplinary study enhances the role of functional diversity in high-mountain ecosystems and its impact on climate change through long-term analysis/climatology. Therefore, it is of great importance for future synergetic research to understand the interactions among aerosols, climate and vegetation. The study is organized as follows: Section 2 provides a description of the study area of Sierra Nevada (SE Spain), including the orography and a map of ecosystems. In addition, this section introduces the aerosol products for ground-based and satellite measurements. Section 3 describes the methodology employed for the analyses. The results are presented in Section 4 including the following subsections: (1) validation analysis, (2) trend analysis of AOD over the complex-orography, (3) interannual summaries of AOD for the period 2001–2022, (4) relationship between AOD and ecosystem types in the Sierra Nevada and (5) relationship between AOD, elevation and ecosystem types. Finally, the conclusions provide a summary of the results and a discussion about the interactions among aerosols, vegetation and climate change. A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 3 2. Data 2.1. Study area The study area is located in Southeastern (SE) Spain, within the highmountain protected area of Sierra Nevada in Granada, where comprises an extensive mountain range (Fig. 1a). The National Park of Sierra Nevada is the largest protected region in the whole of Spain. The core part is preserved as a National Park and is surrounded by the SNNP, which is a complex-orography area where the highest peak is Mulhacen with an altitude of 3478 m and the second-highest is the Veleta summit at 3398 m. The city of Granada is situated in a valley to the southeast of the mountain range, at 680 m above sea level (a.s.l.). The closest AERONET mountain stations are Cerro Poyos and Albergue_UGR at 1809 m and 2800 m a.s.l., respectively (Fig. 1b). Western to southern wind directions prevail in the Sierra Nevada Mountain range with respect to the city due to the so-called mountain breeze (del ´ Aguila et al., 2018) and favors the transport of pollutants from lower altitudes to the high-mountain altitudes (Casquero-Vera et al., 2020). Regarding its climatology, the main external aerosol source regions affecting the study are (1) Europe as a major source of anthropogenic pollution and (2) North Africa as a principal source of natural dust (Cazorla et al., 2017; Lyamani et al., 2005; Guerrero-Rascado et al., 2008, 2009; Valenzuela et al., 2012). Sierra Nevada is a natural laboratory and provides a wide variety of ecosystems that constitutes the Sierra Nevada Biosphere Reserve. The predominant ecosystem types in SNNP (Fig. 1c) are as follows: (1) natural pine tree forest, (2) oak tree forest and pine plantations, (3) highmountain shrubland and (4) agricultural areas, mid-mountain shrubland and grassland (Martínez-L´ opez et al., 2016). This valuable information belongs to the Global Change Observatory of the Sierra Nevada, which was established more than a decade ago (Zamora et al., 2017). 2.2. Aerosol products 2.2.1. AERONET dataset The ground-based dataset of aerosols is obtained from the Aerosol Robotic NETwork (AERONET) (Holben et al., 1998; Dubovik et al., 2000) over three measurement stations at different altitudes: Granada (680 m a.s.l.), Cerro_Poyos (1809 m a.s.l.) and Albergue_UGR (2500 m a. s.l.). Detailed information about the stations is described in Table 1. Granada is an urban environment station while Cerro_Poyos is a mid-altitude station (the predominant ecosystem type is mid-mountain shrubland) and Albergue_UGR is a high-altitude station (the predominant ecosystem types are urban, high-mountain grasslands and shrublands) located on the northern slope of the Sierra Nevada Mountain range (see Fig. 1a). These three AERONET stations are part of the Andalusian Global ObseRvatory of the Atmosphere (AGORA), which consists of several experimental sites on different atmospheric backgrounds in Southern Spain (https://atmosphere.ugr.es/en/about/presentat ion/agora, Laj et al., 2024). All the AERONET stations are equipped with sun photometers model CIMEL CE-318-4. A comprehensive description can be found in Holben et al. (1998). The variable of interest for this study is the AOD at 550 nm and Angstrom Exponent (AE) at Fig. 1. (a) Topographic map of the study area of Sierra Nevada at Southeastern Spain. The three black dots correspond to the three AERONET stations: Granada, Cerro_Poyos and Albergue_UGR; (b) Elevation and distance between the three AERONET stations; (c) Major ecosystem types in the SNNP (P´ erez-Luque et al., 2019). A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 4 440–870 nm for Level 2.0. The AERONET ground-based networks at Granada urban area have been described and used in several previous studies (e.g., Granados-Mu˜ noz et al., 2020; Valenzuela et al., 2017). The periods of available data vary depending on the AERONET station. At Granada urban station, the period goes from 2005 until present with very few interruptions. Regarding the midand high-altitude stations, the measurements are usually performed only during spring, summer and autumn because of the difficult human access to the stations during winter, due to high snow presence. Thus, the Cerro_Poyos station has measurements from 2011 to 2021, normally during periods from June to October, with an interruption in 2020 due to the pandemic. The Albergue_UGR station has measurements only for the summers of 2016 to 2018. Table 1 provides an overview of the number of days with available measurements per station. 2.2.2. MODIS+MAIAC Collection V6.1 For this study, the MCD19A2 V6.1 data product is used, which is the latest version of the MODerate resolution Imaging Spectrometer (MODIS) Terra and Aqua combined Multi-angle Implementation of Atmospheric Correction (MAIAC) Land Aerosol Optical Depth (AOD) gridded Level 2 product produced daily at 1 km resolution. For more information see MAIAC Data User's Guide Collection 6.1. Version 3.1. (Lyapustin, 2022). From the MCD19A2 V6.1 data product (MCD19A2.061), the AOD at 550 nm are used in this study. From now on, AOD will refer to the AOD at the wavelength of 550 nm. In order to obtain high-quality AOD data, some filters have been applied with the bitmask for AOD quality assurance variable (AOD_QA) of MCD19A2.061 product. According to the User's Guide, the best quality AOD is obtained with the combination of two filters: QA. CloudMask =Clear and QA.AdjacencyMask =Clear (Lyapustin, 2022). The data availability goes from February 2000 until present. In the guide of the latest collection, it is mentioned that due to crosscalibration, MAIAC processes MODIS Terra and Aqua jointly as a single sensor (Lyapustin, 2022). Therefore, we are using both sensors indistinctly. 3. Material and methods 3.1. Ground-based and satellite AOD measurements The AERONET stations have been used to validate MODIS+MAIAC from 2005 to 2022. In order to compare all the measurements for the AOD at the same wavelength of 550 nm, the Angstrom Exponent law (Schuster et al., 2006) has been used for computing the AOD from AERONET at 550 nm. The selected grids used for the comparison of MODIS+MAIAC are the minimum 1 ×1 km grid with the AERONET stations at the center of the pixel grid. The values of AERONET taken for the comparison are the average values in the range ±15 min from the MODIS overpasses. The interannual summaries of the AOD products have been processed using the Google Earth Engine (GEE) platform (https://eart hengine.google.com/). The GEE combines a multi-petabyte catalog of satellite imagery and geospatial datasets with planetary-scale analysis capabilities (Gorelick et al., 2017). For our study, we used the web-based Earth Engine code editor to develop algorithms for data collection and processing. We have processed the MODIS+MAIAC V6.1 dataset (MCD19A2.061) as the source for our interannual summaries. 3.2. Validation evaluation For the AOD validation, we have computed the following statistics: the Pearson correlation coefficient (R, Eq. (1)), the square of the Pearson correlation coefficient (R 2 ), square root-mean square error (RMSE, Eq. (2)), mean absolute error (MAE, Eq. (3)), bias (Eq. (4)) and expected error (EE) in Eq. (5), R=  1−∑ N i=1(AODi−AODAER,i)2 ∑ N i=1(AODi−AODAER,i)2 √ √ √ √ √ √ √ √ (1) RMSE = 1 N∑ N i=1(AODi−AODAER,i)2 √ √ √ √(2) MAE =1 N∑ N i=1 AODi−AODAER,i (3) bias =1 N∑ N i=1(AODi−AODAER,i)(4) EE = ±(0.05 +0.15⋅AODAER)(5) where N is the total number of coincident values, AODi corresponds to the AOD for MODIS+MAIAC and AODAER,i corresponds to the AERONET AOD, all of them at 550 nm. The EE represents that the accuracy of the retrievals is given by the envelope of ±(0.05 +0.15⋅AOD)and is within ±1 σ (Levy et al., 2010, 2013). The Global Climate Observing System (GCOS) recommends stricter requirements for quantifying the validation and is also adopted in the evaluation of Aerosol Climate Data Records (CDR) as explained in Popp et al. (2016). The GCOS fraction is the percentage of satellite-retrieved AOD satisfying the GCOS requirement of: GCOS =max(0.04 or 0.1AOD) (Chen et al., 2020). Moreover, we have performed the linear regression in order to obtain information about the uncertainty of the surface reflectance estimation with the intercept and the slope as an indicator of the aerosol model assumptions (Falah et al., 2021). 3.3. Generalized Linear Model To assess the effect of explanatory variables such as the altitude and the ecosystem types on the AOD, we have applied Generalized Linear Models (GLM). Thus, the AOD is considered a function of elevation, latitude, longitude and ecosystems, i.e., AOD =f(elevation, latitude, longitude, ecosystem i ) with i representing each ecosystem type. The goal is to determine the relationship between each ecosystem type and altitude with AOD. The Gaussian distribution for the GLM is considered. 4. Results and discussion 4.1. Validation of MODIS+MAIAC against AERONET Fig. 2 shows the validation analysis between MODIS+MAIAC against AERONET AOD at three stations at different altitudes: the urban station of Granada (680 m a.s.l.) and the two midand high-altitude stations Table 1 Description of the Aerosol Robotic Network (AERONET) ground-based stations. Station Coordinates Altitude (m a.s.l.) Period # days available Type of environment Granada 37.164◦N/ 3.605◦W 680 2005–2022 4290 Urban Cerro_Poyos 37.109◦N/ 3.487◦W 1809 2011–2021 1132 Mid-mountain Albergue_UGR 37.096◦N/ 3.387◦W 2500 2016–2018 291 High-mountain A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 5 called Cerro Poyos (1809 m a.s.l.) and Albergue_UGR (2500 m a.s.l.), respectively. The product used from MODIS+MAIAC is the AOD combined Terra and Aqua of 1 km resolution over land. They show good correlations with AERONET observations, with R values ranging from 0.75 to 0.82, which is in agreement with previous studies (Su et al., 2022; Ye et al., 2022; Yousefi et al., 2020). The RMSE obtained in this study ranges from 0.047 to 0.069 which is lower than that of previous studies and others in the literature such as in East Asia described in Wang et al. (2022). Moreover, the RMSEs obtained in this study are similar to the RMSE =0.06 of previous validation studies (Di Antonio et al., 2023; Lyapustin et al., 2018; Martins et al., 2017). It is also observed that in the midand high-mountain stations, the overestimation of AOD is predominant, with >EE between 15 and 18%, respectively. However, around 80–82% of the AOD retrievals fall within the EE, which is considerably above the one standard deviation (⁓67%) threshold, that is the needed envelope threshold required for a successful validation (Levy et al., 2010). The GCOS fraction for the three stations is ⁓70%. Thus, the statistical variables obtained in this study Fig. 2. Correlation between MODIS+MAIAC AOD against AERONET AOD at 550 nm at three AERONET stations. (upper panel) Scatter plots of AOD for the validation of the three AERONET stations at different heights. The color bar indicates the density of points. Solid gray lines correspond to the EE =±(0.05 +0.15⋅AOD). Dashed gray lines correspond to the linear fit of each dataset. (middle panel) Idem as the upper panel but with the color bar referring to the AE (440–870 nm) of the corresponding AERONET station. The orange and blue solid lines correspond to the linear fit that accomplish AE <1.0 and AE >1.0, respectively. Solid black lines of all figures correspond to the line 1:1. (bottom panel) Probability density functions for each AERONET station, including the fitting parameters of mean value ( μ ) and standard deviation ( σ ). A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 6 indicate that MODIS+MAIAC show good agreement with AERONET observations for all stations. The probability density function of the differences between (MODIS+MAIAC) – AERONET is depicted in Fig. 2 (bottom panel) showing the dependency of the differences on the AOD levels. The Granada station exhibits mean values ( μ ) and standard deviation ( σ ) of −0.0075 and 0.0652, respectively. In contrast, for the mountain stations, the μ values are around 0.02 and σ around 0.04. Furthermore, the results show significant spatial agreement between the MODIS+MAIAC dataset for the three different altitude locations with slopes close to 0.6. The value of the slopes is in agreement with several studies (e.g. Falah et al., 2021). Nevertheless, certain differences can be observed in the AOD values in relationship with the altitude of the station. Overall, the aerosol loading is lower for higher altitude stations for both AERONET and MODIS+MAIAC datasets. The AERONET AE (440–870 nm) serves as an indicator of aerosol particle size, with values below 1 suggesting a relative prevalence of coarse mode aerosols, while higher values indicate a predominance of fine mode aerosols. Based on the scatter plot of Fig. 2 (middle panel), the AOD has been separated into colors corresponding with the AE of each AERONET station for fine particles (AE >1.0) and coarse particles (AE <1.0), i.e., dust-dominated. The latter provides a clear separation of the two regression lines based on the size of the aerosol, given by the AE information. The intercept lines are below 0.05 for all stations. Thus, from the comparison between AERONET and MODIS+MAIAC, it can be observed that fine particles (AE >1.0) have a slope between 0.72 and 0.82 while coarse particles (AE <1.0) have a slope between 0.44 and 0.56. The information of AE, combined with the expected errors of >EE and <EE of AOD, indicate that the major part of overestimated AOD corresponds to the fine-mode (mainly urban pollution), while the underestimated AOD corresponds to coarse-mode aerosols (dust-dominated type). The latter result is similar to that found in Su et al. (2022) for dust-dominated AOD types. These results indicate that the MAIAC algorithm is sensitive to aerosol size (Di Antonio et al., 2023; Rogozovsky et al., 2023). Therefore, MODIS+MAIAC provides persistently high AOD values for fine particles compared to AERONET, whereas the opposite occurs for coarse aerosols. It is worth noticing that as we increase in altitude for the AERONET stations, the slopes change accordingly. Hence, for the highest-altitude station of Albergue_UGR, the differences between the slopes of fine and coarse particles are minimum while for Granada station, the differences are maximum. Ultimately, the validation results are satisfactory, confirming the representativeness of MODIS+MAIAC measurements for the urban station of Granada, and the midand highaltitude mountain stations of Cerro Poyos and Albergue_UGR, respectively. Table 2 describes the statistical analysis performed for the different AERONET stations and the different measurement periods. This analysis is carried out for coincident periods in order to evaluate the differences based on altitude only, regardless of the amount of data. The following conclusions can be extracted from Table 2: - The R values do not vary significantly among periods and Granada station statistics are very similar throughout the different periods. - For coincident periods, the station of Cerro Poyos (mid-altitude location) provides R ≥0.82 which are the highest correlations compared to the other stations and periods. In general, the agreement between MODIS+MAIAC and the different AERONET stations is very good and we can use the former to describe the AOD variability over different altitudes. This result is analyzed in detail in the following section. 4.2. Statistical analysis over the long-term AOD series The long-term series analysis of AOD for MODIS+MAIAC covers the period 2001–2022 for the three locations while AERONET covers the period 2005–2022 at Granada station. The temporal evolution of the monthly averages and the decomposition of these series into trend, seasonality and residuals is shown in Fig. 3. The method employed is the “Seasonal and Trend decomposition using Loess” (STL) which is a robust method for decomposing time series (Cleveland et al., 1990). Monthly averages are computed for the months with at least 10 days of data. Fig. 3a shows that there are frequent gaps in the winter months for the high-altitude station of Albergue_UGR because MAIAC does not retrieve AOD over snow (Lyapustin et al., 2018; Lyapustin, 2022). For the STL analysis, we have applied linear interpolation for missing values to avoid data gaps. The trends shown in Fig. 3b indicate that the AOD for MODIS+MAIAC displays a “quasi” flat trend for the three stations. However, the AOD values are on average ⁓22% higher at Granada than at both mid−/high-altitude stations. On the other hand, the AOD from AERONET displays a downward trend, and yet there are periods with variations that might show a natural cyclical pattern with 3–5 years fluctuations: the AOD dramatically drops from 0.18 to 0.11 in the period 2005–2010, then climbs to 0.15 in 2013, and immediately thereafter the trend decreases until 2018. Finally, a new increasing trend began in 2019. The higher generalized AOD values of AERONET compared to MODIS+MAIAC suggest that the ability of the latter to capture coarse aerosol loadings such as dust or large particles is lower compared to AERONET. This effect is more noticeable in the trends shown in Fig. 3b. In fact, the increasing trend of AOD from 2020 until 2022 is only captured by AERONET due to a sharp increase in Saharan dust intrusions (Cuevas-Agull´ o et al., 2024). The seasonal patterns of AOD displayed in Fig. 3c indicate strong seasonal patterns with a cycle of 12 months for MODIS+MAIAC, with strong peaks occurring in June for Granada and Cerro Poyos and in May for Albergue_UGR. Similarly, negative peaks occur in December for Granada and in November or October annually for the mountain stations. On the other hand, the seasonal patterns of AERONET present strong peaks, generally occurring in June, with semi-annual seasonal patterns with a weak peak in January, annually. However, the patterns range over time, in coincidence with the trend behavior, i.e., from 2005 to 2010 the strong peaks occur in different months with no peak in 2009. Then, from 2010 to 2013 the strong peaks occur in March while from 2013 to 2014 it occurs in May and from then onwards, it occurs in June. The analysis of seasonal patterns reveals that MODIS+MAIAC adequately captures AOD peaks from AERONET at the urban environment of Granada, with the exception of periods characterized by wider data gaps, as illustrated in Fig. 3a. The persistent high AOD values in Granada, as depicted in the temporal series, trend analysis and seasonality (Fig. 3 a-c), align with expectations. Granada is located in an urban environment with multiple Table 2 Statistics of the validation analysis for each AERONET station against MODIS+MAIAC and the corresponding periods of measurements (number of points, N): 2016–2018 (261), 2011–2021 (1251) and 2005–2022 (4759). The statistical values analyzed are: R, R2, MBE and RMSE. Granada (680 m a.s.l.) Cerro Poyos (1809 m a.s.l.) Albergue_UGR (2500 m a.s.l.) 2016–2018 R 0.736 0.837 0.782 R 2 0.542 0.701 0.611 MBE RMSE −0.006 0.062 0.012 0.049 0.024 0.047 2011–2021 R 0.750 0.823 R 2 0.563 0.678 MBE −0.007 0.018 RMSE 0.062 0.050 2005–2022 R 0.754 R 2 0.568 MBE −0.007 RMSE 0.066 A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 7 anthropogenic sources, such as traffic and heating systems, leading to higher aerosol contamination levels. In contrast, the higher altitude stations experience lower aerosol loading due to reduced pollution in those areas. Furthermore, the higher AOD values over the city of Granada are in agreement with local studies at Granada station that found that the topography of Granada basin, surrounded by mountains produced difficult ventilation processed and favored aerosol stagnation (Patr´ on et al., 2017). The residuals shown in Fig. 3d are almost constant over the longterm period except for AERONET, where more outliers appear. These outliers occur more frequently during the winter and summer seasons, suggesting extreme values in real observations. This pattern is likely influenced by higher pollution levels during winter and dust outbreaks during summer. 4.3. AOD characterization at high spatial resolution over a complexorography area In this section, we present the MODIS+MAIAC v6.1 AOD interannual statistic at each pixel for more than two decades comprising 2001–2022 in the SNNP and surrounding areas. The aim is to describe the spatial heterogeneity and temporal variability of the AOD in a complexorography and highly diverse protected area at high-altitude. Therefore, we have analyzed the interannual mean and interannual standard deviation (SD) of the mean annual maps for the period 2001–2022. Fig. 4 shows the spatial patterns of the interannual mean and standard deviation (SD) of the annual mean AOD for the period 2001–2022 over SNNP and its surroundings. The interannual mean (Fig. 4a) shows higher AOD over the valley close to Granada city, flat low altitudes and the western drylands of the Natural Park, with a strong vertical gradient due to elevation. This result aligns with previous research findings for mountainous areas in Nepal, with higher aerosol concentrations found at lower elevations (Dhital et al., 2022). Previous MAIAC collections were affected by land-surface brightness, and were prone to overestimating AOD as surface brightness increased, resulting in a positive bias (Lyapustin et al., 2018). This effect can be appreciated with this collection to a minor extent. However, inside the Natural Park (delimited area in Fig. 4), the areas with the highest AOD values concentrate over altitudes above 2800 m, which are spatially coincident with the ecosystem type of “High-mountain grasslands and rocks” (refer to Fig. 1c). This alignment suggests a plausible relationship between elevated surface reflectance caused by rocks and/ or snow, influencing the retrieval values of AOD. This effect is more noticeable on the very high interannual SD (Fig. 4b) at the high-altitude areas, since the strong interannual variability in snow-cover typical of a Mediterranean mountain strongly affects albedo and, subsequently, AOD retrieval. Moreover, interannual SD revealed some artifacts in MAIAC estimated AOD (small pixel clusters displayed as red/orange Fig. 3. Time series decomposition plots of AOD at 550 nm using the STL technique. The AOD is decomposed into four components. (a) Temporal series of the monthly averages, (b) trends, (c) seasonality and (d) residuals for MODIS+MAIAC v6.1 dataset at the three locations of Granada (green), Cerro_Poyos (orange) and Albergue_UGR (purple) and for AERONET for Granada station (gray). A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 8 Fig. 4. Spatial patterns of the interannual maps of the (a) mean and (b) standard deviation (SD) of the AOD product from MODIS+MAIAC Collection 6.1 for the period 2001–2022 in SNNP complex-orography region. The black line indicates the boundaries of the SNNP while the black dots indicate the AERONET stations. Fig. 5. Interannual monthly maps of the mean AOD product from MODIS+MAIAC Collection 6.1 for the period 2001–2022 in the SNNP for each month. A. del ´ Aguila et al. Atmospheric Research 308 (2024) 107515 9 squares with remarkably high AOD SD, >0.015), which correspond to human land-use changes, solar plants or water areas (dams), whose subsequent changes in land-surface brightness might alter the retrieval algorithm and create artifacts. Fig. 5 shows the spatial patterns of the interannual mean of AOD from MAIAC v6.1 over the last 20 years for each month. In colder months, there are gaps at the highest altitude because pixels completely covered with snow are filtered out. Furthermore, the surroundings of those pixels are noisier because of residual snow and cloud contamination over high elevations, which provides high albedos and eventually reduces the number of available pixels for valid retrievals. In fact, in months like December or May, there might be snow presence to a lower extent at the highest altitude regions leading to errors in the retrievals and hence, resulting in higher AOD values close to the summits. Therefore, there are more noisy pixels coincident with high altitudes in cold months due to the effect of the topography together with cloud and snow contamination. In warmer months, the AODs are generally higher than those in colder months, with smoother spatial patterns. The AODs increase around 0.1 on average compared to colder months, which might be influenced by the Planetary boundary layer height (PBLH) seasonal patterns and the convective layer reaching the Sierra Nevada Mountain at different altitudes (Reche et al., 2022). On the one hand, the altitudinal patterns show enhanced AOD values in the lower valleys and flatter regions and lower AODs at higher altitudes over bare-soil or rocky terrains. The mean PBLH annual mean over the Granada region is 1700 ±500 m a.s.l., although the PBLH shows seasonal cycles, with higher values in summer and spring (Granados-Mu˜ noz et al., 2012). In fact, the maximum convective boundary layer height over Granada in summer and spring is 2000 ±700 and 1600 ±500 m a.s.l, respectively (Moreira et al., 2020). On the other hand, the higher PBLH during the warmer months might favour the aerosol transportation from the city to the mountain by increasing the AOD values over higher elevations during the warmer months. The relationship between MODIS+MAIAC AOD and PBLH over the months is consistent with the results found in a study at the Po Valley (Arvani et al., 2016). There is a noticeable “North-South line” in April and May at the approximate longitude of 3.15 o W, which coincides with a transitional zone of 300 km in the new Collection 6.1 of AOD and may introduce weak spatial gradients across the boundary (Lyapustin, 2022). In summary, the interannual monthly cycles are well captured by the AOD from MAIAC v6.1. In addition, the shape of the mountain range topography is evidenced by the spatial gradient in the AOD distribution for all seasons, with generally lower values at higher altitudes than in low/mid-altitude regions. The high correlations between MODIS+MAIAC and AERONET described in the previous section indicate a good ability of the former to reproduce the AOD monthly variability. 4.4. Relationship between AOD and ecosystems in Sierra Nevada Ecosystems are spatially structured, so environmental variables are often spatially autocorrelated, which can affect the statistical signification of estimates (Alcaraz-Segura et al., 2013; Liras et al., 2008). In this section, we examine the relationship between AOD and different ecosystems present in the SNNP through statistical metrics, significant correlations and annual dynamics. Table 3 shows a summary of the statistical parameters per ecosystem type together with the area covered and altitudinal ranges. The ecosystem types High-mountain meadows, Autochthonous Scot pine, Autochthonous pine over dolomites and Aquatic systems cover such small areas and thus, they have been neglected for the subsequent analysis to avoid biased results. The highest mean AOD value corresponds to High-mountain grasslands and rocks, which aligns with the highest elevated regions above 2700 m a.s.l. where soils are scarcely vegetated and large shiny micaschist outcrops are exposed, which are highly reflective surfaces and hence, might alter MAIAC retrieval by enhancing the AOD values. The lowest mean AOD is found for the Holm oak forest, which is a semideciduous forest with a dense dark green canopy that can trap and avoid dust emissions, but that might also influence the algorithm retrieval due to its low albedo. Fig. 6 shows the interannual monthly evolution of MODIS+MAIAC AOD for the different ecosystems (solid lines) and for the three pixels where the AERONET field stations (dashed lines) are located. The interannual monthly means correspond to the average for each month of every monthly mean for each year. The AOD shows a clear seasonal dynamic, with lower values in colder months and reaching maxima in July, except for the High-mountain grassland ecosystem, whose peak occurs in June. The MODIS+MAIAC AOD at the Granada urban station is above all the ecosystem types and the other two AERONET stations for all months. This result is expected, as there is more pollution in the city compared to the mountain sites. The seasonal pattern of AOD for the ecosystems shows one peak in July whereas in the metropolitan area, there are two peaks: one in May and another one in August, which might be due to the higher frequency of desert dust outbreaks during spring and summer seasons (Cazorla et al., 2017; Navas-Guzm´ an et al., 2013). Thus, generally higher AOD values are detected during spring/summer periods, which is in agreement with previous studies that show a similar dynamic for AOD (e.g. Arvani et al., 2016). Cerro Poyos station is 5 km from an urban area and is also in the middle of the Sierra Nevada border of dolomites, which are characterized by their whiteness, with scarce vegetation and with presence of quarries that often raise dust. Cerro Poyos shows a similar monthly dynamic as Granada, but with lower AOD values, which is consistent with the results found in Section 4.1 and 4.2. Regarding Albergue_UGR, the AOD values are the lowest for all months and there is only one peak Table 3 Summary of statistics of the interannual mean of annual AOD means from MODIS+MAIAC over the ecosystem types in Sierra Nevada. For each ecosystem type the mean ±SD, minimum, maximum, area and altitudinal range are given. Ecosystem type Mean ±SD Min Max Area (km 2 ) Altitudinal range (m) All ecosystems 0.088 ±0.010 0.071 0.128 2794 287–3307 High mountain grasslands and rocks 0.101 ±0.002 0.079 0.128 233 2102–3307 High-mountain meadows 0.089 0.089 0.089 1 1791 High-mountain shrubland 0.090 ±0.002 0.073 0.119 611 1567–2958 Mid-mountain shrubland 0.090 ±0.003 0.074 0.121 708 298–1941 Mid-mountain grasslands 0.087 ±0.003 0.075 0.116 86 293–2203 Autochthonous Scot pine 0.088 ±0.001 0.087 0.089 2 1465–1739 Autochthonous pine over dolomites 0.095 ±0.003 0.090 0.103 6 1043–1374 Pine plantations 0.084 ±0.003 0.071 0.122 639 800–2417 Holm oak forest 0.081 ±0.001 0.072 0.108 196 1079–1951 Pyrenean oak 0.087 ±0.001 0.073 0.103 63 1005–1901 Mountain crops 0.090 ±0.002 0.075 0.112 245 287–1895 Aquatic systems 0.091 ±0.001 0.079 0.111 4 728–1799 A. del ´ Aguila et al.