The Moderate-Resolution Imaging Spectroradiometer (MODIS) installed on the Terra and Aqua satellites is one of the most widely used remote sensing platforms in earth science research (Parkinson, 2003). The Terra and Aqua satellites feature 5-min temporal and 36-channel spectral resolutions, 2330-km viewing swath widths, wide spectral ranges (0.412–14.24 μm), and near-global coverage every 1–2 days. Thus, MODIS not only provides reliable and extensive retrieval data products on aerosols and clouds but also offers an ideal solution for aerosol model development and validation, dynamic analysis of atmospheric pollution, and air quality monitoring with a spatial resolution of 1–10 km (Levy et al., 2007; Remer et al., 2008). Dark target (DT) and deep blue (DB) are the two major algorithms applied to the MODIS instrument for retrieving AOD. DT manifests its excellent performance over high vegetation coverage regions such as farmlands or forests, while DB was developed to retrieve the aerosol properties over Gobi, deserts, or snow-covered areas. The retrieval deviation can be reduced in the estimation of surface albedo and further in AOD over the regions of bright surface with a weaker blue-band reflectance (Hsu et al., 2006). In this study, snow cover and AOD at 550 nm were derived from MODIS-Aqua version 6.1. The resolution of snow cover is 0.05° × 0.05° (lat x lon), and the AOD at 550 nm is 1° × 1° (lat x lon).

The Multiangle Imaging SpectroRadiometer (MISR) onboard Terra was launched on 18 December 1999. The sensors can measure particles with medium spatial resolution and can cover almost the whole world. MISR products have made significant progress in retrieving aerosols. However, the predictability of MISR aerosol products for PM_2.5 concentrations has not improved (Liu et al., 2007). Both MISR and MODIS are onboard the Terra satellite. Statistical comparisons were made with coincident AOD retrieved by MODIS and MISR. The correlation coefficient between MISR and MODIS is approximately 0.9 over the ocean and approximately 0.7 over land (Kahn et al., 2009). Compared to AERONET, marine stations have the highest correlation coefficient (0.9), and dusty sites have the lowest correlation coefficient (0.7) (Kahn et al., 2005). Although MISR has these biases, MISR can still represent the spatial and temporal variations in aerosols. In this study, 550 nm AOD at 550 nm with a horizontal resolution of 0.5° × 0.5° (lat x lon) was derived from the MISR Level 3 Data.

The Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) is produced by the NASA Global Modeling and Assimilation Office (GMAO), which provides datasets dating back to 1980 (Gelaro et al., 2017). Previous studies have evaluated MERRA-2. Compared with models and observations, the MERRA-2 AOD is larger than the simulated AOD, but observations are larger than the MERRA-2 AOD (Randles et al., 2017). Furthermore, previous studies compared observed AOD based on AERONET and found that the MERRA-2 AOD better matches that of AERONET (Randles et al., 2017). Liu et al. (2019) combined the AAI and several satellite datasets to check the reliability of long-term changes in the MERRA-2 dust concentration. Moreover, MERRA-2 underestimates the daily temperature and overestimates the latent heat flux in summer (Draper et al., 2017). Although MERRA-2 still presented these biases, MERRA2 could represent the spatial and temporal variations in aerosols. In this study, we used AOD, dust aerosol optical depth (DAOD), and surface albedo with a horizontal resolution of 0.5° × 0.625°(lat x lon) to analyze the impacts of dust over the TP.

The European Centre for Medium-range Weather Forecasting Reanalysis V5 (ERA5) is a global atmospheric reanalysis dataset (Hersbach et al., 2019) dating back to 1950, with a horizontal resolution of 0.25° × 0.25° (Hersbach et al., 2018). Compared with ERA-Interim, ERA5 was significantly improved. The assimilation method uses the ensemble 4DVar data assimilation scheme with 10 ensemble members, with a horizontal resolution of 31 km × 31 km and a vertical stratification of 137 layers with a top pressure of 1 Pa. The integrated forecast system version was upgraded to Cy41r2, the radiative transfer mode was upgraded to RTTOV-v11, and the temporal resolution was improved to 1 h. The ERA5 reanalysis datasets are better than MERRA-2 in estimating solar irradiance, but there are some biases in cloudy conditions. ERA5 overestimates solar radiation, while MERRA-2 underestimates solar radiation in most areas. In terms of wind speed, ERA5 is superior to MERRA-2 compared to the observation data. In terms of precipitation, ERA-5 performs better on the monthly scale, followed by JRA-55, MERRA-2, and CFS-2 (Taszarek et al., 2021) in terms of thermodynamic theoretical parameters, low-altitude decrement rate, low-altitude wind shear, etc. (Maa et al., 2021). Both ERA5 and MERRA-2 data better represent variables such as temperature, moderateness, mid-tropospheric decrement rate, and mean wind. However, ERA5 has a higher correlation and lower mean error compared to MERRA-2 (Taszarek et al., 2021). In general, ERA5 is better than MERRA-2 data in terms of solar radiation, precipitation, temperature, and wind speed. However, ERA5 excludes aerosol data; we use MERRA-2 data for analyzing aerosols in this article. For the analysis of temperature and radiation, ERA5 data are used. In this study, we used surface temperature and surface net solar radiation (clear) with a horizontal resolution of 0.25° × 0.25° (lat x lon).

By using the correlation analysis method, this article analyzes the correlation relationship between MERRA-2 AOD and MODIS, MISR AOD, verifying the feasibility of MERRA-2 AOD in the study area. The correlation coefficient as defined by Pearson can be used to characterize the correlation between two different variables. Assuming a sample size of n for variables x and y, the correlation coefficient formula is as follows: r = ∑ i = 1 n [ ( x i − x ¯ ) ( y i − y ¯ ) ] ∑ i = 1 n ( x i − x ¯ ) 2 × ∑ i = 1 n ( y i − y ¯ ) 2 (1) where x ¯ and y ¯ represent the average of x and y, respectively.

The significance test for the correlation coefficient r can be performed using the t-test method, assuming that the two variables are not correlated, with the statistic t: t = r n − 2 1 − r 2 (2) n is the sample size and n-2 is the degree of freedom. The significance level is assumed to be α. If t < t_α, the linear correlation is insignificant. If t > t_α, the linear correlation is significant.

Assuming a sample size of n for a certain climate variable xi, and ti represents the time corresponding to xi, establishing a one-dimensional linear regression equation between xi and ti: x ^ i = a + b t i , i = 1,2 , ⋯ , n (3)

In Eq. 3, a is the regression constant and b is the regression coefficient. For the observed data xi and the time ti, the least-squares estimates of the regression coefficient b and the constant a are: { b = ∑ i = 1 n x i t i − 1 n ( ∑ i = 1 n x i ) ( ∑ i = 1 n t i ) ∑ i = 1 n t i 2 − 1 n ( ∑ i = 1 n t i ) 2 a = x ¯ − b t ¯ (4) where x ¯ = 1 n ∑ i = 1 n x i , t ¯ = 1 n ∑ i = 1 n t i

Using the relationship between the regression coefficient b and the correlation coefficient a, finding the correlation coefficient between time ti and the variable xi: r = ∑ i = 1 n t i 2 − 1 n ( ∑ i = 1 n t i ) 2 ∑ i = 1 n x i 2 − 1 n ( ∑ i = 1 n x i ) 2 (5)

The correlation coefficient r indicates the closeness of the linear correlation between the variable x and time t. When r = 0, the regression coefficient b is 0, and the regression line determined by the least-squares estimation is parallel to the x-axis, indicating that the change of x is independent of time t; when r > 0, b > 0, indicating an upward trend of x with time t; when r < 0, b < 0, indicating a downward trend of x with time t. In this article, we mainly use this method to calculate the linear trend of snow and dust aerosols.

This article focuses on calculating the rate of dust change by using statistical methods. The difference between the dust concentrations of two adjacent days (the latter-day minus the previous day) is calculated using dust concentration data from MERRA-2 reanalysis datasets, and the rate of dust change concentration is obtained by dividing the difference by the dust concentration of the day. The main calculation formula is as follows: D u s t r a t e ( i ) = D u s t ( i + 1 ) − D u s t ( i ) D u s t ( i ) (6) where Dust is the atmospheric dust concentration and Dust_rate is the diurnal rate of atmospheric dust concentration. i represents the time series, i represents the current day, and i+1 represents the next day.

The dust concentration is influenced by the season; the effect of seasonal variability of dust concentration needs to be removed for a clearer analysis of the effect of dust on snow. We use the difference transform to remove the seasonal variability. To specifically analyze the effect of high and low values of dust variability on snow, the rates of dust change data from 2009 to 2018 were arranged from smallest to largest in spring, and the high and low values were selected using the 5% and 95% quartiles. In this article, we calculated the 5% and 95% quartiles as −0.63 mg m⁻³ d⁻¹ and 0.65 mg m⁻³ d⁻¹. Therefore, all dates greater than 0.65 mg m⁻³ d⁻¹ are identified as high-value days, while dates less than −0.63 mg m⁻³ d⁻¹ are identified as low-value days. Furthermore, the highest diurnal rates of atmospheric dust concentration and the lowest diurnal rates of atmospheric dust concentration are usually located next to each other. When the diurnal rate of atmospheric dust is high, the dust concentration in the atmosphere gradually increases at this time, mainly because most of dust is released into the atmosphere. When the diurnal rate of atmospheric dust is low, the dust deposits on snow at this time. Thus, the high values represent the atmospheric dust and the low values represent dust-on-snow. High and low-value days are analyzed separately for radiation, temperature, and snow distribution to comprehend the mechanism of the effect of dust to snow.

As a major aerosol over the TP, dust in the atmosphere and snow on the TP remarkably influence the Asian climate. In this article, we choose the most dominant absorbing aerosols (dust aerosols) to estimate the relative contributions of atmospheric dust and dust-on-snow to the change in snow cover over the TP through the use of reanalysis datasets and satellite retrievals. We use three aerosol-related datasets, MERRA-2 reanalysis datasets, MODIS, and MISR retrievals. The MERRA-2 reanalysis datasets include dust-related variables, which are used to investigate the association between dust aerosols and snow cover, and they further analyze the influence of atmospheric dust and dust-on-snow to the change in snow cover over the TP. Due to differences in models and emissions, MERRA-2 AOD shows specific uncertainties in different regions and should be evaluated before use. We need to analyze the applicability of MERRA-2 reanalysis datasets on the TP. To evaluate the feasibility of MERRA-2 reanalysis datasets over the TP, the aerosol products from MERRA-2 are compared with those derived from MODIS and MISR retrievals. Overall, MERRA-2 captures the spatial and temporal distributions of 550 nm AOD over the TP for 2009–2018. The monthly mean AOD distribution from MERRA-2 is similar to that from the MODIS and MISR retrievals. The TP AOD from MERRA-2, MODIS, and MISR has the largest values in April (0.31 ± 0.03, 0.36 ± 0.04, and 0.23 ± 0.03), and December (0.09 ± 0.01, 0.08 ± 0.01, and 0.06 ± 0.01) has the smallest (Figures 1A–C). The MERRA-2 AOD, MODIS AOD, and MISR AOD have the largest values in the northern TP (0.25 ± 0.1, 0.30 ± 0.12, and 0.20 ± 0.08) and gradually decrease from the periphery of the TP to the inside and from north to south over the TP. Meanwhile, the TP AOD from MERRA-2 has been found to be greater than that from MODIS and MISR retrievals. The differences between MERRA-2 and MODIS may be related to cloud contamination effects and emission uncertainty in each model system (Liu and Chen, 2020) (Figures 1D–F). The cloud contamination effects are a source of bias in the MODIS retrieval of AOD, and the emission uncertainty refers to the MERRA-2 reanalysis.

Moreover, the MERRA-2 reanalysis datasets are consistent with the MODIS and MISR retrievals, which could reflect the seasonal and monthly variations in the TP AOD during 2009–2018 (Figure 2). Figure 2 shows that the values of the TP AOD from the MERRA-2 reanalysis datasets, MODIS, and MISR retrievals are 0.20 ± 0.08, 0.16 ± 0.09, and 0.14 ± 0.07 from 2009 to 2018, respectively. The TP AOD from the MERRA-2 reanalysis datasets has the largest values of 0.29 ± 0.03 in spring, followed by 0.23 ± 0.03 in summer, and the lowest AOD of 0.12 ± 0.02 in winter, which are consistent with the MODIS AOD values (0.31 ± 0.05, 0.15 ± 0.04, and 0.12 ± 0.02) and MISR AOD values (0.21 ± 0.02, 0.19 ± 0.03, and 0.08 ± 0.01) (Figure 2A). Although MERRA-2 AOD has been found to be greater than that from MODIS and MISR retrievals, the seasonal and monthly variation characteristics of three datasets are basically consistent. MERRA-2 AOD data can be used to reflect the variation characteristics of AOD in the TP. Figures 2B,C show the comparison between monthly MERRA-2 AOD and MISR AOD, MODIS AOD from 2009 to 2018. The comparison results show better correlations and smaller errors among the MERRA-2 reanalysis datasets, MISR, and MODIS retrievals, with correlation coefficients (R) of 0.90 and 0.76 and root-mean-square errors (RMSEs) of 0.03 and 0.06 (Figures 2B,C). Overall, the MERRA-2 monthly average AOD is in good agreement with MISR and MODIS.

As the main type of aerosol over the TP, dust AOD has the largest value of 0.15 ± 0.02 in spring and the lowest value of 0.03 ± 0.007 in winter, showing gradual decreases from the northern to southern regions of the TP (Figure 3A). The temporal and spatial distributions of DAOD (Figures 3A,C) and dust concentration (Figures 3B,D) are similar to that of AOD over the TP (Figure 1). DAOD (0.15 ± 0.02) and dust concentration (0.2 ± 0.06 mg m⁻³) had the largest values from March to June, while DAOD (0.03 ± 0.007) and dust concentration (0.04 ± 0.01 mg m⁻³) had the lowest values from October to February (Figures 3A,B). The monthly variations in DAOD in the northern TP diverge from those in the southern TP, with the DAOD in the northern TP having the largest values in May (0.28 ± 0.03) and the DAOD in the southern TP having the largest values in May (0.11 ± 0.016) (Figure 3).

The dust occurrence frequency in the northern TP is higher than that in the southern TP. The DAOD (0.25 ± 0.1) and dust concentration (0.35 ± 0.2 mg m⁻³) in the Qaidam Basin are higher than those in other regions over the TP. The difference between the northern TP and southern TP is closely related to the atmospheric circulation and aerosol emissions in the surrounding areas (Figures 3C,F). The results are consistent with those from Xu et al. (2015b), who showed the three-dimensional structure of aerosols over the TP based on MISR and Cloud-aerosol Lidar and Infrared Pathfinder Satellite (CALIPSO) retrievals. They further found that the spatial distribution of dust aerosols in the northern TP is significantly different from that in the southern TP, which can be seen clearly from 6 to 8 km above sea level, especially in spring and summer. Chen et al. (2013) and Yuan et al. (2020) also investigated the meridional transport path of dust from the TD to the northern slope of the TP and found that the westerly jet weakens and moves northward providing dynamic and thermal conditions for meridional transport.

The spatial distribution of dust concentration in spring and winter (Figures 4A,C) is similar to that of MERRA2 DAOD (Figures 1D, 4B,D). The dust concentration in the northern TP (0.5 ± 0.1 mg m⁻³) was higher than that in the southern TP (0.15 ± 0.05 mg m⁻³) in spring from 2009 to 2018. A large dust concentration value (0.6 ± 0.1 mg m⁻³) occurred in the Qaidam Basin, and a small ±value (0.05 ± 0.01 mg m⁻³) was located east of 95°E and south of 30°N on the TP based on MERRA-2 reanalysis datasets (Figure 4A). The spatial distribution in winter is similar to that in spring; due to the frequent dust storms in spring, the dust concentration in winter (0.15 ± 0.02 mg m⁻³) is smaller than that in spring (Figure 4B). In addition, based on Eq. 3, the regression coefficients of monthly data with time series were calculated for spring and winter from 2009 to 2018 to obtain the long-term trend of dust concentration. The long-term trend of dust concentration shows a decreasing trend over the TP, with domain-mean values of −0.012 mg m⁻³month⁻¹ in spring and −0.003 mg m⁻³ month⁻¹ in winter at the 95% significance level. The maximum values of −0.015 mg m⁻³ month⁻¹ in spring and −0.006 mg m⁻³ month⁻¹ in winter are located in the northern TP, which was related to the reduction of dust storms over the TD and Gobi Desert (GD) in recent years (Figures 4C,D).

Snow cover is mainly located in the western TP (Kunlun Mountains) southeast of the TP, especially in the Hengduan Mountains, the Himalayas, and Qilian Mountains, and snow cover is difficult to melt due to the complex terrain. The snow cover in the hinterland of the TP and the Qaidam Basin with a dry climate and scarce precipitation belongs to instantaneous snow, and the snow cover is lower (10% ± 2) (Figure 5C). Both the fraction of coverage and duration of snow cover over the TP are lower than those at high latitudes, which leads to rapid changes in seasonal snow cover over the TP. Due to the influence of local temperature and precipitation, the changes in snow cover are more active in the cold season (Figures 5A,C). The snow cover accumulates gradually during the accumulation period from October to January, which causes the snow cover to have the highest values in January (30% ± 6) and lasts until spring (Figure 5A).

The spatial distributions of surface temperature are opposite to snow cover (Figures 5B,D). The monthly variability of snow on the TP is opposite to surface temperature, with snow reaching minimum values in July and August while surface temperature is at a maximum and snow reaching a maximum in December and January while surface temperature is at a minimum. With increased surface temperature, snow cover will decrease. The surface temperature has larger values in the northern TP, especially the Qaidam Basin, with the maximum value being 11.4°C, and the lowest values mainly occur in the middle, west, and southeast of the TP (−6.3°C ± 0.83, −8.8°C ± 1.15, and −5.0°C ± 1.30, respectively) (Figure 5D). The surface temperature over the TP has the highest value in summer (10.2°C ± 0.69) and the lowest value in winter (−16.7°C ± 1.52), with the highest value in July (11.4°C ± 0.47) and the lowest value in January (−18.6°C ± 1.58) (Figure 5B). The surface temperature has a larger value over the TP in summer; thus, the snow cover of only 15% indicates that most of the snow cover has melted and supplied water to the lake (Zhang et al., 2012), and the minimum snow cover appears in August (12.97% ± 3.92) (Figures 5A,B).

To compare the snow with dust, the snow cover was interpolated to a spatial resolution of 0.5° × 0.625° of MERRA-2. With the lower surface temperature and higher dust concentration over the TP in the winter half-year (December to May), we further investigated the influence of dust aerosols on snow cover over the TP. The distribution of snow cover in spring and winter is consistent with the annual snow cover (Figures 6A,B). In the eastern 90°E of the TP, snow cover increases in spring and decreases in winter. In the western 90°E of the TP, snow cover increases in spring and winter (Figures 6C,D), but the variation trend of dust concentration is opposite to that of snow cover in spring (Figures 4C,D). May dust lead to significant melting of snow in spring and winter, causing most of them at risk of rapid mass loss? We will further discuss the relationship between dust and snow cover. However, due to the accumulation of snow, the relationship between snow cover and dust concentration is not clear (Figures 7C,D).

Snow cover has a larger value (80% ± 12) in the southeastern TP in spring from 2009 to 2018, and the lower values (10% ± 1) are mainly in the center of the TP, including the Kunlun Mountains and Qaidam Basin (Figure 6A). The spatial distribution of snow cover in winter is almost equal to that in spring, while for the entire domain, snow cover in winter is lower than that in spring, and the maximum values (70% ± 8) are mainly located in the southeastern TP. The value of snow cover is 60% ± 4 in Qinghai Lake (Figure 6B). There is an obviously increasing trend of snow cover in the southern TP in spring from 2009 to 2018. The maximum value was mainly located in the southeastern TP, with a maximum value of up to 1.6% month⁻¹; the snow cover in the northeastern TP shows a decreasing trend (−1.2% month⁻¹). However, whether this change is due to dust will be discussed further below. With the deposition of dust aerosols from the TD in the northern TP, dust aerosols mainly influence fresh snow in the northern TP (Figure 6C). The variation trend of snow cover in winter is significantly different from that in spring. Most areas of the TP show decreasing trends and are mainly located in the west, southeast, and northeast of the TP, with a maximum value of up to −2% month⁻¹. Could dust have caused the decrease in snow cover in regions of maximum snow cover? We will use the rate of dust to analyze. The snow cover shows an increasing trend in the southwestern TP, which is mainly caused by the influence of precipitation (Figure 7D).

Due to the influence of precipitation and temperature, the correlation between snow and dust is uncertainty. Therefore, positive correlation will occur in the case of low precipitation, and inverse phase deposition of dust in snow will occur when the dust concentration is large, while the influence of dust in the atmosphere and dust in snow is different, so the correlation results are different. Therefore, we cannot get an obvious signal from the climatic state, so we use the rate of dust change to analyze.

To explain the influence of dust on snow over the TP, this study analyzed the temporal and spatial correlations between dust concentration and snow cover. Dust has a significant correlation with snow cover in most regions of the TP at an 85% significance level of the t test. In particular, dust has a greater influence on snow cover duration in the regions of maximum snow cover (Figure 7). In spring, there is a significant negative correlation between dust and snow cover in the northern 35°N and eastern 90°E of the TP, with the highest correlation coefficients being −0.6. However, in winter, there is a significant positive correlation between dust and snow cover in the TP, with the highest correlation coefficients being 0.5. In the northwest and the hinterland of the TP, there is a significant positive correlation between dust and snow cover, with values of 0.4 in spring and 0.5 in winter, which are mainly influenced by precipitation more than dust aerosols (Figures 7C,D).

In spring, there is a significant negative correlation between dust and snow cover in the northern 35°N and eastern 90°E of the TP. Surface temperature and precipitation are two direct factors that can affect snow, while dust usually affect the snow by affecting surface temperature. The influence of precipitation will lead to a non-significant correlation in some areas. Therefore, we select the northern region of the TP (35–40°N, 87–105°E) and calculate the rate of dust and the variability characteristics of precipitation and select the high and low values of dust rate in the case of precipitation less than 0.0049 mm, so as to reduce the influence of precipitation on snow and then analyze the influence of dust on snow in the northern region of the TP (35–40°N, 87–105°E).

We select the northern region of the TP (35–40°N, 87–105°E) and calculate the rate of dust and the variability characteristics of precipitation and select the high and low values of dust rate in the case of precipitation less than 0.0049 mm, so as to reduce the influence of precipitation on snow and then analyze the influence of dust on snow in the northern region of the TP (35–40°N, 87–105°E). Dust physical processes, including dust emissions, long-term transport, and dry/wet deposition, are closely related to the climatic effects of dust aerosols. The study of the climate effect of dust is always very complex. It is difficult to separately analyze the effect of dust in the atmosphere and snow for observations, which leads to the lack of observational facts in the climate effects of dust aerosols over the TP. We distinguished between atmospheric dust and dust-on-snow by calculating the diurnal rate of dust concentration in the northern region of the TP (35–40°N, 87–105°E) in this study (Figure 8). The highest (red dots, >0.65 mg m⁻³ d⁻¹) and the lowest diurnal rates of atmospheric dust concentration (blue dots, <−0.63 mg m⁻³ d⁻¹) are usually located next to each other in Figure 8. In particular, the dust concentration is higher during high events (red dots) due to local dust emissions and transport from other sources. However, dust aerosols settle quickly to snow due to dry/wet deposition, which leads to dust concentrations in the atmosphere decreasing rapidly in low events (blue dots). The pattern of surface dust concentration in high and low events is almost the opposite in the northern TP (Figures 8, 9). Therefore, the high and low values of the diurnal rates of atmospheric dust represent atmospheric dust and dust-on-snow over the TP, respectively. Furthermore, we investigate the influencing mechanism of atmospheric dust and dust-on-snow on snow cover over the TP based on high and low events.

Based on high and low values of dust variability, using radiation data from MERRA-2 and meteorological element data from ERA5, we analyze the characteristics of changes in solar radiation, surface temperature, and snow cover. Through analyzing the distribution characteristics of radiation, surface temperature, and snow cover in the northern region of the TP (35–40°N, 87–105°E) based on the high and low values of the rate of dust, we can analyze the influence of atmospheric dust and dust-on-snow on snow cover in the northern region of the TP (35–40°N, 87–105°E). The solar radiation anomaly and surface temperature anomaly are mainly obtained by calculating the average values of solar radiation and surface temperature in the high and low values to subtract the average values from 2009 to 2018 in spring. Dust can affect solar radiation directly, while the variation of solar radiation directly affects surface temperature. The radiation anomaly and surface temperature anomaly can be interpreted as being caused by dust. The precipitation is lower in the high and low events; thus, we can ignore the influence of precipitation on snow over the northern TP (Figures 9C,D).

Atmospheric dust reduces the surface net solar radiation (−3.84 W m⁻²) through absorbing shortwave radiation, decreases the surface temperature (−2.27°C), and finally increases the snow cover (−1.04%). The high events occur mainly from March to May, which are mainly from local dust emissions and external transport. During high events, surface dust concentration anomaly has negative values (−0.026 mg m⁻³), and the minimum value is −0.05 mg m⁻³. During this time, surface dust aerosols are emitted to the atmosphere. As one of the major LAAs, dust aerosols can significantly absorb solar radiation. Due to the influence of atmospheric dust, the surface net solar radiation anomaly has lower values of −3.84 W m⁻² in the northern TP (Figure 10C), which caused a negative surface temperature anomaly of −2.27°C in the northern TP during the high events (Figure 11A). The surface temperature can significantly influence snow cover, and the snow cover anomaly shows higher values of 3% in the eastern 90°E of the north TP and lower values of 2% in the western 90°E of the north TP. The spatial distribution of surface albedo is similar to that of snow cover; a positive surface albedo anomaly mainly occurs in the northern TP (0.0004).

The dust deposited on snow reduces the surface albedo by −0.004, induces surface warming (0.42°C), and reduces snow cover (−2.00%) by rapid snowmelt in the northern TP. The positive surface dust concentration anomaly is high over the northern TP (0.057 mg m⁻³), and the negative values of −0.01 mg m⁻³ are located on the southern 37°N and eastern 98°E of the northern TP during low events (Figure 9B). Dust deposited on snow can reduce surface albedo; thus, the surface albedo anomaly has lower values of −0.004 in the northern TP, with minimum values down to −0.016. (Figure 10B). The spatial distribution of the surface net solar radiation anomaly is contrary to that of the surface albedo, and the larger values of 3 W m⁻² are located in the south of northern TP (Figure 10D). The surface temperature has a larger value of 2°C, mainly in the west of the northern TP (Figure 11B). The spatial distribution of snow cover is almost opposite to that of surface temperature, while for the entire snow cover domain, the lower values of snow cover of −2.00% are mainly located north of the TP (Figure 11D). The significant decreasing snow cover over the TP is induced due to the snow darkening effect. Dust-on-snow increases the surface solar shortwave radiation by decreasing the surface albedo, which further strengthens the surface temperature, weakens the snow cover, and further increases dust emissions. Hence, TP dust can create a positive feedback loop, which decreases snow cover. In the same period, snow cover can partially block the upward longwave radiation from the ground and affect the ground thermal regime. When the snow cover decreases, the upward longwave radiation from the ground increases, which weakens the surface temperature. The effect of snow darkening is higher than the thermal regime of snow cover, which finally increases surface temperature and decreases snow cover. Figure 12B shows the significant increase in the surface latent heat flux by the dust-on-snow. This is due to the increased soil moisture induced by a rise in the amount of snowmelt over the northern TP in spring. The warming of the northern TP also increases the regional surface sensible heat flux in Figure 12C. The total surface heat flux (latent + sensible heat flux) shows a larger value over the northern TP (4.29 W m⁻²). Hence, dust-on-snow over the TP can warm the TP and enhance its thermal effects by increasing the surface latent heat flux and sensible heat flux, thereby affecting the Asian climate.