ENSO and Light-Absorbing Impurities and Their Impact on Snow Albedo in the Sierra Nevada de Santa Marta, Colombia

: Snow albedo is an important variable in the coupled atmosphere-earth system at the global level. Moreover, studying its behavior allows us to know the state of the cryosphere. The Sierra Nevada de Santa Marta (SNSM) is a glacier area and the northernmost tropical (10.82 ◦ N, 73.75 ◦ W) region in South America. It has a height of up to 5775 m.a.sl., which is the second highest mountain in the world near the marine coast. We analyzed variations in snow albedo related to snow cover, snowfall, temperature, light-absorbing impurities such as blank carbon (BC), organic carbon (OC) and dust, and El Niño—Southern Oscillation (ENSO) phenomenon through 20 years (2000–2020). We mainly use daily data from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra and Aqua NASA satellites. Results showed through correlations that snow albedo has decreased due to Land Surface Temperature (55%, p < 0.001), a positive phase of ENSO (42%, p < 0.001) and dust (37%, p < 0.01) in the SNSM. Additionally, a dust negative e ﬀ ect was more evident on the southern side (up to 49%, p < 0.001) of the SNSM. Backward trajectories by the NOAA HYSPLIT model suggest that dust sources would be soil erosion in the surrounding region. Results can help recognize the inﬂuence of ENSO and dust in the glacier decrease of the SNSM. the NDSI showed a positive correlation with albedo of up to 45% ( p < 0.001) for the November-December (ND) period in the Southeast area. Furthermore, SP showed a correlation with SA of up to 75% ( p < 0.001) for the ND period in the Northern area. SA showed a negative correlation of 36% ( p < 0.001), 42% ( p < 0.001), and 41% ( p < 0.001) for LST, ENSO ( + ), and dust, respectively, in the ND period in Southwestern area. Our data indicate that the greatest negative impact on SA was produced by dust in the Southeast area with 48% ( p < 0.001) during the period from May to June (MJ) and by ENSO ( + ) in the Southwestern area with 42% ( p < 0.01) for the ND period. Overall, it was observed that BC and OC did not have a single negative e ﬀ ect on SA, it seems to have e ect three + + Dust)


Introduction
The cryosphere is the frozen part of the hydrosphere and the lithosphere due to its physical interactions with the atmosphere. The cryosphere (as snow and ice) plays an important role in biogeochemical cycles and in the global energy balance [1][2][3][4]. Thus, any qualitative and quantitative change in the physical properties and extent of the cryosphere can affect global air circulation, air and ocean temperatures, river flow in mountains, ocean current patterns, and sea level [5][6][7][8][9][10]. In fact, this indicates that studying the cryosphere is important to better understand these ecosystems and as an indicator of climate change at a global and local level [11,12].
Snow albedo is an important variable in determining the amount of solar radiation adsorbed on the cryosphere. Furthermore, it is defined as a relationship between the incoming and reflected solar radiation by a snow-covered surface. For example, fresh snow can almost completely reflect incoming solar radiation (0.95), but can decrease to 0.45 due to metamorphism [13][14][15]. Snow albedo variations are mainly favored by surface temperature, snowfall, snow age, and snow impurities [14,[16][17][18][19]. Light-absorbing impurities (LAI) on snow reduce snow albedo and absorb more solar radiation Geosciences 2020, 10, 437 3 of 21

Analyzed Area
SNSM is the northernmost tropical cryosphere area in Colombia and South America (10.82 • N, 73.75 • W), as shown in Figure 1 [60]. It is geologically associated with the Tropical Andes [65]. Moreover, SNSM is the highest coastal mountain system in the world (5775 m.a.sl.). Its highest peaks are less than 45 km from the coastline in the Caribbean Sea [66]. Since 1964, it is a protected area called Sierra Nevada de Santa Marta National Natural Park, covering 3830 km 2 . Its highest areas are inhabited by the Koguis, Arhuacos, Wiwas, and Kankuamos indigenous communities, who also administer and help preserve the SNSM from the Kogui-Malayo-Arhuaco reservation as a sub-area within the mentioned national natural park. These communities limit access to areas of snow and glaciers, as they are considered sacred areas for their culture [67]. Furthermore, the SNSM serves as a source of water to 36 rivers, which are the main source of fresh water to approximately 3.4 million inhabitants in the departments of Magdalena, La Guajira, and Cesar from North Colombia [60, 68,69].

Analyzed Area
SNSM is the northernmost tropical cryosphere area in Colombia and South America (10.82° N, 73.75° W), as shown in Error! Reference source not found. [60]. It is geologically associated with the Tropical Andes [65]. Moreover, SNSM is the highest coastal mountain system in the world (5775 m.a.sl.). Its highest peaks are less than 45 km from the coastline in the Caribbean Sea [66]. Since 1964, it is a protected area called Sierra Nevada de Santa Marta National Natural Park, covering 3830 km². Its highest areas are inhabited by the Koguis, Arhuacos, Wiwas, and Kankuamos indigenous communities, who also administer and help preserve the SNSM from the Kogui-Malayo-Arhuaco reservation as a sub-area within the mentioned national natural park. These communities limit access to areas of snow and glaciers, as they are considered sacred areas for their culture [67]. Furthermore, the SNSM serves as a source of water to 36 rivers, which are the main source of fresh water to approximately 3.4 million inhabitants in the departments of Magdalena, La Guajira, and Cesar from North Colombia [60, 68,69].

Data Analysis
We analyzed a data set from March 2000 to February 2020. We selected the areas with a height greater than 4000 m.a.s.l. (Error! Reference source not found.) because, according to previous studies, this is where snowfall occurs and glaciers are located [57]. We mainly use MODIS daily data (https://modis.gsfc.nasa.gov/data/) with a special resolution of 0.5 and 1 km, as shown in Error! Reference source not found.. The latest MODIS data collection (V6) includes several enhancements to the previous version including, but not limited to, elimination of Terra sensor degradation issues and improvements to atmospheric calibration [70]. In addition, snow products such as MYD10A1 and MOD10A1 were compiled with algorithms that only include the highest quality atmospherically

Data Analysis
We analyzed a data set from March 2000 to February 2020. We selected the areas with a height greater than 4000 m.a.s.l. (Figure 1) because, according to previous studies, this is where snowfall occurs and glaciers are located [57]. We mainly use MODIS daily data (https://modis.gsfc.nasa.gov/data/) with a special resolution of 0.5 and 1 km, as shown in Table 1. The latest MODIS data collection (V6) includes several enhancements to the previous version including, but not limited to, elimination of Terra sensor degradation issues and improvements to atmospheric calibration [70]. In addition, snow products such as MYD10A1 and MOD10A1 were compiled with algorithms that only include the highest quality atmospherically corrected observations [71]. Previous studies indicate that the general error of the albedo product in snow (MOD10A1 and MYD10A1) can vary from 1% to 10% for good observations with low atmospheric disturbances [43]. In satellite imagery analysis, the steep slopes of complex terrain can increase errors. In this sense, this type of error is likely higher in the SNSM due to complex terrain. However, changes in snow albedo can still be quantified, so we opted for a temporal change analysis per pixel that is expected to mitigate the influence of topography to some extent by avoiding direct intercomparing of pixels influenced by different slopes-aspects [17,39]. The MODIS instrument products do not generate data in the presence of clouds in order to minimize the effect of missing data under these conditions. We combined the data recovered from the Terra and Aqua satellites, by means of the daily pixel by pixel comparison to replace the failed data or overcast sky. This was observed during the pass of the Terra satellite but not by the Aqua satellite and vice versa. Since there was data for both, the one with the best quality was selected [72]. This is because these satellites have the MODIS instrument on board and pass through the studied area daily, but, at different times, increase the probability of taking cloud-free observations. For example, when only using Terra (MOD10A1), data were retrieved for 69% of the analyzed time period. However, when combined with data retrieved from Aqua (MYD10A1), the cloud-free data increased to 83%. Snow precipitation (SP) were estimated (per pixel) for the rains retrieved from the GPM satellites (https://gpm.nasa.gov/missions/GPM) and LST below 0 • C (MOD11A and MYD11A1) [17,41]. For instance, whenever a pixel registered a data of rainfall and a temperature below 0 • C at the same time, then it is assumed that it is SP [17,41]. We used snow cover as Normalized difference snow index (NDSI), which is an index that is connected to the presence of snow in a pixel and is a more accurate description of snow detection. Then, snow cover was detected using the NDSI ratio of the difference in VIS and SWIR reflectance. NDSI = ((band 4 − band 6)/(band 4 + band 6)). Thus, a pixel with NDSI > 0.0 is considered to have some snow present. At the same time, a pixel with NDSI <= 0.0 is a snow free land surface [91]. Daily data sets of SA, NDSI, SP, and LST were analyzed using box-and-whisker plots. With this, it is possible to observe the variability and skewness of the data and, if so, its positive or negative direction. In this way, the trend of the data can be estimated, either increasing or decreasing in the period analyzed. In addition, we analyzed the SA monthly variations with the bi-monthly ENSO index (MEI) (https://psl.noaa.gov/enso/mei/), and with the reanalysis data of BC, OC, and dust from MERRA-2 (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/) [92]. SA mean monthly values were correlated by means of Pearson's coefficient (r) with mean monthly values of NDSI, LST, SP, warm phase of ENSO as ENSO (+), ENSO cold phase as ENSO (−), BC, OC, and dust. Finally, we explore the potential sources of dust, as the most influential LAI in the SA, which may be reaching the snowy areas and reducing snow albedo. Then, we show the map of the degree of soil erosion in Northern Colombia [93]. Additionally, we used the web version of the NOAA HYSPLIT trajectory model driven by meteorological reanalysis data from Global Data Assimilation System (GDAS, https://www.ready.noaa.gov/HYSPLIT_traj.php). The multiple backward of air parcels arriving to SNSM have been computed for 24-h of Caribbean air mass transport [94]. Snow precipitation (SP) were estimated (per pixel) for the rains retrieved from the GPM satellites (https://gpm.nasa.gov/missions/GPM) and LST below 0 °C (MOD11A and MYD11A1) [17,41]. For instance, whenever a pixel registered a data of rainfall and a temperature below 0 °C at the same time, then it is assumed that it is SP [17,41]. We used snow cover as Normalized difference snow index (NDSI), which is an index that is connected to the presence of snow in a pixel and is a more accurate description of snow detection. Then, snow cover was detected using the NDSI ratio of the difference in VIS and SWIR reflectance. NDSI = ((band 4 − band 6)/(band 4 + band 6)). Thus, a pixel with NDSI > 0.0 is considered to have some snow present. At the same time, a pixel with NDSI <= 0.0 is a snow free land surface [91]. Daily data sets of SA, NDSI, SP, and LST were analyzed using boxand-whisker plots. With this, it is possible to observe the variability and skewness of the data and, if so, its positive or negative direction. In this way, the trend of the data can be estimated, either increasing or decreasing in the period analyzed. In addition, we analyzed the SA monthly variations with the bi-monthly ENSO index (MEI) (https://psl.noaa.gov/enso/mei/), and with the reanalysis data of BC, OC, and dust from MERRA-2 (https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/) [92]. SA mean monthly values were correlated by means of Pearson's coefficient (r) with mean monthly values of NDSI, LST, SP, warm phase of ENSO as ENSO (+), ENSO cold phase as ENSO (-), BC, OC, and dust. Finally, we explore the potential sources of dust, as the most influential LAI in the SA, which may be reaching the snowy areas and reducing snow albedo. Then, we show the map of the degree of soil erosion in Northern Colombia [93]. Additionally, we used the web version of the NOAA HYSPLIT trajectory model driven by meteorological reanalysis data from Global Data Assimilation System (GDAS, https://www.ready.noaa.gov/HYSPLIT_traj.php). The multiple backward of air parcels arriving to SNSM have been computed for 24-h of Caribbean air mass transport [94].

Daily Variation of Snow Albedo, NDSI, Snow Precipitation, and Temperature
Error! Reference source not found. shows boxplots for daily SA, NDSI, SP, and LST datasets. SA showed a mean of 41.4% while a mean of 42.6%, indicating that the data is skewed towards higher values (mean > median) with outliers of up to 100% shown in Error! Reference source not found.a. This skewness toward positive values is related to the occurrence of higher mean daily SP of 7.46 mm with extreme values or outliers of up to 83 mm (Error! Reference source not found.c). However, the snow cover, represented by the NDSI (Error! Reference source not found.b) showed a skewness towards lower values with a mean of 0.07 and a median of 0.1 (mean < median). This decrease in snow cover would be influenced by the snowmelt caused by the increase in temperature with a median of 0.97 °C, a mean of 1.31 °C, a minimum value of −8.45 °C, and a maximum value of 9.6 °C (Error! Reference source not found.d).  Figure 3 shows the mean monthly variation from March 2000 to February 2020 of SA and SP. Figure 3a show a high consistency between the monthly behavior of the SP and SA from March 2000 to February 2020. It is appreciated that the SNSM has a bi-modal cycle of rainfall with peaks in May and October of each year (Figure 3b). In addition, it seems to influence an increase of the SA, which presents the maximum monthly mean during the same months of May and October. Meanwhile, Figure 4a displays the monthly average variation for SA and the bi-monthly multivariate ENSO index (MEI) from 2000 to 2020. It shows that, for most years, high snow albedo values were related to the negative phases (ENSO cold phase) and vice versa.

Impact of ENSO and LAI on SA
Error! Reference source not found. shows the mean monthly variation from March 2000 to February 2020 of SA and SP. Error! Reference source not found.a show a high consistency between the monthly behavior of the SP and SA from March 2000 to February 2020. It is appreciated that the SNSM has a bi-modal cycle of rainfall with peaks in May and October of each year (Error! Reference source not found.b). In addition, it seems to influence an increase of the SA, which presents the maximum monthly mean during the same months of May and October. Meanwhile, Error! Reference source not found.a displays the monthly average variation for SA and the bi-monthly multivariate ENSO index (MEI) from 2000 to 2020. It shows that, for most years, high snow albedo values were related to the negative phases (ENSO cold phase) and vice versa. Error! Reference source not found.b,c showed that the values of BC and OC show that the higher values were related to reductions in SA. However, no significant negative relationships of BC and OC were observed in SA. Error! Reference source not found.d shows the monthly behavior of snow albedo and dust surface mass concentration. It is shown that high snow albedo values were associated with low concentrations of surface dust. While the combined effect of BC, OC, and dust seems to have a significant effect in reducing the SA of the SNSM (Error! Reference source not found.e). Thus, considering the rain cycles observed in Error! Reference source not found.b, and the   Figure 4d shows the monthly behavior of snow albedo and dust surface mass concentration. It is shown that high snow albedo values were associated with low concentrations of surface dust. While the combined effect of BC, OC, and dust seems to have a significant effect in reducing the SA of the SNSM (Figure 4e). Thus, considering the rain cycles observed in Figure 3b, and the studied sub-areas shown in Figure 1, we estimate in Table 2 the r coefficient between SA and the analyzed variables. The SA variables showed significant correlations (positive or negative), according to the expected physical behavior. For example, the NDSI showed a positive correlation with albedo of up to 45% (p < 0.001) for the November-December (ND) period in the Southeast area. Furthermore, SP showed a correlation with SA of up to 75% (p < 0.001) for the ND period in the Northern area. SA showed a negative correlation of 36% (p < 0.001), 42% (p < 0.001), and 41% (p < 0.001) for LST, ENSO (+), and dust, respectively, in the ND period in Southwestern area. Our data indicate that the greatest negative impact on SA was produced by dust in the Southeast area with 48% (p < 0.001) during the period from May to June (MJ) and by ENSO (+) in the Southwestern area with 42% (p < 0.01) for the ND period. Overall, it was observed that BC and OC did not have a single negative effect on SA, but it seems to have a combined effect when the sum of the three (Bc + Oc + Dust) was considered. example, the NDSI showed a positive correlation with albedo of up to 45% (p < 0.001) for the November-December (ND) period in the Southeast area. Furthermore, SP showed a correlation with SA of up to 75% (p < 0.001) for the ND period in the Northern area. SA showed a negative correlation of 36% (p < 0.001), 42% (p < 0.001), and 41% (p < 0.001) for LST, ENSO (+), and dust, respectively, in the ND period in Southwestern area. Our data indicate that the greatest negative impact on SA was produced by dust in the Southeast area with 48% (p < 0.001) during the period from May to June (MJ) and by ENSO (+) in the Southwestern area with 42% (p < 0.01) for the ND period. Overall, it was observed that BC and OC did not have a single negative effect on SA, but it seems to have a combined effect when the sum of the three (Bc + Oc + Dust) was considered.
Geosciences 2020, 10, x FOR PEER REVIEW 9 of 22 Error! Reference source not found. shows the monthly mean and variation of the snow albedo for the whole SNSM area and for the subareas shown in Error! Reference source not found.. As mentioned before, snow albedo in the SNSM is related to snow precipitations that have two annual peaks (Error! Reference source not found.b). In general, the whole area presented maximum values of monthly mean snow albedo (46.5% in May and 46% in October) for each year (Error! Reference source not found.a). The highest monthly values of SA presented were 48.6%, 46.34%, and 45% in May for the northern, southwestern, and southeast areas, respectively. Generally, recovery of SA during the months of May is strongly related to SP, which shows significant positive correlational values between 53% (p < 0.001) and 64% (p < 0.001) for MJ periods (Error! Reference source not found.). In addition, the positive influence of NDSI on SA was more evident in the second annual peak presented in the months of September to December, which showed a significant correlation between 33% (p < 0.01) and 58% (p < 0.001) for the September to October (SO) period, as previously appearing in Error! Reference source not found..
The largest negative monthly variations of SA were in June and December for all areas analyzed (as shown in Error! Reference source not found.). The positive phase of the ENSO also shows a negative correlation with SA mainly in the second half of the year (June to December). Furthermore, dust was significantly negatively correlated with SA in SNSM. It showed correlations of 37% (p <  Table 2. Empirical results of the Pearson correlation coefficient (r) between snow albedo and the variables analyzed. Where *, ** and *** denote significant correlations at p < 0.05, p < 0.01, and p < 0.001, respectively.   Figure 5 shows the monthly mean and variation of the snow albedo for the whole SNSM area and for the subareas shown in Figure 1. As mentioned before, snow albedo in the SNSM is related to snow precipitations that have two annual peaks (Figure 3b). In general, the whole area presented maximum values of monthly mean snow albedo (46.5% in May and 46% in October) for each year (Figure 5a). The highest monthly values of SA presented were 48.6%, 46.34%, and 45% in May for the northern, southwestern, and southeast areas, respectively. Generally, recovery of SA during the months of May is strongly related to SP, which shows significant positive correlational values between 53% (p < 0.001) and 64% (p < 0.001) for MJ periods ( Table 2). In addition, the positive influence of NDSI on SA was more evident in the second annual peak presented in the months of September to December, which showed a significant correlation between 33% (p < 0.01) and 58% (p < 0.001) for the September to October (SO) period, as previously appearing in Table 2.  The largest negative monthly variations of SA were in June and December for all areas analyzed (as shown in Figure 5). The positive phase of the ENSO also shows a negative correlation with SA mainly in the second half of the year (June to December). Furthermore, dust was significantly negatively correlated with SA in SNSM. It showed correlations of 37% (p < 0.01), 41% (p < 0.001) and 43% (p < 0.001) in the ND period for the whole area, southwestern area, and southeastern area, respectively. Assuming that the dust sources are nearby to the area analyzed, Figure 6 shows the degree of soil erosion around SNSM, where it can be seen that, in most of this region on the Colombian side, the degree of soil erosion is very high or high. This likely causes dust particles from disturbed and loose soils to be dragged by the winds and reach the snowy areas of the SNSM. When exploring the sources of dust in the snow and glaciers of SNSM, we simulated a backward trajectory using the web version of the NOAA HYSPLIT model via the Real-time Environmental Applications and Display system (READY) [94,95]. We run 24 h backward trajectories for each of the months which the greatest impact of dust was observed in SA (Table 2), namely, September, October, November, and December. The results show that superficial air masses pass through the eroded area shown in Figure 6.

Area
Thus, these air masses would be lifting and dragging the dust to the snowy areas of the area studied in the SNSM, entering mainly from the northeast side (Figure 7).   [93]. Each degree has been established according to the level of degradation of the A horizon soil layer. The A horizon is known as the top layer of the mineral soil horizons, often referred to as topsoil. This layer is important because it contains dark decomposed organic matter (humus). Then, this figure represents the degree of soil erosion by total loss of the A horizon (very high), loss greater than 75% of the A horizon (high), loss between 25% and 75% of the A horizon (moderate), and loss of less than 25% of the A horizon  [93]. Each degree has been established according to the level of degradation of the A horizon soil layer. The A horizon is known as the top layer of the mineral soil horizons, often referred to as topsoil. This layer is important because it contains dark decomposed organic matter (humus). Then, this figure represents the degree of soil erosion by total loss of the A horizon (very high), loss greater than 75% of the A horizon (high), loss between 25% and 75% of the A horizon (moderate), and loss of less than 25% of the A horizon (low) [96]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).
Geosciences 2020, 10, x FOR PEER REVIEW 13 of 22 (low) [96]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Influence of NDSI and Snow Precipitation on Snow Albedo
Flanner and Zender [97] showed that snowfall provides fresh snow with small snow grains, increasing the total effective surface area of snow grains, and resulting in increased snow albedo. Other studies displayed results in the same direction and also linked high values of snow albedo to

Influence of NDSI and Snow Precipitation on Snow Albedo
Flanner and Zender [97] showed that snowfall provides fresh snow with small snow grains, increasing the total effective surface area of snow grains, and resulting in increased snow albedo. Other studies displayed results in the same direction and also linked high values of snow albedo to a longer time of snow cover permanence [17,18,40,[98][99][100][101][102]. In the inner tropics, precipitation shows seasonal changes due to the Intertropical Convergence Zone that ensures two wet seasons each year [103]. Furthermore, studies showed an increasing trend in rainfall, both on an annual scale and during the rainy season, after the mid-20th century, mainly north of 11 • S [104,105]. This is consistent with our findings. For instance, for almost all the analyzed area, SP and NDSI showed higher significant positive correlations with SA during the April to June and September to December periods (Table 2).

Impact of Temperature and ENSO on Snow Albedo
Rabatel et al. [57] analyzed glaciers of the tropical Andes of Colombia, Venezuela, Ecuador, Peru, and Bolivia. They found that the higher frequency of ENSO (+) events since the late 1970s, along with a warming of the troposphere over the tropical Andes, could explain much of the shrinkage of these tropical glaciers. They found that the monthly mass balance in glaciers was up to 3.5 times more negative during ENSO (+) events than in an average month in Colombia. Recent studies found similar data in tropical Andean glaciers such as Conejeras Glacier in Colombia. It shows great sensitivity to temperature changes and especially to the evolution of the ENSO(+) phenomenon with great loss of mass and area during these warm events, suffering a reduction of 37% between 2006 and 2017 [106]. Additionally, the average summer temperature near the surface is 0.7 to 1.3 • C higher during ENSO(+) than during ENSO(−), which improves the melting glacier surface [107]. These studies reinforce our results. In effect, Figure 2 shows an LST increase and significant negative correlations between ENSO(+) and SA (see Figure 5). As mentioned before, SA in the SNSM is related to SP that has two annual peaks (Figure 3b). In general, the whole area presented maximum values of monthly mean of SA (46.5% in May and 46% in October) for each year (Figure 5a). The highest monthly SA means presented values of 48.6%, 46.34%, and 45% in May for the northern, southwestern, and southeast areas, respectively. Generally, recovery of SA during the months of May is strongly related to SP, which shows significant positive correlational values between 54% (p < 0.001) and 64% (p < 0.001) for MJ periods ( Table 2). In addition, the positive influence of NDSI on SA was more evident in the second annual peak presented in the months of September to December, which showed a significant correlation between 33% (p < 0.01) and 58% (p < 0.001) for the SO and ND periods, as previously presented in Table 2. Thus, our analysis shows how the increase in frequency of ENSO(+) and its temperature risings are associated with the glacier mass loss identified in the SNSM between 2000 and 2018 [61].

LAI Effect on Decreasing of Snow Albedo
BC, OC, and dust are the main LAI particles that snow albedo affects through the darkening of these surfaces. The consequent absorption of more solar radiation, therefore, increase the snow melting rate [26,29,31,34,37]. Recent studies carried out in India analyzed the inorganic ions, OC, BC, dust, and sea salt through concentrations of PM 2.5 Surface Concentrations from MERRA-2. They found correlations of 0.96 and 0.6 between estimates and measurements for the monthly and daily averaged concentrations, respectively [84]. Thus, this shows that MERRA-2 data are useful in long-term air quality studies. Dust has been considered as a particle that scatters light in the atmosphere, but a particle that absorbs light only when it is deposited on snow or ice [6]. In fact, in many parts of the world, it has been identified that dust decreases the snow albedo [32,49,[108][109][110][111]. Studies conducted in the San Juan Mountains, CO, USA, showed that dust in snow has direct and indirect radiative impacts. Direct absorption by dust contributes approximately 80% of the total radiative forcing, while thickening of the grain of snow represents approximately 20%, and, consequently, a decrease in snow albedo [112]. These studies are consistent with our results. For example, our results revealed a statistical highly significant correlation (p < 0.001) that relates the snow albedo decrease with the dust increase of up to 41% and 43% for the sub-areas located at the southwestern and southeast, respectively (Table 2). Droughts have increased and land alternations in regional sources are related to the increases in dust deposition on mountain snow cover [113], enhancing absorbed solar radiation and its snowmelt rates [18,114]. Our results suggest that soil erosion in the surroundings of the SNSM ( Figure 6) and the dragging of exposed dust from its disturbed soil is blown by wind to the snowy studied areas (Figure 7), which is greater in the southern areas (as shown in Table 2) where the slope favors the atmospheric dust deposition. Other studies also point in this direction, since the reduction of snow cover duration is amplified where the mountainous areas receive dust from disturbed lands [108]. Therefore, under global warming and associated desertification by soil erosion, this process could threaten water resources fed by thawing in arid and semi-arid regions. Our results also displayed that BC and OC would not have a significant effect on the reduction of SA in the SNSM. This may be due to the fact that the main source of BC and OC are forest burns that have occurred sporadically with the duration of a few days in this region [115]. Therefore, the monthly averages of BC and OC since MERRA 2 were not sensitive to these changes.
Our findings show strong correlation between dust and the decrease in snow albedo in the SMSN during the last two decades (2000-2020). However, it presents some limitations. Daily concentrations of BC, OC, and dust were not analyzed. Especially, this did not allow us to observe the short effect that biomass burning would have in increasing LAI particles, such as BC and OC in snow albedo decrease in the SNSM. In addition, due to the limitations to access the snow and glacial areas imposed by the local indigenous communities [67], it has not been possible to collect snow and ice samples to analyze the presence of LAI particles. Future research can analyze the daily scale events of BC and OC increases and their effects on snow albedo variations in the SNSM. Furthermore, the possible effect that the arrival of the Saharan dust intrusion that are presented each year would have in snow albedo decrease of the studied area could be studied.

Conclusions
It is known that snow albedo changes are associated with climatic variables and negative feedback from light-absorbing atmospheric aerosols. Therefore, our research analyzed the relationship of snow cover, snowfall, temperature, and light-absorbing impurities such as black carbon, organic carbon, and dust on snow albedo in the SNSM during the last two decades (2000-2020). As expected, our studies showed that increases in NDSI snow cover and snowfall improved snow albedo in the SNSM, especially in the months of May and October of each year. However, increases of temperature and ENSO (+) frequency were statistically highly correlated with snow albedo by up to −55% and −42% (p < 0.001), respectively. Furthermore, our findings suggest that the presence of light-absorbing impurities such as dust in snowy areas would increase the reduction of snow albedo, especially in the sub-areas located in the south of the SNSM. Therefore, the negative correlation is between 34% (p < 0.01) and 49% (p < 0.001), which is caused mainly by dust from soil erosion in the nearby region. Overall, this study displayed that variations in snow cover, snowfall, temperature, and dust have had a significant correlation on the behavior of snow albedo in the SNSM. Thus, considering the importance of snowy areas through tropical mountains, this study can help decision-makers to better understand the changes under climate change scenarios.

Abbreviations
The following abbreviations are used in this manuscript.