Long-Term Variations of Global Solar Radiation and Atmospheric Constituents at Sodankylä in the Arctic

: An empirical model of global solar irradiance (EMGSI) under all sky conditions was developed by using solar radiation and meteorological parameters at Sodankylä. The calculated hourly global solar irradiance is in agreement with that observed at the ground during 2008–2011 and at the top of the atmosphere (TOA). This model is used to calculate the global solar irradiance at the ground and its attenuation in the atmosphere due to absorbing and scattering substances in 2000–2018. The sensitivity test indicates that the responses of global solar irradiance to changes in water vapor and scattering factors are nonlinear and negative, and global solar irradiance is more sensitive to changes in scattering (expressed by the scattering factor S/G, S and G are diffuse and global solar radiation, respectively) than to changes in water vapor. Using this empirical model, we calculated the albedos at the TOA and the surface, which are in agreement with the satellite-retrieved values. A good relationship between S/G and aerosol optical depth (AOD) was determined and used to estimate AOD in 2000–2018. An empirical model for estimation of tropospheric NO 2 vertical column density (VCD) was also developed and used to calculate tropospheric NO 2 VCD in 2000–2018. During 2000–2018, the estimated global solar irradiance decreased by 0.92%, and diffuse irradiance increased by 1.28% per year, which is ascribed to the increases of S/G (1.73%) and water vapor (0.43%). Annual surface air temperature increases by 0.07 ◦ C per year. Annual mean loss of global solar irradiance caused by absorbing and scattering substances and total loss are 1.94, 1.17 and 3.11 MJ m − 2 , respectively. Annual mean losses of absorbing and scattering global solar irradiance show negative and positive trends, respectively, and the annual total loss increases by 0.24% per year. Annual mean losses due to absorption were much larger than those due to scattering. The calculated albedos at the TOA are smaller than at the surface. The calculated and satellite-retrieved annual albedos decrease at the TOA and increase at the surface. During 2000–2018, annual means of the AOD and the tropospheric NO 2 VCD increased by 8.23% and 0.03% per year, respectively.


Introduction
The surface temperatures in the Arctic increase much faster than the global average, a phenomenon known as Arctic Amplification. The most significant increase occurs in the winter [1][2][3]. The Arctic warming is related to the imbalance of energy budget at the ground and the top of the atmosphere (TOA), sea-ice-albedo feedback, lapse-rate feedback, temperature inversion, enhanced greenhouse effect [1] and the reasons for the climate change in the Arctic region are still unclear.
The Sun provides important energy to the atmospheric gases, liquids and particles (GLPs) and initiates their changes in the three phases through chemical and photochemical reactions (CPRs). Therefore, solar radiation transfers in and interacts with the atmospheric GLPs. The attenuation in the atmosphere and reflections at the ground and the TOA are necessary to be investigated; these physical and chemical processes control/interact (VOCs) and OH radicals. (2) Scattering, total scattering of global solar radiation by GLPs, is expressed as e −S/G . Then, an empirical model of global solar irradiance at a horizontal plane was applied in the Arctic for all sky conditions [12]: Gcal = A 1 e −kWm × cos(Z) + A 2 e −S/G obs + A 0 (1) where G and S were hourly global and diffuse horizontal solar irradiance (MJ m −2 ), respectively. Coefficients A 1 and A 2 express the individual absorption and scattering processes at the TOA, and A 0 the reflection of global solar irradiation at the TOA. A 1 , A 2 and A 0 represent the averaged values at a mean distance from Earth to Sun in the studying time period, i.e., four years (2008)(2009)(2010)(2011) in this study. The sum of A 1 , A 2 and |A 0 | is the total solar irradiation at the TOA that equals to solar constant (I o = 1367 W m −2 , using the sum of A 1 , A 2 and |A 0 | multiplied by 1,000,000/3600). High-quality data during January 2008 to December 2011 were used to develop a representative empirical model that represents the realistic relationship between physical and chemical processes in the atmosphere: measurements in the early morning and late evening were excluded to avoid the influence of large cosine errors of radiation sensors at high Z (i.e., data obtained for Z <75 o , corresponding to radiation data from 6:00 to 14:00 local time were used). As a result, 3962 hourly data (sample number n = 3962) from April to August in 2008-2011 were used for model development, and the usage rate of the observed data was 53.24%. Hourly mean solar irradiance and meteorological parameters (e.g., T, RH and E) were calculated for each hour, from 6:00 to 14:00, and these hourly averages were calculated for daily and monthly averages.
The coefficients (A 1 , A 2 , A 0 and MJ m −2 ) in Equation (1) were determined by using a multi-parameter fit to the measured hourly mean global solar irradiance. The results are shown in Table 1, i.e., the coefficients, the coefficient of determination (R 2 ), the mean absolute value of relative error (δ) between calculated and measured G, the mean absolute deviation (MAD, in exposure unit, MJ m −2 and in percentage of mean measured value, %), and the root mean square error (RMSE, in exposure unit and in percentage of mean measured value). The scatter plot of calculated vs. observed G is shown in Figure 1. The slope of the linear regression of G cal on G obs (forced through 0) is 0.98 with R 2 of 0.81, which is different from the R 2 in Table 1, because different equations were used (linear regression and Equation (1)). Generally, the estimated G is in agreement with the measured G under all sky conditions. The albedo (A 0 /(A 1 + A 2 + |A 0 |)) at the TOA was 24.1% in 2008-2011. Table 1. The coefficients and constants, coefficient of determination (R 2 ), average and maximum of the absolute relative bias (δavg and δmax (%)), normalized mean square error (NMSE), and standard deviations of calculated and observed solar global irradiances (σ cal and σ obs , respectively). The mean bias error (MAD, MJ m −2 and %) and the root mean square error (RMSE, MJ m −2 and %) (n = 3962).   There was a strong correlation between G and absorbing and scattering terms (R = 0.91 at the confidence level α = 0.001). A stronger correlation also existed between G and the scattering term e −S/G (R = 0.84) than G and absorption term e −kWm ×cos(Z) (R = 0.44), and a weak correlation between absorption term and the scattering term (R = 0.10), indicating that the absorbing and scattering processes can be described separately in general. The RMSE (0.22, Table 1) was in good agreement with the mean RMSE of 0.22 obtained by using 7 empirical models with better simulations out of 105 models under all sky conditions [17], displaying reliable simulations by using the empirical model, and larger than that from using empirical models under clear sky conditions [18], indicating the model simulations are influenced by clouds and aerosols [12][13][14][15][16][17][18].
Measurements of global solar irradiance from 4:00 to 16:00 during 1 April 2008 to 31 August 2011 were used to validate the empirical model of global solar irradiance. The mean absolute relative bias was 24.56%, and normalized mean square error (NMSE) was 0.04. RMSE was 0.20 (MJ m −2 ) and 17.05% (n = 7442). It is reported that the uncertainties for numerous solar radiation models are <20% [26]. The scatter plot of calculated vs. observed global irradiance is given in Figure 2. The estimated monthly average of G was also in agreement with the observed with the relative bias by 4.25% for the average and 9.57% for the maximum, and the RMSE values were 0.05 (MJ m −2 ) and 4.11% (Figure 3). The annual averages of estimated and the observed G showed a similar variation with the relative bias by 3.6% for the average and 5.5% for the maximum (Figure 4). Both the estimated and measured G showed increasing trends by 0.1% and 1.6% per year, respectively during 2008-2011.   According to above evaluations, the estimated hourly, monthly and annual mean global irradiance corresponded with the observed values.

Global Solar Irradiance during 2000-2018
During 2000-2018, annual averages of observed global and diffuse solar irradiance (n= 92,941) were 0.59 and 0.29 MJ m −2 , thus, the direct horizontal irradiance was a little larger than the diffuse irradiance at the surface, average of S/G was 0.64, and averages of surface air temperature and relative humidity were 0.76 • C (ranged from −39.90 to 31.50 • C) and 80.31%, respectively. Applying the empirical model of global solar irradiance, hourly global irradiance was calculated for Sodankylä during 1 January 2000 to 31 December 2018. The input data were observed hourly global and diffuse solar irradiance for S/G and water vapor pressure at the ground. Global solar irradiance in December was too low and not considered in the later analysis, and n decreased to 69,171. The calculated and observed hourly global solar irradiance varied similarly, and calculated values were larger than the observed by 18.11% for the average in 2000-2018. It is reasonable that the empirical model represents the global irradiance and its relationships with absorbing and scattering processes at relative clean atmospheric conditions (i.e., relative low S/G at 0.50 in 2008-2011 in the model development). The calculated and observed monthly global irradiance, diffuse irradiance and S/G are shown in Figure 5. Global solar irradiance in June or July was higher over clean atmospheric conditions, e.g., 2001 to 2003, 2007 and 2011 at low S/G (<0.60), but lower over high GLP loading, e.g., 2008, 2009 and 2017 at high S/G (>0.60). Scattering factor S/G was lower during April to September and higher in October to March. During 19 years, the monthly averages of the calculated global irradiance decreased by 0.01% and diffuse irradiance increased by 0.11% per year, respectively. They were associated with the increases of S/G by 0.14% and water vapor by 0.06%. Air temperature and relative humidity increased by 0.83% (corresponding to 0.01 • C) and 0.23% for monthly average, respectively. On the average, the annual air temperature increased about 2.09 • C during 2000-2018, also indicating the Arctic warming [27].
Long term interannual variations were also investigated during 2000-2018 ( Figures 6 and 7). The calculated annual global solar irradiance decreased by 0.92% and diffuse irradiance increased by 1.28% per year, respectively. Both were associated with the increases of S/G by 1.73% and water vapor by 0.43% per year. Annual air temperature increased by 4.64% (corresponding to 0.07 • C), and relative humidity increased by 0.23% per year, respectively. The annual air temperature showed a higher correlation with water vapor (R = 0.79) than with relative humidity (R = 0.20). A good correlation was determined in annual water vapor pressure (E) and air temperature (T) and relative humidity (RH) as E = 0.293 × T + 0.024 × RH + 3.994 (R = 0.828, n = 19, α = 0.001).  Based on above as well as the following integrated results, the increase of air temperature is evidently associated with the increased GLPs, as well as water vapor. If the output of the Sun remains, solar global radiation at the ground would decrease and air temperature increase with the increase of atmospheric substances in the future.

The Losses of Global Solar Irradiance in the Atmosphere during 2000-2018
The global solar irradiance losses caused by absorbing and scattering GLPs (G LA and G LS ) were calculated by using A 1 (1-e −kWm ×cos(Z)) and A 2 (1-e −S/G ), respectively. The total loss G L is G LA + G LS . The monthly losses due to absorption G LA dominated the total loss G L and showed clear seasonal variations, lower in April-September and higher in October-March, whereas G LS did not have evident seasonal variation and most peaks appeared in October-February. During January 2000 to December, 2018, (1) monthly G LA decreased slightly by 0.003%, associated with the increase of water vapor by 0.06%; (2) monthly G LS increased by 0.06%, associated with the increase of S/G by 0.14%; and (3) monthly G L increased by 0.02% (Figure 8). The contributions of monthly mean absorbing and scattering losses (R LA and R LS ) to monthly mean total loss were 63.33% (ranged from 53.95 to 71.99%) and 36.67% (28.01-46.05%), respectively, during January 2000 to December 2018 (Figure 9), corresponding to the monthly averages of water vapor at 6.48 hPa (0.57-14.33) and S/G at 0.55 (0.30-0.94). Generally, R LA was lower in April-September and higher in October-March, whereas R LS varied inversely compared to R LA (i.e., most peaks appeared in October-February). The corresponding annual mean absorbing and scattering losses (R LA and R LS ) to annual mean total loss were 63.32% (60.09-65.48%) and 36.68% (34.52-39.04%), respectively.  To better understand the characteristics of global solar radiation, monthly averages of global solar irradiance and its controlling factors (i.e., E and S/G), meteorological variables, including air temperature, wind speed (v) and square of wind speed, were analyzed, together with their correlations during April-September.
During April-September in 2000-2018, the observed global solar irradiance also decreased a little (0.002%), and the observed diffuse irradiance increased by 0.11%. Air temperature increased by 0.05% (corresponding to 0.010 • C), and relative humidity increased 0.02% for monthly average. The calculated monthly mean global irradiance was larger than the observed by 17.64%, associating with low S/G (0.49) and water vapor pressure (8.83 hPa). The calculated monthly mean global irradiance decreased by 0.07%, corresponding to very little decrease of observed G, which was caused by the increases of S/G by 0.13% and water vapor by 0.03%. The monthly losses G LA , G LS and G L decreased by 0.0002%, 0.078% and 0.029%, respectively, associating with the increases of E and S/G.
The calculated annual global solar irradiance decreased by 0.80% and the observed diffuse irradiance increased by 1.37% per year respectively ( Figure 10). Both were associated with the increases of S/G by 1.55% and water vapor by 0.28% per year. Annual air temperature increased by 0.45% (corresponding to 0.04 • C), and relative humidity increased by 0.14% per year, respectively ( Figure 11). On the average, the annual air temperature increased 0.78 • C per year during April-September higher than 0.07 • C per year during a whole year in 2000-2018. The annual losses G LA decreased by 0.01%, associating with the increase of E; G LS increased by 0.88%, associating with the increase of S/G; total loss G L increased by 0.32% ( Figure 12). G L showed a similar variation with G LS , but complicated variation with G LA . There was a strong correlation between G L and G LS (R = 0.99) and a weak correlation between G L and G LA (R = −0.07). The contributions of annual mean absorbing and scattering losses (R LA , R LS ) to annual mean total loss were 63.32% (ranged from 60.96 to 65.48%), 36.68% (34.52-39.04%), respectively, corresponding to the annual averages of water vapor at 8.83 hPa (7.67-9.70) and S/G at 0.49 (0.42-0.60). "Annual mean" air temperature in April-September was 10.14 • C (8. 19-11.92). It is evident that G L was mainly contributed by absorption (G LA ) , but strongly modulated by scattering (G LS ), indicating the absorbing and scattering GLPs play different roles in the losses of global irradiance.   Comparing the annual contributions of the energy losses due to absorbing and scattering GLPs, the ratios of R LA /R LS were 1.72 and 1.62 for a whole year and for "April-September" at Sodankylä (their averaged S/G were 0.55 and 0.49), and 3.31 and 10.11 for S/G <0.80 and S/G ≥0.80 at Qianyanzhou (refer to as QYZ, a subtropical Pinus forest in China), respectively. It indicates that the high GLP loads results in larger absorbing loss than scattering loss at these two forest regions, i.e., the absorbing substances attenuate much global radiation than scattering substances, especially for QYZ with very high GLP loads.

Sensitivity Analysis
The response of global solar irradiance to changes of absorbing or scattering substances (i.e., E or S/G) was calculated by using Equation (1) (n = 3962), when another factor kept at its original values. The results are shown in Figure 13 and Table 2.  Table 2. Changing rates of global solar irradiance (%) caused by the changes of one factor (%), while another factor was kept at its original level. Hourly global solar irradiance at the surface increased/ decreased with the decrease/increase of the water vapor, implying that the increased absorbing GLPs result in much loss of global irradiance in the atmosphere and lower irradiance at the surface. Similarly, global solar irradiance increased/decreased with the decreases/increases of the atmospheric scattering GLPs (represented by S/G), indicating that the increased scattering GLPs lead to much loss of global irradiance in the atmosphere. Global solar irradiance was more sensitive to changes in the scattering factor (S/G) than the absorbing factor (E), reflecting that the changes of scattering GLPs (such as aerosols, SOA, clouds) have stronger effects on the global solar radiation than absorbing GLPs, e.g., the mean ratio between the response rate of G to S/G and that to E at different changing rates was about 7.1. The responses of global solar radiation to changes in both E and S/G were negative and nonlinear ( Figure 13 and Table 2).

Estimation and Comparison of Albedos at the TOA and the Surface
Albedos at the TOA and the surface are key factors influencing the energy balance, weather, and climate [28][29][30]. The Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) Edition 4.1 [31,32] provides the monthly shortwave flux and incoming solar flux on the TOA for all-sky, clear (for cloud free) sky and clear sky over 1 • ×1 • region (https://ceres.larc.nasa.gov/products.php?product=EBAF-Product accessed on 9 April 2021).
The annual albedos averaged from April to September at the TOA during 2008-2011 were calculated by using the up solar irradiance divided by down solar irradiance and the coefficients in Equation (1): The albedo at the surface was calculated by using the following equation: albedo at the surface = A 0 ×Atte/(A 1 e −kWm ×cos(Z) The attenuation rate (Atte) of the scattering term was used to estimate the attenuation of the reflection A 0 at the surface. In general, the albedos at the TOA and the surface obtained from the monthly satellite-retrieved radiance and using Equations (2)-(4) showed similar interannual variations under all skies during 2008-2011 (Supplementary Materials Figure S1). The satellite-retrieved albedos at the TOA (0.414) were larger than those calculated (0.268) with relative bias of 35.26% and absolute error of 0.146 for the 4-year average. The surface downward shortwave irradiance was underestimated by the satellite with mean relative bias of 42.61%, compared to these measured at Sodankylä. So, it was corrected by using the measured monthly shortwave irradiance, and the corrected albedo (0.290) was in closer agreement with the calculated by a smaller relative bias of 7.71% (Supplementary Materials Figure S2).
The calculated and satellite-retrieved annual surface albedos varied similarly with relative bias of 20.00% and absolute error of 0.05 for the 4-year average (Supplementary Materials Figure S1). The calculated and satellite-retrieved annual albedos at the TOA and the surface showed decreasing tendencies during 4 years. Both the calculated and satellite-retrieved albedos showed the similar characteristics, i.e., the albedos at the TOA were larger than the surface.
The above methods were applied to calculate the albedos at the TOA and surface in 2000-2018 (Supplementary Materials Figure S2). During April-September in 2000-2018, the calculated and satellite-retrieved annual albedos at the TOA displayed decreasing trends by 0.58% and 0.07%, respectively, and the corresponding annual albedos at the surface showed increasing trends by 1.24% and 0.10% per year, respectively. As a reference, the atmospheric substances (i.e., S/G) increased by 1.15% per year, water vapor increased by 0.10% per year, and the observed global solar irradiance increased by 0.27% per year during April-September from 2000 to 2018.
To get the accurate albedo, the satellite-retrieved albedo is needed to be corrected by using surface measured global solar irradiance. Considering that the satellite-retrieved upward and downward surface global irradiance was underestimated by 77.58% and 73.99%, respectively, for 2008-2011 and 2000-2018, no corrections were made for satelliteretrieved albedo at the surface. However, this introduced small errors for the albedo comparison. The time and space match is also the main reason for the errors in the comparison. It is reported that the error of the albedo retrieved by MODIS data in the shortwave band is 85.9% [33].
The decrease of albedo at the TOA in 2000-2018 may be caused by (1) the increases of absorbing and scattering GLPs in the atmosphere (represented by E and S/G), including NO 2 (Section 3.8) and absorbing aerosols; (2) the increased emissions of biogenic volatile organic compounds (BVOCs) caused by the increases of air temperature and PAR [34,35], as T and PAR increase with the increase of G simultaneously; (3) the enhanced direct absorption and indirect use of UV and visible radiation by BVOCs and their oxidation products in CPRs reacting with OH radicals [36].
The change of global solar radiation at the ground in 2000-2018 was associated with energy in (1) absorption and scattering caused by all kinds of GLPs in the atmosphere, including gases (CO 2 , CH 4 , N 2 O, NO 2 , SO 2 , O 3 , HCHO, etc.), water vapor, aerosols (black carbon, SOA, etc.), clouds and (2) the reflection by the Earth-atmosphere system at the TOA. All these factors contribute the solar dimming and/or brightening and are suggested to be considered together, not only these commonly recognized gases (e.g., CO 2 , CH 4 , N 2 O, and water vapor), aerosols (e.g., BC), air pollutants and clouds [37][38][39][40][41][42].
The albedos at the TOA and the surface can be obtained from the empirical model in this study (Equations (5)-(7)), as well as obtained from popular radiative transfer models that require more atmospheric inputs, including aerosol, cloud and water properties [43], and not from the current empirical models introduced in the introduction.

Further Evaluation of the Empirical Model
To evaluate the performance of the empirical model of global solar irradiance under all sky conditions, firstly, the RMSE values were calculated by using the hourly values and were 0.315, 37.70% for 2000-2007 (n = 29,356); 0.353, 45.38% for 2008-2011 (n = 14,525); and 0.329, 42.12% for 2012-2018. These RMSE values were a little larger than that in model development (Section 2.2). It is reasonable, because the measurement errors were larger under these situations. This and other empirical models also show lower performance for cloudy and overcast skies than clear skies [12][13][14][15][16][17][18][19][20]44], indicating that the high atmospheric GLPs (aerosols, clouds, etc.) cause large simulation uncertainties and the model simulation for cloudy sky conditions should be improved in the future.
It also reveals that the RMSE values as well as observation errors decreased with the increase of the observation samples (n). In more detail, the RMSE values were in the range  (1), the estimated total global irradiance at the TOA, was compared with the solar constant I 0 . To further investigate the performance of the empirical model, S/G was divided into different ranges and other parameters (G, S, E, etc.) were also divided with the same ranges as the S/G values. The ratio between the calculated G at the TOA and I 0 , relative error between calculated and measured G at the ground, RMSE were calculated and are shown in Table 3. Table 3. The ratio between the calculated G at the TOA and I 0 , relative error (δavg, δmax) and RMSE (MJ m −2 and %) between calculated and measured G at the ground. It provides that the ratio of calculated G at the TOA to I 0 decreased from 1.216 to 0.999 when S/G decreased from S/G ≤ 0.80 to S/G ≤ 0.30, and increased with the S/G further decrease after S/G ≤0.30, i.e., the empirical model can obtain the best simulation of global solar irradiance (ratio = 0.999) under the clean atmospheric condition (S/G ≤0. 30), and the sufficient sample number was also an important requirement (n = 1322). As the ratios were increased when S/G continue to decrease from S/G ≤0.20 and the sample numbers were smaller than 1322. Generally, the much clean of the atmosphere, the better simulation of G was manifested. There were relatively large estimation errors (δavg, δmax and RMSE) for the empirical model under high GLP loads (S/G ≥0.50), because the larger errors in measurement and model simulation influenced by aerosols, clouds, rain, fog, etc.

S/G Range
In summary, the simulations of the empirical model were in good agreement with the observations, according to the above performance at the ground and at the TOA.

Estimation of Aerosol Optical Depth (AOD)
Aerosol optical depth (AOD) from March 2013 to October 2017 was obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) [45] onboard NASA (National Aeronautics and Space Administration)'s Terra satellite. The daily AOD Dark Target (DT) and Deep Blue (DB) merged L2 product was used to calculate monthly AOD values.
The DB algorithm was developed for use over bright surfaces and all land types [46][47][48]. The MODIS/Terra AOD C6.1 data were averaged in a radius of 50 km around the study site [12,49].
The scattering of global solar radiation in the atmosphere is mainly contributed by aerosols, clouds, and water. In a sample way, hourly AOD, hourly mean S/G and water vapor pressure (E) from May 2013 to August 2016 were selected for the following analysis. A strong correlation between S/G and E and AOD was found (R = 0.874, n = 284, α = 0.001), suggesting an empirical equation: The statistical metrics were δavg = 6.93% and δmax = 30.76%. Then, this empirical equation was used to calculate AOD, and the δavg = 25.55% and δmax = 56.50%. The standard deviations were 0.07 and 0.08 for satellite-retrieved and calculated AODs, respectively.
To validate the empirical equation for estimating AOD, observational data of solar radiation and water vapor from April 2013 to September 2017 when the satellite AOD values are available were used. The estimated and satellite-retrieved monthly AOD varied in a similar pattern (Supplementary Materials Figure S3), their relative bias ranged from 9.95 to 262.31%, and the maximum of AOD difference was 0.11. Annual mean calculated and satellite-retrieved AODs were also in reasonable agreement with relative bias ranged from 2.92 to 33.39%, and their difference ranged from 0.002 to 0.009 (Supplementary Materials Figure S4). The calculated and satellite-retrieved annual AOD showed similar decreasing trend by 3.48% and 4.93% per year, respectively.
Then, the empirical equation of S/G and AOD was applied to estimate AOD from January 2000 to December 2018, using the observational data at Sodankylä. Monthly AODs from April to September were counted for annual average. Considering the satellite AODs were retrieved usually at low S/G, S/G < 0. 45  The correlation between S/G and E and AOD (R = 0.874) was larger than that between E and AOD (R = 0.523, n = 284), revealing the scattering substance can be mainly separated and expressed as water vapor and aerosols in the atmosphere. The estimated monthly and annual mean AODs varied similarly to the satellite-retrieved AODs. This sample method (Equation (5)) can be used to calculate AOD.
The uncertainties of AOD estimates are mainly from the different time (continuous estimates using Equation (5) and not continuous from the satellite) and space coverage caused by ground and satellite observations, clouds amount in over pass time for the satellite and during a day (i.e., for estimation of AOD without satellite measurement).
To investigate the reason of the increases of the calculated and observed AODs after 2013, we analyzed the monthly averaged total frequent occurrence of seven and three types of aerosols below 12 km derived from the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) [50], as well as air temperature in daytime during 2006-2018 (Supplementary Materials Figures S6-S9). The seven aerosol classifications are clean marine, dust, polluted continental/smoke, clean continental, polluted dust, elevated smoke, and dusty marine; and three types of aerosols, namely polluted continental/smoke, clean continental and elevated smoke, were selected to reflect continental and smoke aerosols to some extent.
The averaged AOD in April-September 2014 was higher than in 2013 and in 2015 well corresponding to the frequent occurrences of 3 and 7 type aerosols in August, but, not well with that during April-September (Supplementary Materials Figure S6). Monthly calculated and observed AODs peaked in August, 2014 and other years were in line with the monthly frequent occurrences of 3 and 7 type aerosols (Supplementary Materials Figure S7). The ratios of frequent occurrence of 3 to 7 type aerosols were larger in June, July, August with the average of 0.72, 0.69 and 0.78, and smaller in April, May and September with 0.57, 0.57 and 0.58 during 2007-2018, implying that the important continental and regional aerosol sources, e.g., transportation, smoke and oxidation products of BVOCs. The higher AOD in June-August was well associated with higher air temperature in July-August, together with the larger ratios of 3 to 7 type aerosols, indicating BVOC oxidation and their products in summer play dominant roles in regional aerosol formation. The growth rate of nucleation-mode particles has a clear maximum in summer [51] can be used as evidence.
The averaged wind speeds were 2.42 and 2.47 m s −1 , and the south wind was dominated at Sodankylä in each month during April-September in 2002-2012 and 2013-2018. BVOC emissions are mainly driven by air temperature [34,52], their oxidation is an important source of small size aerosols (see References [51,53,54] and references therein), considering the high air temperature usually appears in July, August, especially in 2014 (Supplementary Materials Figures S8 and S9), and the larger contributions of 3 to 7 type aerosols in June, July and August. Taiga forest is the largest terrestrial biological community in the world, covers between 50-degrees latitude north and the Arctic Circle, accounting for 29% of the world's forest coverage. In addition, biomass burning can lead to increase of BVOC emissions [35]. So, Taiga forest together with biomass burning (e.g., boreal forest fires and agricultural fires) may play important contributions in background aerosol at Sodankylä. The contribution of smoke from the boreal forest and agricultural fires is well documented ([50] and references therein]). The averaged AOD in April-September peaked in 2011, 2014 and 2016, and it also corresponded well to the higher air temperature in July-August. Based on above integrated results, the calculated AOD may be used as a reference to study long-term variations of AOD and SOA production in 2000-2018.

Estimation of Tropospheric NO 2 Vertical Column Densities
Tropospheric NO 2 vertical column densities (VCDs) from April 2007 to November 2016 were derived from GOME2. The overpass data were calculated for the center of the satellite pixel which includes Sodankylä. The space difference was up to 50 km from the station.
To develop an empirical model for estimation of tropospheric NO 2 , observational data of tropospheric NO 2 VCD, global solar radiation and water vapor during April 2008 to November 2011 were selected. Three key processes associated with global solar radiation transfer in the atmosphere were considered: (1) Absorption of NO 2 in UV and visible regions: it is expressed by e −k'cm , k' is the NO 2 absorption coefficients in UV and visible regions, and 13.01×10 4 Pa cm −1 [12,55]; c is tropospheric NO 2 VCD (cm −2 ), and m is optical air mass. (2) The absorption term: absorption and use of global solar radiation by GLPs, not including NO 2 . It is expressed as e −kWm ×cos(Z). The meaning of this term is discussed in detail in Reference [12]. The hourly global solar irradiance considering NO 2 absorption in all sky conditions was calculated by using the following empirical model: G = C 1 e −k cm + C 2 e −kWm × cos(Z) + C 3 e −S/G obs + C 0 This equation is converted for the estimation of NO 2 VCD, as follows: e −k cm = D 1 G + D 2 e − kWm × cos(Z) + D 3 e −S/G obs + D 0 (7) Considering that tropospheric NO 2 VCD was retrieved under low cloud fraction and the solar radiation data with S/G <0.90, NO 2 VCDs within 1 standard deviation were used in the empirical model development (Equation (6)). A strong correlation was found between G and three key terms, the correlation coefficient was 0.956 (n = 216, α = 0.001). The coefficients D i (i = 1, 2, 3) and the constant D 0 were −0.006, 0.481, 0.006, and 0.639, determined by using a multiparameter fit hourly NO 2 VCD from 8:00 to 19:00 during May 2008-September 2011. The main results are shown in Table 4. The calculated and observed NO 2 VCDs are presented in Supplementary Materials Figure S9. It is clear that most standard deviations of calculated NO 2 VCDs are less than the satellite derived (Supplementary Materials Figure S9) and similar results is shown in Table 4, σ cal < σ obs . In addition, model simulation of annual NO 2 in 2008-2011 was also within 1 standard deviation of the observed in empirical model development (Supplementary Materials Figure S10). It reflects that the reasonable estimations are determined by using the simultaneous NO 2 VCDs in empirical model development. To validate the empirical model of NO 2 VCD, the tropospheric NO 2 VCD from 8:00 to 19:00 during April 2008 to November 2011 (n = 1445) was calculated by using Equation (7) and the observational data. The calculated NO 2 VCD was underestimated by 47.09% with NMSE =1.76. It is caused by the data-selection criteria in model development, i.e., S/G <0.90 and 1 standard deviation of NO 2 VCD were used. In order to improve the simulation of hourly NO 2 VCD, the constant C 0 was adjusted to 0.535, and the relative bias reduced to −0.99%, NMSE = 0.73 (n = 1445). The relative biases ranged from −0.05 to 10.33% for annual NO 2 VCD in 2008-2011.
Then, these new coefficients were used to estimate tropospheric NO 2 VCDs from 8:00 to 19:00 during April 2000 to November 2018. The annual mean NO 2 VCDs were also obtained from May to October (Supplementary Materials Figure S10). The calculated and satellite-retrieved annual NO 2 VCD increased by 3.83% and 3.23% per year during 2007 to 2016, respectively. Generally, the empirical model displayed a reasonable estimation of hourly NO 2 VCD in annual average and its long-term variation. Tropospheric NO 2 VCD showed a slow increasing trend by 0.03% per year in 2000-2018.
Generally, the estimation of NO 2 VCDs were in agreement with the satellite-retrieved ones for annual average and long-term variation trend during 2007-2017, and the calculated NO 2 VCD had smaller uncertainties than satellite-retrieved values. It should be mentioned that the estimated NO 2 VCDs represent these for the conditions, S/G <0.90, 8:00-19:00 and within 1 standard deviation. The space and time match (continuously obtained from using the Equation (7) and not continuously from the satellite) at ground and satellite observations would introduce estimation errors of NO 2 VCDs. High clouds were not considered in satellite-retrieved NO 2 , but it was not for the estimation from using empirical model (i.e., S/G <0.90), especially all NO 2 VCDs analyzed from 8:00 to 19:00. Larger NO 2 VCDs appeared at S/G >0.90 were not used in empirical model development, causing the low annual NO 2 VCD estimations. The satellite observed mean annual NO 2 VCD increased 0.20×10 15

Empirical Model and Estimation of Global Solar Irradiance
An empirical model describing the energy relationships between global solar irradiance and its key processes was developed, based on analyzing representative radiation and meteorological variables measured in natural atmospheric conditions. Global solar radiation is attenuated by absorption and scattering, which are expressed by the absorbing and scattering terms, respectively. It can simulate the global solar irradiance both at the ground and at the TOA, the absorbing and scattering irradiance associated with the GLPs, as well as albedos at the TOA and the surface. This study is to comprehensively investigate the long-term variations of global solar irradiance, the losses of absorbing and scattering energy, albedos at the surface and the TOA, atmospheric constituents (e.g., NO 2 and AOD), and their integrated roles in the Earth-atmosphere system. The empirical model of global solar irradiance is an application of the previous empirical models in the UV, visible and short wavelength bands under all sky conditions [12].
Global solar radiation is attenuated by the atmospheric GLPs, influences the regional vertical and horizontal air movements; for example, the higher the GLP loads, the lower the wind speed. During 2000-2018, wind speed at Sodankylä increased by 0.01% for monthly averages (January-December) and 0.14% for annual averages (Supplementary Materials Figures S11 and S12) Supplementary Materials Figure S11 shows the monthly mean wind speed and the monthly mean S/G at Sodankylä during January 2001 to December 2018, and the evident behavior of decreasing wind speed and increasing S/G. Supplementary Materials Figure S13a,b shows monthly wind speed and monthly atmospheric substance for S/G interval at 0.05. These observations reveal that low wind speed associated with the high atmospheric substance loads (S/G). In other words, the increase of atmospheric GLPs, due to direct emission and production from CPRs [59], and transportation from other regions caused much loss of global radiation in the atmosphere, less energy reaching at the surface and long-wavelength radiation heating the atmosphere, and then the decrease of wind speed. Similar phenomenon and mechanism are also displayed at a subtropical region in the northern hemisphere [12].
Global solar radiation arriving at the surface and their absorbing and scattering losses showed obvious monthly and interannual variations (Section 3.3) at this Arctic region, and drives the atmospheric movement on regional scale. It is important to thoroughly investigate radiation transfer, absorbing and scattering processes and their distribution, reflections, and the interactions in radiation-GLPs in the whole atmospheric column.

Contributions of Absorbing and Scattering
The contribution of absorption (G LA ) to the total loss of global solar irradiance (G L ) was more important than that of scattering (G LS ) (i.e., annual averages of R LA and R LS were 63.33% and 36.67% in 2000-2018), meaning that the absorbing substances attenuate more global solar radiation than the scattering substances. In other words, the scattering GLPs play much less important roles in the short wavelength region, compared to the corresponding annual contributions (R LA and R LS ) of 35.3% and 67.7% in the UV region and 4.7% and 95.3% in the visible region in North China during 2004-2006. The UV radiation excites more GLPs to participate in CPRs, and more energy is consumed in the UV region (35.30%) than in the visible region (4.70%) [36,60]. The similar feature that R LA is much larger than R LS is also found in Qianyanzhou subtropical Pinus forest [12], revealing the same mechanism is existed at these two sites. The independence of absorbing and scattering terms in the UV, visible and short wavelength regions helps to describe individual absorbing and scattering process and energy. The absorbing term expresses total absorbing through different mechanisms in the UV, visible and near-infrared regions, and the scattering term expresses the total scattering due to all scattering GLPs [12]. In short, solar radiation in different wavelength regions exhibited different absorbing and scattering characteristics and distributions.

Long-Term Variations and Interactions in Global Solar
Irradiance and Its Losses, Albedos at the TOA and the Surface, NO 2 , AOD and S/G Solar radiation as an important energy input to the Earth-atmosphere system governs and influences air movements, GLP changes and their conversions between gas, liquid and particle phases. All these interactions and dynamic changes between solar energy and GLPs determine the regional climate and climate change.
Compared to the annual averages in 2000-2012, AOD increased double in 2013-2018 (Section 3.7). It was associated with the simultaneous increases of S/G, diffuse radiation and air temperature, and the decrease of global solar radiation. The main reasons are deduced as follows: (1) the increases of local BVOC emissions with the increase of air temperature [34,35], and the further increase of SOA formation from BVOC oxidations [53,54,61,62]; (2) a little increase of tropical NO 2 VCD (1.28%) favorites the chemical and photochemical production of aerosols; (3) long-distance transport of GLPs from the midlatitude regions, e.g., Southeast Asia [63][64][65]; and (4) much higher S/G at QYZ than at Sodankylä in 2013-2016. Similarly, the annual mean albedo at the surface and NO 2 VCD also increased by 7.58% and 1.28% during January-December in 2013-2018, compared to in 2000-2012, well corresponding to the increases of AOD and S/G. Therefore, all above processes resulted in a small decrease of G (0.69%) and increase of diffuse radiation (22.38%) in 2013-2018, compared to 2000-2012.
It should be discussed about the NO 2 empirical model, the absorption role of NO 2 in the UV and visible regions, its actual energy contribution in solar global radiation is quantified by using the observational data, no matter how small NO 2 is. For example, a similar UV energy balance method is used to simulate surface O 3 , considering the relationship in O 3 , NO, NO 2 , absorbing and scattering term, and UV radiation [66]. The actual roles of O 3 , NO and NO 2 are described in a UV balance equation, similar to this study, and surface O 3 variation and the photochemical interactions in O 3 and NO, NO 2 are well revealed; for example, O 3 responses to the changes of NO, NO 2 and water vapor correspond well with the known photochemical mechanisms, reflecting how the signals of NO 2 can be captured and displayed in the empirical model, though NO 2 absorption in the UV region and NO 2 concentration is smaller than O 3 . If NO 2 VCDs larger than 1 standard deviation (e.g., 2σ and 3σ) are used in the further development of the empirical model, larger NO 2 VCDs will be captured, especially for these after 2011. In other words, the development of the empirical model depends on the objectives of the investigation.
Based on the above results and other related studies in the mid-latitude and Arctic regions, the comprehensive investigations in the physical, chemical and photochemical, biological processes for the Earth-biosphere-atmosphere system on regional and global scales should be taken into consideration in the future.

Comparisons of Global Solar Radiation and Its Related Parameters at Two Sites in the Northern Hemisphere
To understand solar global radiation and the interaction between solar global radiation and atmospheric variables in the northern hemisphere, it is benefit to analyze and compare solar global radiation and its loss, other driving factors at two sites in the northern hemisphere. Qianyanzhou, a Pinus plantation at Taihe county, Jiangxi province, China (26 • 44 48 N, 115 • 04 13 E), was selected as a representative site in the mid-latitude region.
The annual averages of monthly global solar irradiance and other variables were calculated for Sodankylä and QYZ sites (in all sky conditions) during 2013-2016 (Table 5). The ratios of all parameters between Sodankylä and QYZ are also given in Table 5. The calculated global solar irradiance at the surface was 54.23% lower at Sodankylä than QYZ, the total loss of global solar irradiance was 63.08% larger at Sodankylä than QYZ, together with the albedo at the TOA was larger 24.05% at Sodankylä than QYZ, resulting in air temperature was −19.65 • C lower at Sodankylä than at QYZ. The annual averages of E and S/G were much lower at Sodankylä than that at QYZ, but, the longer optical length at Sodankylä than QYZ is a key reason to cause the large losses (G LA , G LS , G L ) of global solar irradiance. Table 5. Annual averages of monthly global solar irradiance and other variables at Sodankylä (refer to as Sod) and QYZ sites during 2013-2016, and the ratios of all parameters between Sodankylä and QYZ (alb and sur denote albedo and surface, respectively). Global solar irradiance at the surface decreased at QYZ in 2013-2016, which is mainly caused by the increases of S/G (emissions of GLPs, GLP products from BVOC oxidation and CPRs, etc.) and reflection at the TOA [12]. During 2013-2016, the annual mean loss of G LA , G LA and G L at QYZ were (1) 1.64, 0.49 and 2.13 MJ m −2 (corresponding to 455.22, 135.50 and 590.72 W m −2 ) for S/G <0.80 conditions, respectively; and (2) 1.69, 0.16 and 1.85 MJ m −2 for S/G ≥0.80, respectively. These results were obtained by similar calculation as this study, the other information related to BVOCs is reported in Reference [35]. Compared annual averages of Sodankylä during 2000-2018 (Section 3.2) and QYZ during 2013-2016, the ratios of G LA , G LS and G L between Sodankylä to QYZ were 1.18, 2.38 and 1.46, respectively, for S/G <0.80 at QYZ, and 1.15, 7.31 and 1.68, respectively, for S/G ≥0.80 at QYZ. Thus, much global solar radiation, about 46.01% and 68.11%, was attenuated by GLPs higher at Sodankylä than QYZ for these two atmospheric conditions. In addition, more scattered energy was also lost at Sodankylä than QYZ.
The annual contributions of absorbing and scattering losses to annual total loss were 76.77% and 23.23% under the situations of E at 22.11 hPa and S/G at 0.73, and 91.08% and 8.92% under the situations of E at 21.60 hPa and S/G at 0.88 for QYZ. Therefore, affecting factors (i.e., E and S/G) at different levels resulted in the different absorbing and scattering processes and energy lost at two sites. Much global solar energy was scattered in the atmosphere at Sodankylä than QYZ with a factor of 4.55, revealing that the clean atmosphere (lower E and S/G) strongly enhanced the scattering at Sodankylä than QYZ. There were large differences in annual absorbing and scattering losses and their ratios (G LA /G LS ) were 6.19 at QYZ, and 1.58 at Sodankylä, revealing much more absorbing GLPs than scattering GLPs at these two sites, and this phenomenon/effect was more evident at QYZ than Sodankylä, implying the impacts of their production from regional biogenic and anthropogenic activities. Thus, the calculated albedo at the TOA was larger at Sodankylä than at QYZ, and was similar for the albedo at the surface. The similar annual albedos at the surface were close during 4 years, indicating the reasonable estimations at these two forest regions in the north hemisphere. The larger annual albedo at Sodankylä than at QYZ reveals that the atmospheric substances were more scattering and less absorbing at Sodankylä, conversely, the atmospheric substances were more absorbing at QYZ (Table 4). These differences in the albedo at the TOA are mainly caused by the different regional atmospheric constituents, i.e., cleaner atmosphere at Sodankylä with low water vapor and other GLPs (E and S/G).
Due to much energy was lost in the atmosphere and larger albedo at the TOA at Sodankylä than at QYZ, the annual mean air temperature, 2. In general, the annual mean ratio of solar absorbing loss divided by S/G and then by air temperature (T) (Q LA /(S/G)/T) varied slightly at different S/G intervals (S/G < 0.65) during 19 years (Supplementary Materials Figure S14). To decrease the influence of random errors, air temperatures less than one standard deviation and other corresponding variables were removed in the analysis (n = 98). As the sample numbers of monthly ratios of G LA /(S/G)/T were small (i.e., 1 and 2 for S/G intervals 0.65-0.70 and 0.71-0.75, respectively), and these data were not taken into the calculation. Therefore, averaged ratio (G LA /(S/G)/T) was 0.32 MJ • C −1 for most representative atmospheric conditions. Finally, the relationship between G LA /(S/G) and air temperature (T) was determined at different S/G intervals, and it ranged from 0.30 to 0.75 (R = 0.68, α= 0.05, n = 98): It reveals that the energy of G LA for the united GLPs is partially converted to thermal energy, which increases the air temperature, and the other part is consumed in photochemical reactions without relation with air temperature, and is a constant (0.045). A relative equilibrium state exited in the interactions between absorbing solar energy divided by atmospheric GLPs and air temperature, using 19-year data, and the increase of absorbing solar energy of 0.273 MJ m −2 heated the atmospheric GLPs and increased the air temperature by about 1 • C. The calculated absorbing solar energy, thermal and photochemical energy increased with the increases of S/G (Supplementary Materials Figure S15), and their mean averages were 1.79, 1.76 and 0.02 MJ m −2 (corresponding to 495.70, 489.11 and 6.58 Wm −2 ), respectively. The mean contributions of thermal and photochemical energy to absorbing energy were 98.60% and 1.40%, respectively. The small photochemical energy (ranged from 4.14 to 8.90 W m −2 ) is deduced as UV and visible absorption by all kinds of atmospheric substances.
The total solar absorbing energy (including direct absorption and indirect utilization in UV, visible and near-infrared wavelength regions, such as the known absorbing gases, as well as VOCs, organic compounds, black carbon, etc.) is lost and utilized in the atmosphere should be considered in the study of regional climate and climate change.
Similar to Sodankylä, the relationship between G LA /(S/G) and air temperature was determined at different S/G intervals at QYZ, and it ranged from 0.70 to 1.00 (n = 40, R = 0.914, α = 0.05): For unite atmospheric GLPs, the increase of absorbing energy of 0.087 MJ m −2 resulted in the increase of temperature by about 1 • C at QYZ, i.e., which was a factor of 1/3 to that of Sodankylä, and the photochemical energy was 0.040 at QYZ, close to that at Sodankylä. It is clear that much absorbing energy lost in the atmosphere by a factor of 3.14 at Sodankylä than at QYZ results in the large increase of air temperature, and it may be an important reason for the Arctic amplification [1][2][3]. Comparing the annual changing rates (%) in 2000-2018 at Sodankylä and in 2013-2016 at QYZ [12], annual mean S/G increase rate is 1.73% at Sodankylä and 6.08% at QYZ and their ratio is 0.28, whereas annual mean air temperature increases at 2.09 • C for Sodankylä and 0.07 • C for QYZ and their ratio is 28.63, implying that the atmosphere at QYZ has larger total heat capacity than Sodankylä, because of the high atmospheric substances (S/G). Or, the low total heat capacity of the atmosphere at Sodankylä is also benefit to a larger increase in air temperature. Therefore, it is vital to control anthropogenic emissions of GLPs and secondary GLP products through CPRs for controlling regional climate warming.
The calculated absorbing solar energy, thermal and photochemical energy at QYZ also increased with the increase of S/G, and their mean averages were 1.79, 1.76 and 0.03 MJ m −2 (corresponding to 496.68, 487.70 and 9.00 Wm −2 ), respectively. The mean contributions of thermal and photochemical energy to absorbing energy were 98.17 % and 1.83%, respectively. The photochemical energy was also very low and ranged from 7.19 to 10.82 W m −2 . Similarly, for calculated solar global irradiance (G cal ) at the ground at Sodankylä, the corresponding ratio of G cal /(S/G)/T decreased with the increase of S/G by 7.66% during 19 years, indicating that no equilibrium state exists in the interaction between solar global radiation at the ground and the GLPs (Supplementary Materials Figure S16).
Similarly, a sensitivity analysis (Supplementary Materials Figure S17) was also performed for QYZ as Sodankylä, using a similar model as Equation (1) and observational data at QYZ [12,35]. At both Sodankylä and QYZ sites, global solar irradiance was more sensitive to the changes in the scattering factor than absorbing factor, and this phenomenon was more prominent at Sodankylä than at QYZ, as the ratio of the response rate to scattering factor to the ratio of the response rate to absorbing factor was 7.10 at Sodankylä and 1.57 at QYZ, indicating that the atmosphere was much more cleaner and drier at Sodankylä than QYZ, as S/G and E were 0.50 and 8.55 at Sodankylä and 0.71 and 23.96 at QYZ, and (S/G)/E was 0.06 at Sodankylä and 0.03 at QYZ. The atmospheric absorbing and scattering GLPs and their changes influence the changes and distributions of global solar radiation at the surface and in the atmosphere (Section 3).
Comparing the responses of G to absorbing and scattering factors between Sodankylä and QYZ, G was more sensitive to changes in the absorbing factor at QYZ than at Sodankylä, which is caused by the higher absorbing GLPs (represented by E, 22.11 hPa) at QYZ than at Sodankylä (E = 8.56 hPa); but, G was more sensitive to changes in the scattering factor at Sodankylä than QYZ, which is associated with the low solar altitude (i.e., long optical length) and low GLPs at Sodankylä (S/G = 0.50) than high solar altitude and high GLPs at QYZ (S/G = 0.73).

Conclusions
Solar energy (including global solar radiation and its losses in the atmosphere, the reflection at the TOA and the surface), the atmospheric constituents (including H 2 O, gases and aerosols) and their long-term variations are essential to thoroughly understand the physical and chemical processes in the atmosphere, biological processes, and regional climate and climate change at Sodankylä and the Arctic. An empirical model of global solar irradiance (EMGSI) under all sky conditions was developed by using observational data in 2008-2011, at Sodankylä. The calculated hourly global solar irradiance was in agreement with the observed at the ground and the TOA. This model was applied to calculate the global solar irradiance at the ground and its attenuation in the atmosphere from January 2000 to December 2018. The results showed clear monthly and interannual variations. The sensitivity test indicated that global solar irradiance was more sensitive to changes in scattering than absorbing. The responses of global solar irradiance to changes in absorbing and scattering factors were nonlinear and negative. The albedos at the TOA and the surface were calculated by using this empirical model, and the simulated albedos were in agreement with the corresponding satellite-retrieved values.
A good relationship between S/G and AOD was determined and applied to estimate AOD in 2000-2018. An empirical model for estimation of tropospheric NO 2 VCD was developed and used to estimate long-term variations of tropospheric NO 2 VCD in 2000-2018. Negative correlation between the scattering factor and wind speed was found. The enhanced AOD in 2014, especially in August, was probably attributed to Taiga forest and forest fires/biomass burning, i.e., BVOC emissions and BVOC oxidation play important roles in aerosol formation in the Sodankylä region.
During 2000-2018, the estimated global solar irradiance decreased by 0.92%, and diffuse irradiance increased by 1.28% per year. They were associated with the increases of S/G by 1.73% and E by 0.43% per year. Annual air temperature increased by 0.07 • C. Annual mean losses of global solar irradiance in the atmosphere caused by absorbing and scattering GLPs and total loss were 1.94, 1.17 and 3.11 MJ m −2 (539.82, 323.86 and 863.68 Wm −2 ), respectively. Annual losses of global solar irradiance due to absorbing and scattering substances decreased by 0.02% and increased by 0.72%, respectively, and total loss increased by 0.24% per year. The contributions of annual mean absorbing and scattering losses to total loss were 63.32% and 36.68%, respectively, i.e., more energy was lost in the atmosphere by absorption. The calculated albedos were smaller at the TOA than at the surface. Simulated and satellite-retrieved annual albedos decreased by 0.58% and 0.07% per year at the TOA, and increased by 1.24% and 0.10% per year at the surface, respectively. In 2000-2018, annual mean AOD and tropospheric NO 2 VCD increased by 8.23% and 0.03% per year, respectively.
Global solar radiation and its loss, and other variables at Sodankylä and QYZ in the northern hemisphere were compared. Evident differences were found and were beneficial for our better understanding of the interactions between solar radiation and GLPs, atmospheric movement, regional climate and climate change.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/atmos12060749/s1, Figure S1: Annual mean albedos averaged from April to September in 2008-2011 calculated and satellite retrieved under all sky conditions at Sodankylä station. TOA cal, TOA sat and TOA sat cor denote the calculated, satellite retrieved and corrected satellite retrieved albedos at the TOA, surface cal and surface sat denote the calculated and satellite retrieved albedos at the surface, Figure S2: Annual mean albedos averaged from April to September in 2000-2018 calculated and satellite retrieved under all sky conditions at Sodankylä station, Figure S3: Monthly AOD calculated versus satellite retrieved (AOD cal and AOD sate) with error bars indicating standard deviation of satellite retrieved AOD at Sodankylä station, Figure S4: Annual mean AOD calculated versus satellite retrieved with error bars indicating standard deviation of calculated AOD at Sodankylä station, Figure S5: Annual mean AOD calculated versus satellite retrieved with error bars indicating standard deviation of calculated AOD at Sodankylä station, Figure S6: Total frequent occurrence of 7 and 3 types of aerosols during April-September (left) in August (right), Figure S7: Monthly sums of frequent occurrence of 7 and 3 types of aerosols during April-September, Figure S8: Averaged air temperatures during July-August at Sodankylä, Figure S9: NO 2 VCDs calculated verses satellite retrieved (NO 2 cal and NO 2 sate) with error bars indicating standard deviation of satellite retrieved NO 2 VCD at Sodankylä station, Figure S10: Annual NO 2 VCDs calculated versus satellite retrieved with error bars indicating standard deviation of satellite NO 2 VCDs during 2000 to 2018, Figure  S11: Monthly wind speed and monthly atmospheric substance over Sodankylä during January 2000-December 2018, Figure S12: Annual wind speed (v) and monthly atmospheric substance over Sodankylä during 2000-2018, Figure S13: Monthly wind speed (v) and monthly atmospheric substance (S/G) for x-axis using S/G interval at 0.05 (left) and their scatter plot (right) during January 2000-December 2018, Figure S14: Ratio of solar absorbing loss to S/G and then to air temperature (T) at different S/G intervals, Figure S15: The absorbing solar energy, thermal and photochemical energy (G LA , thermal and photochemical) at different S/G intervals, Figure S16: Ratio of calculated global solar irradiance to S/G and then to air temperature (T) at different S/G intervals, Figure S17