Clear-Sky Surface Solar Radiation and the Radiative E ﬀ ect of Aerosol and Water Vapor Based on Simulations and Satellite Observations over Northern China

: The distribution and trend of clear-sky surface solar radiation (SSR) and the quantitative e ﬀ ects of aerosol and water vapor are investigated in northern China during 2001–2015 using radiation simulations and satellite observations. Clear-sky SSR in northern China is high in summer and low in winter, which is dominated by astronomical factors and strongly modulated by the seasonal variations of radiative e ﬀ ects of aerosol (ARE) and water vapor (WVRE). The larger variation of WVRE than ARE indicates that water vapor plays a more important role in moderating the seasonal variation of clear-sky SSR. Clear-sky SSR shows an overall decreasing trend of –0.12 W / m 2 per year, with decrease more strongly than –0.60 W / m 2 per year in west-central Shandong and increase (about 0.40 W / m 2 ) in south-central Inner Mongolia. The consistency of spatial distribution and high correlation between clear-sky SSR and ARE trend indicate that the clear-sky SSR trend is mainly determined by aerosol variation. Dust mass concentration decreases about 16% in south-central Inner Mongolia from 2001 to 2015, resulting in the increase in clear-sky SSR. In contrast, sulfate aerosol increases about 92% in west-central Shandong, leading to the decreasing trend of clear-sky SSR.


Introduction
Solar radiation reaching the Earth's surface, also known as the surface solar radiation (SSR), is the primary energy source for life on the planet [1,2]. It drives the majority of the terrestrial processes (e.g., the carbon cycle and hydrological cycle) and plays a crucial role in a large number of sectors (e.g., agriculture and solar energy production) [3][4][5]. In the past decade, the role of solar energy in the global energy transformation has increased rapidly. The percentage of solar power in global electricity generation has increased from 0.1% in 2008 to 2.2% in 2018 [6].
SSR is controlled by many factors such as solar geometries and atmospheric constituents [7,8]. Solar geometries including the Earth-Sun distance and the solar zenith angle determine solar radiation reaching the top of the atmosphere and drive the seasonal and diurnal variations of SSR [9,10]. Solar radiation is modulated by atmospheric constituents as it passes through the atmosphere. Cloud has been regarded as a key modulator of solar radiation in the atmosphere. Previous studies have found that the reduction in solar radiation from the 1960s to the 1980s at many sites coincides with the increase of cloud cover [11,12]. However, cloud change cannot always explain the variation a combined dataset, combining the DT and DB data. The resultant merged AOD field is more spatially complete than either the DB or DT AOD fields [40], which provides extended spatial coverage and is important for investigating aerosol loading over areas with complex topography and underlying surface in North China. We use the combined Dark Target and Deep Blue product from Terra MODIS Level-3 monthly data in 1 • × 1 • resolution. Ruiz-Arias (2013) assessed the Level-3 MODIS AOD against more than 500 AERONET ground stations around the globe and the root mean square error is 0.14 [41].
Single scattering albedo (SSA) used for solar radiation modeling come from the Kinne aerosol climatology, which is developed based on merging of data from Aerosol Comparisons between Observations and Models (AEROCOM) global modeling and ground measurements from the Aerosol Robotic Network (AERONET) [42][43][44]. The Kinne aerosol climatology contains monthly SSA data in 1 • × 1 • resolution with RMSE of 0.06 in comparison with SSA from AERONET [42].
The water vapor data is obtained from the ERA-Interim data, which is a global atmospheric reanalysis data provided by the European Centre for Medium-Range Weather Forecasting (ECMWF). The data provides meteorological parameters such as surface temperature, soil temperature, and relative humidity [45]. The ERA-Interim data used here is the monthly water vapor with a spatial resolution of 0.25 • × 0.25 • . The RMSE of water vapor is 2.8 kg/m 2 in comparison with radiosonde observations [46].
The solar radiation data in ground observation stations is from the daily data set of basic elements of meteorological radiation in China provided by the National Meteorological Information Center. SSR from 10 stations during 2001-2015 is used to validate the modeled SSR in this study, and the locations of these stations are shown in Figure 1.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 21 combined dataset, combining the DT and DB data. The resultant merged AOD field is more spatially complete than either the DB or DT AOD fields [40], which provides extended spatial coverage and is important for investigating aerosol loading over areas with complex topography and underlying surface in North China. We use the combined Dark Target and Deep Blue product from Terra MODIS Level-3 monthly data in 1°×1° resolution. Ruiz-Arias (2013) assessed the Level-3 MODIS AOD against more than 500 AERONET ground stations around the globe and the root mean square error is 0.14 [41]. Single scattering albedo (SSA) used for solar radiation modeling come from the Kinne aerosol climatology, which is developed based on merging of data from Aerosol Comparisons between Observations and Models (AEROCOM) global modeling and ground measurements from the Aerosol Robotic Network (AERONET) [42][43][44]. The Kinne aerosol climatology contains monthly SSA data in 1°×1° resolution with RMSE of 0.06 in comparison with SSA from AERONET [42].
The water vapor data is obtained from the ERA-Interim data, which is a global atmospheric reanalysis data provided by the European Centre for Medium-Range Weather Forecasting (ECMWF). The data provides meteorological parameters such as surface temperature, soil temperature, and relative humidity [45]. The ERA-Interim data used here is the monthly water vapor with a spatial resolution of 0.25°×0.25°. The RMSE of water vapor is 2.8 kg/m 2 in comparison with radiosonde observations [46].
The solar radiation data in ground observation stations is from the daily data set of basic elements of meteorological radiation in China provided by the National Meteorological Information Center. SSR from 10 stations during 2001-2015 is used to validate the modeled SSR in this study, and the locations of these stations are shown in Figure 1. Aerosol data from the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) is also used. MERRA-2 product is an atmospheric reanalysis provided by the NASA Global Modeling and Assimilation Office (GMAO) [47,48]. In this study, the surface aerosol mass concentration data with a spatial resolution of 0.625°× 0.5° from the "tavgM_2d_aer_Nx" dataset is used to analyze the influence of aerosol compositions on clear-sky SSR variation. Aerosol data from the second Modern-Era Retrospective analysis for Research and Applications (MERRA-2) is also used. MERRA-2 product is an atmospheric reanalysis provided by the NASA Global Modeling and Assimilation Office (GMAO) [47,48]. In this study, the surface aerosol mass concentration data with a spatial resolution of 0.625 • × 0.5 • from the "tavgM_2d_aer_Nx" dataset is used to analyze the influence of aerosol compositions on clear-sky SSR variation. In this study, the clear-sky SSR is simulated using MAGIC [49], which is based on radiation simulated by a Radiative Transfer Model (RTM) and stored in look-up tables to improve the computing efficiency. It assumes plane atmosphere and uses the modified Lambert function for each wavelengths under cloud-free conditions [49].The RTM model libRadtran [50] has been used for the extensive analysis of the interaction between radiation and the atmosphere and the computation of the LUTs. All RTM calculations for LUTs establishment are performed for the US standard atmosphere, which provides vertical information on the total number density of atmospheric constituents. The clear-sky SSR is simulated with input of aerosol properties from MODIS product and Kinne aerosol climatology, water vapor from ERA-Interim water vapor and surface albedo from the Surface and Radiation Budget/Clouds and the Earth's Radiant Energy System (SARB/CERES) [51]. The clear-sky SSR assuming conditions without aerosol and without water vapor is then simulated by adjusting the extinction of aerosol and water vapor to zero, respectively, remaining other inputs unchanged. All the simulations are conducted in time interval of 15 min and resolution of 0.5 • × 0.5 • . The instantaneous ARE (ARE int ) and WVRE (WVRE int ) is calculated as follows: where F clear is the clear-sky solar flux including aerosol and water vapor, F no−aerosol is the clear-sky solar flux without aerosol and with water vapor, and F no−wv is the clear-sky solar flux with aerosol and without water vapor. Then, the daily ARE (ARE daily ) and WVRE (WVRE daily ) are calculated by the integration of the instantaneous data over 24 h: The trends of clear-sky SSR, ARE, and WVRE are calculated based on their monthly mean values. In order to eliminate the great influence of seasonal cycle, the data is first deseasonalized by subtracting the seasonal cycle obtained from the harmonic regression [52,53]. The trends (ω) of clear-sky SSR, ARE, and WVRE are then calculated by applying the least-square fit to the deseasonalized time series. The significance of trend is evaluated using the ratio of |ω/σ ω |, where σ ω is the standard deviation of the trend [54]. The trend is regarded as significant at the 95% confidential level when |ω/σ ω | is larger than 2.

Validation of Simulated Clear-Sky SSR
The reliability of simulated clear-sky SSR is validated by comparison with the ground observations. As SSR in stations are observed under all-sky conditions, SSR observations should be distinguished between cloudy and clear-sky conditions. An empirical threshold optimization method is used to select clear-sky conditions from the SSR observations. If the observed daily SSR value exceeds 89% of the clear-sky solar radiation simulated by the MAGIC, a given day can be identified as a clear-sky condition [55]. The threshold for clear-sky situations refers to Bartok (2017) [55], who selected the empirical optimum threshold by testing different thresholds against the clear-sky detection and interpolation algorithm of Long and Ackerman (2000) [56] using the Baseline Surface Radiation Network (BSRN) data. As the consequence, the optimal threshold showing the higher average between clear-sky and all-sky matches is detected to be 89%. Table 1 and Figure 2 shows the statistical analysis between observed and simulated clear-sky SSR at 10 stations. The overall slope of fitted straight line between observations and simulations at 10 stations is 1.03, the determination coefficient (R 2 ) is 0.98, the mean absolute error (MAE) is 8.13 W/m 2 (3.62%), and the root mean square error (RMSE) is 10.60 W/m 2 (4.17%). The accuracy of simulation from this study is better than results from Otunla (2019) [57] (RMSE: 15-26%) and Zhang and Wen (2014) [58] (RMSE: 9.6%). R 2 in all the stations are higher than 0.96, and MAE range from 6.17-12.98 W/m 2 . RMSE at nine stations are less than 7% in the Tianjin station. From above, the comparison shows a very good agreement, which confirms the reliability of the simulation in this study.

Spatial Pattern of Clear-Sky SSR and Associated Factors in North China
As shown in Figure 3, the 15-year average of clear-sky SSR in North China presents a clear distribution of high values in the northwest and low values in the southeast. A low-value center with clear-sky SSR less than 215 W/m 2 is found in the neighboring region of Shandong, Hebei, and Tianjin, which corresponds to core of the North China Plain (NCP). To the southeast of the low-value center, the clear-sky SSR increase slightly to about 225 W/m 2 in south and east of Shandong. And to the northwest of the low-value center, the clear-sky SSR increase dramatically to more than 240 W/m 2 in Shanxi, northwest Hebei and south-central Inner Mongolia. Yu et al. (2020) studied the clear-sky solar radiation over arid and semi-arid areas in china, and their results shows that the clear-sky SSR over the south-central Inner Mongolia is about 250 W/m 2 , which is consistent to our results.

Spatial Pattern of Clear-Sky SSR and Associated Factors in North China
As shown in Figure   Comparing spatial distribution of clear-sky SSR in Figure 3 and altitudes in Figure 1, the locations with low clear-sky SSR correspond to those with low altitudes. The low-altitude plain regions are usually characterized with dense population and developed industry, and a large amount of anthropogenic aerosol is emitted by industrial and domestic activities [59]. Comparing spatial pattern of clear-sky SSR in Figure 3 and ARE in Figure 4a, regions with low clear-sky SSR correspond to regions with strong attenuation effect of aerosol on solar radiation. High ARE (larger than -40 Comparing spatial distribution of clear-sky SSR in Figure 3 and altitudes in Figure 1, the locations with low clear-sky SSR correspond to those with low altitudes. The low-altitude plain regions are usually characterized with dense population and developed industry, and a large amount of anthropogenic aerosol is emitted by industrial and domestic activities [59]. Comparing spatial pattern of clear-sky SSR in Figure 3 and ARE in Figure 4a, regions with low clear-sky SSR correspond to regions with strong attenuation effect of aerosol on solar radiation. High ARE (larger than −40 W/m 2 ) are found in the neighboring region of Shandong, Hebei and Tianjin, which are much higher than the national mean ARE in China (−15.70 W/m 2 ) [60]. However, low ARE (smaller than −15 W/m 2 ) are found in north Hebei and south-central Inner Mongolia. High terrain altitude induces low total atmospheric column and large clear-sky SSR, also contributing to the spatial distribution of SSR [61].
The WVRE (Figure 4b) shows an obvious distribution of high values in the southeast and low values in the northwest as the distribution of water vapor is influenced by the distance to the sea, underlying surface, latitude, and topography [34,62]. The strongest attenuation effect of water vapor on solar radiation (larger than −40 W/m 2 ) occur in north Henan and south Shandong. While the weakest (around −30 W/m 2 ) occur in south-central Inner Mongolia due to high altitude and latitude, dry climate and long distance to the sea [63]. The spatial difference of extinction of water vapor on solar radiation contributes to the spatial distribution of clear-sky SSR. However, the spatial difference of WVRE is much smaller than ARE, and the spatial distribution of clear-sky SSR is more strongly influenced by aerosol.
W/m 2 ) are found in the neighboring region of Shandong, Hebei and Tianjin, which are much higher than the national mean ARE in China (-15.70 W/m 2 ) [60]. However, low ARE (smaller than -15 W/m 2 ) are found in north Hebei and south-central Inner Mongolia. High terrain altitude induces low total atmospheric column and large clear-sky SSR, also contributing to the spatial distribution of SSR [61]. The WVRE (Figure 4b) shows an obvious distribution of high values in the southeast and low values in the northwest as the distribution of water vapor is influenced by the distance to the sea, underlying surface, latitude, and topography [34,62]. The strongest attenuation effect of water vapor on solar radiation (larger than -40 W/m 2 ) occur in north Henan and south Shandong. While the weakest (around -30 W/m 2 ) occur in south-central Inner Mongolia due to high altitude and latitude, dry climate and long distance to the sea [63]. The spatial difference of extinction of water vapor on solar radiation contributes to the spatial distribution of clear-sky SSR. However, the spatial difference of WVRE is much smaller than ARE, and the spatial distribution of clear-sky SSR is more strongly influenced by aerosol. . The seasonal variation of solar radiation at a given location is mainly driven by the astronomical factors such as the tilt of the Earth's rotation axis with respect to the Earth's orbital plane around the Sun and the resulting seasonality of insolation and the Earth-Sun distance. Figure 6 shows the monthly variation of regional average of clear-sky solar radiation at the surface and the top of the atmosphere (TOA) in northern China. The regional average of clear-sky SSR increases from 127.58 W/m 2 in January to the peak of 323.24 W/m 2 in May and then decreases to a low level in December, showing a single-peak change. The seasonal pattern of clear-sky SSR is similar to that at the TOA, indicating that the seasonal variation of clear-sky SSR is primarily determined by the solar radiation reaching the TOA, which is modulated by the astronomical factors.  . The seasonal variation of solar radiation at a given location is mainly driven by the astronomical factors such as the tilt of the Earth's rotation axis with respect to the Earth's orbital plane around the Sun and the resulting seasonality of insolation and the Earth-Sun distance. Figure 6 shows the monthly variation of regional average of clear-sky solar radiation at the surface and the top of the atmosphere (TOA) in northern China. The regional average of clear-sky SSR increases from 127.58 W/m 2 in January to the peak of 323.24 W/m 2 in May and then decreases to a low level in December, showing a single-peak change. The seasonal pattern of clear-sky SSR is similar to that at the TOA, indicating that the seasonal variation of clear-sky SSR is primarily determined by the solar radiation reaching the TOA, which is modulated by the astronomical factors.     In summer, dramatic increases of clear-sky SSR (about 50 W/m 2 ) is found in the northwest part of northern China, which is due to the increased solar radiation reaching the TOA as a result of short Earth-Sun distance in summer [10]. However, clear-sky SSR show little increase in the boundary of Shandong, Hebei, and Henan. In these regions, ARE show a significant increase to more than −60 W/m 2 , indicating a much stronger attenuation effect of aerosol on solar radiation in summer. WVRE also shows a large enhancement from spring to summer. In autumn, a large decrease in clear-sky SSR is found and the differences among regions are reduced in northern China. The decrease is mainly dominated by the astronomical factors as the radiative effects of aerosol and water vapor on solar radiation are weakened simultaneously. In winter, clear-sky SSR decrease to lower than 170 W/m 2 in most regions of northern China and the lowest values (about −130 W/m 2 ) are found in the boundary of Tianjin, Shandong, and Hebei. The attenuation effect of aerosol and water vapor on solar radiation also reach the valley in winter. ARE values are even lower than −7 W/m 2 in parts of Hebei and Inner Mongolia, and WVRE values are lower than −20 W/m 2 in most of northern China.

Seasonal Variation of Clear-Sky SSR and Driving Factors in North China
The seasonal variation of clear-sky SSR is dominated by the solar radiation reaching the TOA, as discussed above. However, it should be noted that there is some difference between the monthly pattern of solar radiation at the surface and that at the TOA. The peak at the surface shows in May, and solar radiation are strongly attenuated in June, July, and August. Figure 9 shows the monthly variation of ARE and WVRE. ARE increases from January as temperature increases, which is mainly due to the increase of gas-particle transformation and hygroscopic growth of aerosol particles with temperature and humidity [64,65]. In spring, sandstorms occur frequently in northwest China, and dust aerosol is transported to northern China, which increases aerosol loading and its extinction on solar radiation [33]. Dust aerosol is characterized with high absorption in low wavelengths and results in low SSA in this season, which also contributes the increase in ARE. Pani et al. (2018) studied the ARE over an urban atmosphere in northern peninsular southeast Asia (PSEA) in 2014, and showed ARE ranged between −45.3 to −103.4 W/m 2 from March to April [66], which is larger than that (lower than −30 W/m 2 ) in our study during the same period. The difference mainly come from different aerosol loading and properties. Influenced by increased aerosol emission from biomass burning from February to April, AOD was larger than 1 and SSA is lower than 0.9 in some regions of PSEA, which contributed to the higher ARE [66,67]. The maximum ARE (−38.90 W/m 2 ) occurs in July, and then ARE decreases sharply in August and September. The decrease of AOD due to the wet deposition by heavy rain causes a decrease in ARE in August and September [68]. Moreover, large SSA values (0.94) are found in these months, and the decrease of aerosol absorption also contributes to the decrease in ARE. The minimum ARE (−11.40 W/m 2 ) occurs in December, when the Asian winter monsoon brings clean and dry air to northern China.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 21 Figure 7 and Figure 8 show the spatial distribution of seasonal mean ARE and WVRE in northern China. ARE and WVRE in spring are lower than -20 W/m 2 and -40 W/m 2 in the northwest regions of northern China, respectively, and about -45 W/m 2 and -50 W/m 2 in west Shandong, southeast Hebei, and parts of north Henan. In summer, dramatic increases of clear-sky SSR (about 50 W/m 2 ) is found in the northwest part of northern China, which is due to the increased solar radiation reaching the TOA as a result of short Earth-Sun distance in summer [10]. However, clear-sky SSR show little increase in the boundary of Shandong, Hebei, and Henan. In these regions, ARE show a significant increase to more than -60 W/m 2 , indicating a much stronger attenuation effect of aerosol on solar radiation in summer. WVRE also shows a large enhancement from spring to summer. In autumn, a large decrease in clear-sky SSR is found and the differences among regions are reduced in northern China. The decrease is mainly dominated by the astronomical factors as the radiative effects of aerosol and water vapor on solar radiation are weakened simultaneously. In winter, clear-sky SSR decrease to lower than 170 W/m   The seasonal variation of clear-sky SSR is dominated by the solar radiation reaching the TOA, as discussed above. However, it should be noted that there is some difference between the monthly pattern of solar radiation at the surface and that at the TOA. The peak at the surface shows in May, and solar radiation are strongly attenuated in June, July, and August. Figure 9 shows the monthly variation of ARE and WVRE. ARE increases from January as temperature increases, which is mainly due to the increase of gas-particle transformation and hygroscopic growth of aerosol particles with temperature and humidity [64,65]. In spring, sandstorms occur frequently in northwest China, and dust aerosol is transported to northern China, which increases aerosol loading and its extinction on WVRE are low in winter and even lower than ARE in January and February. However, the extinction of solar radiation by water vapor increases dramatically from March to July as the Asian summer monsoon brings warm and humid air [69,70]. As shown in Figure 10, water vapor increases simultaneously. WVRE reaches the maximum of -70.88 W/m 2 in July, which is about 1.82 times that of ARE. In August and September, though a large decrease is found in WVRE, the decrease in ARE is more obvious, and the ratios of WVRE to ARE increase to 2.14 and 2.24, respectively. The high values of WVRE and ARE in June, July, and August are responsible for the inconsistent seasonal variation of solar radiation at the surface and TOA. In the following months, WVRE decreases dramatically and reaches to the minimum of -13.84 W/m 2 in December as a result of the dry air condition caused by the winter monsoon [69,70]. As shown in Figure 9, the monthly variation of WVRE is much more drastic than ARE, and water vapor plays a more important role in moderating the seasonal variation of clear-sky SSR. WVRE are low in winter and even lower than ARE in January and February. However, the extinction of solar radiation by water vapor increases dramatically from March to July as the Asian summer monsoon brings warm and humid air [69,70]. As shown in Figure 10, water vapor increases simultaneously. WVRE reaches the maximum of -70.88 W/m 2 in July, which is about 1.82 times that of ARE. In August and September, though a large decrease is found in WVRE, the decrease in ARE is more obvious, and the ratios of WVRE to ARE increase to 2.14 and 2.24, respectively. The high values of WVRE and ARE in June, July, and August are responsible for the inconsistent seasonal variation of solar radiation at the surface and TOA. In the following months, WVRE decreases dramatically and reaches to the minimum of −13.84 W/m 2 in December as a result of the dry air condition caused by the winter monsoon [69,70]. As shown in Figure 9, the monthly variation of WVRE is much more drastic than ARE, and water vapor plays a more important role in moderating the seasonal variation of clear-sky SSR. WVRE are low in winter and even lower than ARE in January and February. However, the extinction of solar radiation by water vapor increases dramatically from March to July as the Asian summer monsoon brings warm and humid air [69,70]. As shown in Figure 10, water vapor increases simultaneously. WVRE reaches the maximum of -70.88 W/m 2 in July, which is about 1.82 times that of ARE. In August and September, though a large decrease is found in WVRE, the decrease in ARE is more obvious, and the ratios of WVRE to ARE increase to 2.14 and 2.24, respectively. The high values of WVRE and ARE in June, July, and August are responsible for the inconsistent seasonal variation of solar radiation at the surface and TOA. In the following months, WVRE decreases dramatically and reaches to the minimum of -13.84 W/m 2 in December as a result of the dry air condition caused by the winter monsoon [69,70]. As shown in Figure 9, the monthly variation of WVRE is much more drastic than ARE, and water vapor plays a more important role in moderating the seasonal variation of clear-sky SSR.  Figure 11 shows the monthly clear-sky SSR time series and the deseasonlized clear-sky SSR trend in northern China from 2001 to 2015. The clear-sky SSR in North China shows an overall decreasing trend of -0.12 W/m 2 per year but is not significant at the 95% confidential level. Figure 12 shows the spatial distribution of clear-sky SSR trends in northern China from 2001 to 2015. The trends among different regions vary greatly. There is on significant trend over Beijing, north Hebei, and south Shanxi. While positive trend with annual growth rate larger than 0.30 W/m 2 are found in south-central Inner Mongolia. On the other hand, a significant decreasing trend (larger than −0.60 W/m 2 per year) occurs in west-central Shandong, in some places of which trend up to -0.80 W/m 2 per year is even found. Clear-sky SSR in these regions decrease about 4% in 15 years.  Figure 11 shows the monthly clear-sky SSR time series and the deseasonlized clear-sky SSR trend in northern China from 2001 to 2015. The clear-sky SSR in North China shows an overall decreasing trend of -0.12 W/m 2 per year but is not significant at the 95% confidential level. Figure 12 shows the spatial distribution of clear-sky SSR trends in northern China from 2001 to 2015. The trends among different regions vary greatly. There is on significant trend over Beijing, north Hebei, and south Shanxi. While positive trend with annual growth rate larger than 0.30 W/m 2 are found in southcentral Inner Mongolia. On the other hand, a significant decreasing trend (larger than -0.60 W/m 2 per year) occurs in west-central Shandong, in some places of which trend up to -0.80 W/m 2 per year is even found. Clear-sky SSR in these regions decrease about 4% in 15 years.        Comparing clear-sky SSR trend in Figure 12 and ARE trend in Figure 13a, the areas where clearsky SSR decrease/increase have a good match with the areas where ARE is enhanced/weakened. The areas with the largest decrease of clear-sky SSR and the largest enhancement of ARE are located in west-central Shandong, while the areas with the largest increase of clear-sky SSR, and the largest weakening of ARE is found in south-central Inner Mongolia. The good consistency in the spatial distribution of clear-sky SSR and ARE trend suggests that the variation of ARE may be the dominated factor for the long-term variation of clear-sky SSR. Comparing clear-sky SSR trend in Figure 12 and ARE trend in Figure 13a, the areas where clear-sky SSR decrease/increase have a good match with the areas where ARE is enhanced/weakened. The areas with the largest decrease of clear-sky SSR and the largest enhancement of ARE are located in west-central Shandong, while the areas with the largest increase of clear-sky SSR, and the largest weakening of ARE is found in south-central Inner Mongolia. The good consistency in the spatial distribution of clear-sky SSR and ARE trend suggests that the variation of ARE may be the dominated factor for the long-term variation of clear-sky SSR. Figure 14a shows the comparison of clear-sky SSR trend and co-located ARE trend. The slope of linear regression is about 0.981 and R 2 between clear-sky SSR and ARE trend is 0.995, which indicates that the clear-sky SSR trend is mainly determined by the variation of attenuation effect on solar radiation by aerosol.

Long-Term Trend of Clear-Sky SSR and the Roles of Aerosol and Water Vapor in Northern China
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 21 Figure 14a shows the comparison of clear-sky SSR trend and co-located ARE trend. The slope of linear regression is about 0.981 and R 2 between clear-sky SSR and ARE trend is 0.995, which indicates that the clear-sky SSR trend is mainly determined by the variation of attenuation effect on solar radiation by aerosol.   Figure 14b shows the comparison between clearsky SSR trend and WVRE trend. The trend of WVRE varies slightly as trend of clear-sky SSR changes greatly, and R 2 between them is lower than 0.07, indicating that the variation of extinction by water vapor has little impact on the trend of clear-sky SSR. Figure 15 shows    Figure 14b shows the comparison between clear-sky SSR trend and WVRE trend. The trend of WVRE varies slightly as trend of clear-sky SSR changes greatly, and R 2 between them is lower than 0.07, indicating that the variation of extinction by water vapor has little impact on the trend of clear-sky SSR. Figure 15 shows        To figure out the roles of natural and anthropogenic aerosol emissions in clear-sky SSR variation, the annual variation of typical aerosol compositions such as dust and sulfate from MERRA2 is investigated. As shown in Figure 17, dust mass concentration at the surface decreases dramatically in North China, indicating a decreasing influence of dust aerosol. Dust mass concentration in south-central Inner Mongolia decreases about 16% from 2001 to 2015. As a result, the extinction of dust particles on solar radiation decreases and more solar radiation reaches the surface, causing the increase of clear-sky SSR in this region. In contrast, sulfate mass concentration in this region is in low level and shows little variation due to the undeveloped industry and sparse population here. Dust mass concentration in west and central Shandong also shows a great decrease. However, sulfate mass concentration increases about 92% from 2001 to 2015. The dramatic increase in sulfate aerosol leads to the decreasing trend of clear-sky SSR in this area.

Discussion
To assess the influence of input parameters on SSR simulation, we conducted sensitivity tests for SSR simulation by adding a perturbed value for each parameter inputted in the MAGIC model. The sensitivity test is based on two steps-first, SSR was simulated using the measured values as input parameters, obtaining the non-perturbed results. Second, we simulated SSR by adding a perturbed value for an input parameter and keeping other parameters unchanged, thereby obtaining a perturbed result. This step is repeated to obtain perturbed results for all the input aerosol parameters, such as AOD, SSA, and water vapor. A similar strategy has been adopted for uncertainty estimation in previous studies [71]. The perturbed values of AOD, SSA, and water vapor values were selected based on the uncertainties of these parameters as shown in Section 2. AOD was perturbed by ±0.14, SSA was perturbed by ±0.06, and water vapor was perturbed by ±2.80 kg/m 2 . Table 2 shows the absolute and percentage changes between the non-perturbed and perturbed results. An increase of 0.14 in AOD from the measurements causes an SSR decrease of 3.71%. By contrast, a decrease of 0.14 in AOD results in a decrease with similar magnitudes. The largest SSR change in this analysis is associated with SSA uncertainty. An increase of 0.06 in SSA provides an increase of 4.54% for SSR. For water vapor, an increase of 2.80 kg/m 2 in water vapor causes a 0.64% decrease in SSR.

Conclusions
In this study, we analyzed the distribution and trend of clear-sky SSR and the quantitative effects of aerosol and water vapor over northern China based on clear-sky SSR, ARE, and WVRE simulated by the MAGIC. The comparison between simulated and observed clear-sky SSR shows a R 2 of 0.98,

Discussion
To assess the influence of input parameters on SSR simulation, we conducted sensitivity tests for SSR simulation by adding a perturbed value for each parameter inputted in the MAGIC model. The sensitivity test is based on two steps-first, SSR was simulated using the measured values as input parameters, obtaining the non-perturbed results. Second, we simulated SSR by adding a perturbed value for an input parameter and keeping other parameters unchanged, thereby obtaining a perturbed result. This step is repeated to obtain perturbed results for all the input aerosol parameters, such as AOD, SSA, and water vapor. A similar strategy has been adopted for uncertainty estimation in previous studies [71]. The perturbed values of AOD, SSA, and water vapor values were selected based on the uncertainties of these parameters as shown in Section 2. AOD was perturbed by ±0.14, SSA was perturbed by ±0.06, and water vapor was perturbed by ±2.80 kg/m 2 . Table 2 shows the absolute and percentage changes between the non-perturbed and perturbed results. An increase of 0.14 in AOD from the measurements causes an SSR decrease of 3.71%. By contrast, a decrease of 0.14 in AOD results in a decrease with similar magnitudes. The largest SSR change in this analysis is associated with SSA uncertainty. An increase of 0.06 in SSA provides an increase of 4.54% for SSR. For water vapor, an increase of 2.80 kg/m 2 in water vapor causes a 0.64% decrease in SSR.

Conclusions
In this study, we analyzed the distribution and trend of clear-sky SSR and the quantitative effects of aerosol and water vapor over northern China based on clear-sky SSR, ARE, and WVRE simulated by the MAGIC. The comparison between simulated and observed clear-sky SSR shows a R 2 of 0.98, MAE of 8.14 W/m 2 (3.62%), and RMSE of 10.60 W/m 2 (4.17%), confirming the reliability of the simulation in this study.
During 2001-2015, clear-sky SSR shows a clear distribution of high values in the northwest and low values in the southeast. A low-value center with clear-sky SSR less than 215 W/m 2 is found in the boundaries of Shandong, Hebei, and Tianjin, which is characterized with low altitude plains, developed industry, and high anthropogenic emissions. Regions with low clear-sky SSR correspond to areas with strong extinction effect of aerosol on solar radiation, indicating that aerosol loading plays an important role in the spatial distribution of clear-sky SSR.
Clear-sky SSR in northern China shows a similar seasonal pattern with solar radiation at the TOA with high value (304.90 W/m 2 ) in summer and low value (137.62 W/m 2 ) in winter, which is mainly modulated by astronomical factors. The seasonal variation of WVRE is much more drastic than ARE, and water vapor plays a more important role in moderating the seasonal variation of clear-sky SSR.
The Using fixed atmospheric vertical distribution is one source of error in this study, and we will carry out studies adopting observations of atmospheric vertical profiles in SSR simulations to reduce error in the future work.