Spatio-Temporal Variation in AOD and Correlation Analysis with PAR and NPP in China from 2001 to 2017

: Atmospheric aerosols can elicit variations in how much solar radiation reaches the ground surface due to scattering and absorption, which may a ﬀ ect plant photosynthesis and carbon uptake in terrestrial ecosystems. In this study, the spatio-temporal variations in aerosol optical depth (AOD) are compared with that in photosynthetically active radiation (PAR) and net primary productivity (NPP) during 2001–2017 in China using multiple remote sensing data. The correlations between them are analyzed at di ﬀ erent scales. Overall, the AOD exhibited a northeast-to-southwest decreasing pattern in space. A national increasing trend of 0.004 year − 1 and a declining trend of − 0.007 year − 1 of AOD are observed during 2001–2008 and 2009–2017. The direct PAR (PAR dir ) and di ﬀ use PAR (PAR dif ) present consistent and opposite tendency with AOD during two periods, respectively. The total PAR (PAR total ) shows a similar variation pattern with PAR dir . In terms of annual variation, the peaks of AOD coincide with the peaks of PAR dif and the troughs of PAR dir , indicating that aerosols have a signiﬁcant positive impact on PAR dir and a negative impact on PAR dif . Furthermore, the PAR dir has a stronger negative association with AOD than the positive correlation between PAR dif and AOD at national and regional scales, indicating that PAR dir is more sensitive to aerosol changes. The NPP has higher values in the east than in the west and exhibits a signiﬁcant increasing trend of 0.035 gCm − 2 day − 1 after 2008. The NPP has a negative correlation ( − 0.4–0) with AOD and PAR dif and a positive correlation (0–0.4) with PAR dir in most areas of China. The area covered by forests has the highest NPP-PAR correlation, indicating that NPP in forests is more sensitive to the PAR than is the NPP in grasslands and croplands. This study is beneﬁcial for interpreting the aerosol-induced PAR impact on plant growth and for predicting plant production on haze days.


Introduction
Aerosols-solid or liquid particles suspended in the air with diameters ranging from 0.001 to 100 µm-are one of the main sources of air pollution and play an important role in the atmospheric environment [1]. The amount of aerosols can be measured by aerosol optical depth (AOD), which is the integrated extinction coefficient over a vertical column of a unit cross-section [2]. Besides the influence on climate and public health, aerosol particles are believed to influence solar radiation by scattering and absorption and especially affect photosynthetically active radiation (PAR) (400-700 nm), which possibly affects plant growth and production [3,4]. Monitoring the variation of AOD is significant not only for the detection of solar radiation change but also for better understanding of the aerosol-induced radiation impact on plant growth.

Aerosol Data
The monthly combined Deep Blue and Dark Target (DB/DT) AOD dataset from 2001 to 2017 are derived from the MODIS C6 level-3 products (MOD08_M3). The MOD08_M3 product with 1° resolution were downloaded from NASA Level-1 and the Atmosphere Archive and Distribution System Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/search/). The combined DB/DT dataset is available starting from the fall of 2014. This dataset is retrieved based on both DT and DB algorithms and has improvements in the Deep Blue algorithm for AOD retrieval over the bright surface [6,17]. The combined dataset preserves the quality of AOD retrievals and has a wider spatial coverage over land by improving the approaches of pixel selection, cloud mask and surface reflectance estimation [18]. The accuracy of DB/DT AOD is validated and has an RMSE (root-mean-square error) of 0.133-0.169 against Aerosol Robotic Network (AERONET) observations [18].

PAR Data
The monthly direct PAR (PARdir) and diffuse PAR (PARdif) data under a clear sky at a 1° × 1° grid are obtained from the SYN1deg edition 4.1 product of Clouds and the Earth's Radiant Energy System (CERES). The total PAR (PARtotal) are calculated by the sum of PARdir and PARdif. The CERES PAR products are computed from the fluxes at the top of the atmosphere (TOA) and surface for the visible spectrum (400-700 nm) based on a radiation transfer model that takes clouds, temperature, water

Aerosol Data
The monthly combined Deep Blue and Dark Target (DB/DT) AOD dataset from 2001 to 2017 are derived from the MODIS C6 level-3 products (MOD08_M3). The MOD08_M3 product with 1 • resolution were downloaded from NASA Level-1 and the Atmosphere Archive and Distribution System Distributed Active Archive Center (DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/search/). The combined DB/DT dataset is available starting from the fall of 2014. This dataset is retrieved based on both DT and DB algorithms and has improvements in the Deep Blue algorithm for AOD retrieval over the bright surface [6,17]. The combined dataset preserves the quality of AOD retrievals and has a wider spatial coverage over land by improving the approaches of pixel selection, cloud mask and surface reflectance estimation [18]. The accuracy of DB/DT AOD is validated and has an RMSE (root-mean-square error) of 0.133-0.169 against Aerosol Robotic Network (AERONET) observations [18].

PAR Data
The monthly direct PAR (PAR dir ) and diffuse PAR (PAR dif ) data under a clear sky at a 1 • × 1 • grid are obtained from the SYN1deg edition 4.1 product of Clouds and the Earth's Radiant Energy System (CERES). The total PAR (PAR total ) are calculated by the sum of PAR dir and PAR dif . The CERES PAR products are computed from the fluxes at the top of the atmosphere (TOA) and surface for the visible spectrum (400-700 nm) based on a radiation transfer model that takes clouds, temperature, water vapor from Aqua and Terra geostationary satellites and reanalysis databases as input data [19,20]. The CERES Ed4A products are available starting in 2017, and the downward fluxes are improved by updating the inputs (cloud properties, temperature and humidity profiles, surface albedo and aerosol) [21]. The monthly mean downward fluxes from the SYN1deg product were demonstrated to have good accuracy with a bias of 3 Wm −2 compared with the ground observations at 85 globally distributed sites [22]. The CERES PAR data are available on the website of https://ceres-tool.larc.nasa.gov/ord-tool/jsp/SYN1degEd41Selection.jsp. Because the AOD data are available only under cloudless conditions, the PAR dir and PAR dif of clear sky conditions are used in this study.

NPP Data
The NPP is the stored carbon after subtracting the maintenance and growth respiration from the total amount of organic matter fixed by plants [23,24]. It can be utilized as a measure of crop yield, grass biomass, and forest productivity [25]. The NPP data used in this study are obtained from the monthly and eight-day PSN (net photosynthesis) products (MOD17A2_M_PSN and MOD17A2H), which are level 4 MODIS products and are resampled to a spatial resolution of 1 • × 1 • . The NPP of MOD17 is calculated as carbon mass per unit area per day [25,26]. The algorithm of MOD17A2 NPP relies on the basic light-use efficiency (LUE) model. The MOD17A2 NPP data are computed from the PAR, fraction of absorbed photosynthetically active radiation (fAPAR) and leaf area index (LAI) obtained from another MODIS product (MOD15A2), meteorological data from the Data Assimilation Office (DAO) and biome-wise coefficients of biome properties lookup table (the BPLUT) [23].

Land Cover Data
The land cover data at a 500 m spatial resolution is achieved from the MODIS Land Cover Climate Modelling Grid Product (MCD12Q1, Collection 6). The dataset of the International Geosphere-Biosphere Programme (IGBP) classification is used to aggregate information into three broad vegetation classes in this study-forests, grasslands and croplands. The evergreen needleleaf forest, evergreen broadleaf forest, deciduous needleleaf forest, deciduous broadleaf forest and mixed forest are aggregated to Forests. Grasslands are grasslands dominated by herbaceous annuals. Croplands are composed of croplands and mosaics of croplands and cropland/natural vegetation [27]. Because the AOD, PAR and NPP data have a spatial resolution of 1 • × 1 • while the land cover data has a 500 m spatial resolution, the AOD, PAR and NPP data were oversampled from 1 • × 1 • to 500 m using spline interpolation to be consistent with the land cover data and to avoid the mixture of different vegetation classes.

Seasonal Kendall Trend Test
The Mann-Kendall method is a non-parametric statistical test method of trend detection for time series data [28,29]. This method can estimate the tendency regardless of whether data follow a normal distribution and is less influenced by outliers compared with linear regression [30,31]. The original hypothesis H 0 states that the time series data have no trend, and the alternative hypothesis H 1 is that a trend exists in the dataset. However, the test result may be influenced by the seasonality of the time series data [32]. As shown in the boxplots of Figure 2, the AOD, PAR dir , PAR dif , PAR total and NPP exhibit evident seasonal characteristics, affecting the significance of long-term trend detection. Therefore, the Seasonal Kendall method is used in this study for the trend tests, and this method is insensitive to the seasonality of the data. The Seasonal Kendall method tests the trend in each season separately and then combines the results into one overall test [31]. As an adaption of the Seasonal Kendall test, although the Regional Kendall test considers the multiple locations across the entire area, the spatial clustering and autocorrelation could be challenging issues. The test statistic for season i (i = 1, …,12) Si is defined as [32]: where is the number of samples in season i, and are the jth and kth sample in season i, sign( − ) is −1, 0 or 1 when is smaller, equal or larger than , respectively. The expectation of Si (E( )) is 0 and the variance ( ) is calculated by Equation (2): where is the sample size of a given tie. The distribution of Si is normal under the limit as → ∞, then Seasonal Kendall statistic S is defined as the sum of Si (∑ S ). The total variance Var(S) is derived as E(S) = 0 The data are assumed to be independent so that the covariance terms are equal to 0. The standardized test Z-statistic indicates the significance of the trend. When |z| > Z / (α is the Assuming the time series sample X is made up of seasonal subsamples Xi (i = 1, 2, . . . 12), the null hypothesis H 0 ' for the Seasonal Kendall is that the subsample Xi consists of independent and identically distributed random variables, while the alternative hypothesis H 1 ' is that an increasing or decreasing tendency exists in one or more season [32,33].
The test statistic for season i (i = 1, . . . ,12) S i is defined as [32]: where n i is the number of samples in season i, x ij and x ik are the jth and kth sample in season i, sign x ij − x ik is −1, 0 or 1 when x ij is smaller, equal or larger than x ik , respectively. The expectation of Si (E(S i )) is 0 and the variance Var(S i ) is calculated by Equation (2): Remote Sens. 2020, 12, 976 6 of 20 where t i is the sample size of a given tie. The distribution of Si is normal under the limit as n i → ∞, then Seasonal Kendall statistic S is defined as the sum of Si ( 12 i=1 S i ). The total variance Var(S) is derived as The data are assumed to be independent so that the covariance terms are equal to 0. The standardized test Z-statistic indicates the significance of the trend. When |z| > Z 1−α/2 (α is the significance level and is 0.05 in this study), the time series data are considered to have an increasing or decreasing trend at the significance level of 0.05. The term Z 1−α/2 is the percentile of the normal distribution of 1 − α/2.
The non-parametric robust estimate of the trend (β i ) in season i can be calculated by the median value.
The overall trend β is the median of the slope values (β i ). The trend is increasing when β > 0 and decreasing when β < 0. Figure 3 compares the 17-year annual mean distributions of AOD, PAR and NPP in China from 2001 to 2017. The annual mean AOD over China spatially varies in the range of 0-1.08. The AOD level in the southeastern coastal area is higher than that in the northwestern area, which is inversely correlated with the terrain variation. This is because lower terrain areas have more air pollutant emissions from anthropogenic activities, while areas with high elevation are sparsely populated and have high vegetation cover. The highest AOD level (>0.75) occurs in Beijing, Hebei, Henan, Shandong, Shanghai, Jiangsu, Sichuan Basin and Hubei provinces with a dense population and developed industry. The aerosols in these areas arise mainly from anthropogenic emissions. In addition, the mountains surrounding the Sichuan Basin inhibit the dispersion of aerosols [30]. The high AOD level in the west appears in the Taklimakan Desert of southern Xinjiang, and this high level may be caused by sand dust [5]. The regions with low aerosol loadings are distributed in the Tibetan Plateau, north of Xinjiang and Inner Mongolia, as these areas are underdeveloped and have a sparse population density.

Annual Distribution
The annually averaged PAR dir and PAR dif in Figure 3b,c show exactly opposite distributions. The PAR dir across China ranges from 34 Wm −2 to 108 Wm −2 , with higher values in the northwest and lower values in the southeast, which is generally opposite the observed trend of the AOD distribution. The PAR dif (17-60 Wm −2 ) shows a similar spatial pattern as the AOD, which increases from northwest to southeast. The Tibetan Plateau has the highest PAR dir and lowest PAR dif due to the high altitude (over 4000 m), high sunlight ranges (1600-3400 h/year), low water vapor contents and thin air [34,35]. The reduction of PAR in Tibetan Plateau in the process of radiation transmission is smaller than that in another area because of a shorter pathway through the atmosphere. Moreover, the PAR is less influenced by the absorption and scattering of aerosols due to the thinner air in the Tibetan Plateau [30]. The lowest PAR dir occurs in Heibei, Henan, Shandong, Jiangsu, Anhui, Hubei and Hunan provinces, Remote Sens. 2020, 12, 976 7 of 20 which are also the regions with the most aerosol pollution, indicating that aerosols may inhibit PAR dir . The highest PAR dif is dominant in southern China, including Guangdong, Guangxi, Hunan, Guizhou, and Jiangxi, where the AOD is larger than 0.4 and the water vapor content is higher than that in northern China. The PAR total decrease from southwest to northeast. The lowest level occurs in the southeast due to the lower PAR dir and PAR dif . The larger values of PAR total in the south of China can be attributed to large PAR dif . The largest PAR dir in Tibetan Plateau determines the largest PAR total there.
The averaged NPP in China ranges from 0-4 gCm −2 day −1 , with a mean value of 1.26 gCm −2 day −1 . The annual mean NPP (Figure 3d) shows higher values in the southeast than in the northwest. The highest NPP level above 2.5 gCm −2 day −1 , is mainly distributed in tropical humid regions, including Zhejiang, Fujian, Guangdong and western Yunnan. The white area represents no values, where there is essentially no vegetation cover. Eastern Inner Mongolia and Qinghai-Tibetan Plateau show the lowest NPP values under 0.6 gCm −2 day −1 , which may be due to the harsh climate and the grassland cover [12]. However, for individual provinces, the NPP does not show an obvious spatial pattern with AOD or PAR. For example, the North Plain and Sichuan Basin have the highest AOD, but the NPP values there are at the medium level. This discrepancy could be attributed to other factors, such as vegetation type, temperature, precipitation, and human irrigation. The regions with high NPP coincide with the forest distribution (Figure 3a), especially in summer. Yunnan has the highest NPP due to the higher cover of evergreen broad-leaved forests and tropical rainforests than that found in other regions in China.  The averaged NPP in China ranges from 0-4 gCm −2 day −1 , with a mean value of 1.26 gCm −2 day −1 . The annual mean NPP (Figure 3d) shows higher values in the southeast than in the northwest. The highest NPP level above 2.5 gCm −2 day −1 , is mainly distributed in tropical humid regions, including Zhejiang, Fujian, Guangdong and western Yunnan. The white area represents no values, where there is essentially no vegetation cover. Eastern Inner Mongolia and Qinghai-Tibetan Plateau show the lowest NPP values under 0.6 gCm −2 day −1 , which may be due to the harsh climate and the grassland cover [12]. However, for individual provinces, the NPP does not show an obvious spatial pattern with Remote Sens. 2020, 12, 976 8 of 20 AOD or PAR. For example, the North Plain and Sichuan Basin have the highest AOD, but the NPP values there are at the medium level. This discrepancy could be attributed to other factors, such as vegetation type, temperature, precipitation, and human irrigation. The regions with high NPP coincide with the forest distribution (Figure 3a), especially in summer. Yunnan has the highest NPP due to the higher cover of evergreen broad-leaved forests and tropical rainforests than that found in other regions in China.

Seasonal Variation
The seasonal variation characteristics of AOD, PAR and NPP in China are shown in Figure 4. Although the magnitudes vary in different seasons, the spatial patterns of AOD, PAR and NPP in each season are generally consistent with their respective annual mean distributions. As shown in Figure 4a, greater aerosol loadings are observed during summer and spring than in autumn and winter. The maximum AOD occurring in spring might be caused by large amounts of dust and biomass burning. High water vapor content possibly contributes to the second-highest AOD in summer [5]. Therefore, the seasonal variations in AOD closely change with emission sources, climatic factors and anthropogenic activities.
The PAR are found to be largest in summer due to the smallest solar zenith, and the values are second-highest in spring and much lower in autumn and winter. The PAR dir and PAR dif generally show opposite and similar spatial patterns with AOD in each season, respectively, which is consistent with the comparison of annual average distribution.
Similar to the AOD seasonal variation characteristics, the NPP is highest in summer, followed by that in spring and autumn, and lowest in winter. The NPP in the southeast is higher than that in the northwest in all seasons. The highest NPP in each season, especially in spring and summer, appears in the area covered by forests.  [5]. In general, the PAR dir and PAR dif show similar and opposite tendencies, respectively, with AOD in two periods, indicating the positive impact of aerosols on PAR dif and negative impact on PAR dir . Moreover, when comparing the AOD and PAR annual variation, the averaged AOD reaches its peaks in 2003, 2006, 2008, 2011 and 2014, which coincide with the peaks of PAR dif and the troughs of PAR dir and PAR total . Note that the trend magnitude of PAR dir (−0.515 Wm −2 before 2008 and 0.238 Wm −2 after 2008) is almost two times as much as that of PAR dif (0.249 Wm −2 before 2008 and −0.096 Wm −2 after 2008). The long-term trend and its spatial distribution of PAR total are consistent with PAR dir because the PAR dir has a dominant impact on PAR total due to higher value and larger variation magnitude than PAR dif . From 2001 to 2008, the NPP does not show a significant trend in most areas. The NPP shows an increasing trend of 0.035 gCm −2 day −1 at the 95% confidence level in most of China during the second period. The trend of NPP in the east (0.05-0.15 gCm −2 day −1 ) is larger than that in the west. Note that NPP has a greater increase from 2012 to 2013 than in other years. This partially because the increase of both PAR dif and PAR dir causes a large rise of PAR total . each season are generally consistent with their respective annual mean distributions. As shown in Figure 4a, greater aerosol loadings are observed during summer and spring than in autumn and winter. The maximum AOD occurring in spring might be caused by large amounts of dust and biomass burning. High water vapor content possibly contributes to the second-highest AOD in summer [5]. Therefore, the seasonal variations in AOD closely change with emission sources, climatic factors and anthropogenic activities. The PAR are found to be largest in summer due to the smallest solar zenith, and the values are second-highest in spring and much lower in autumn and winter. The PARdir and PARdif generally show opposite and similar spatial patterns with AOD in each season, respectively, which is consistent with the comparison of annual average distribution.

The Correlations between AOD, PAR and NPP
The pixel-wise correlations (R) between AOD, PAR and NPP are calculated using de-seasonalized data by removing the seasonal mean values. The correlation and significance maps are depicted in Figure 6. As shown in Figure 6a,b, the AOD has negative correlations with the PAR dir and PAR total and positive correlations with the PAR dif (significant level p < 0.05) in all of China except the Tibetan Plateau. The AOD over eastern China was found to have stronger negative correlations with PAR (−0.8 < R < −0.4 for PAR dir and PAR total , and 0.4 < R < 0.8 for PAR dif ) than that in the western area (−0.4 < R < −0.2 for PAR dir and PAR total , and 0.2 < R < 0.4 for PAR dif ). This result implies that aerosols have more impacts on the PAR in eastern China than on the PAR in the west due to higher aerosol concentrations. The correlation between AOD and PAR over the Tibetan Plateau is insignificant. This result is because the aerosol loading there is low during the entire year ( Figure 4) due to the high elevation and thin air; thus, the PAR changes in this area are more attributable to variation in the solar zenith angle rather than to that in aerosol [35]. The AOD and NPP show a significant negative correlation (−0.6 < R < −0.2) over the middle area of China. The PAR dir and PAR total PAR generally have positive correlations with NPP while the diffuse components are negatively correlated with NPP. Correlations with p < 0.05 between PAR and NPP (0 < R PARdir_NPP < 0.6, 0 < R PARtotal_NPP < 0.6 and −0.6 < R PARdir_NPP < 0) are observed in most of China except the northeast and Tibetan Plateau.

The Correlations between AOD, PAR and NPP
The pixel-wise correlations (R) between AOD, PAR and NPP are calculated using deseasonalized data by removing the seasonal mean values. The correlation and significance maps are depicted in Figure 6. As shown in Figure 6a,b, the AOD has negative correlations with the PARdir and PARtotal and positive correlations with the PARdif (significant level p < 0.05) in all of China except the Tibetan Plateau. The AOD over eastern China was found to have stronger negative correlations with PAR (−0.8 < R < −0.4 for PARdir and PARtotal, and 0.4 < R < 0.8 for PARdif) than that in the western area (−0.4 < R < −0.2 for PARdir and PARtotal, and 0.2 < R < 0.4 for PARdif). This result implies that aerosols have more impacts on the PAR in eastern China than on the PAR in the west due to higher aerosol concentrations. The correlation between AOD and PAR over the Tibetan Plateau is insignificant. This result is because the aerosol loading there is low during the entire year ( Figure 4) due to the high elevation and thin air; thus, the PAR changes in this area are more attributable to variation in the solar zenith angle rather than to that in aerosol [35]. The AOD and NPP show a significant negative Because the surface irradiance and the spectral distribution depend significantly on the solar zenith (θ), the effects of aerosols on PAR and NPP are analyzed for narrow θ ranges. Figure 7 presents the scatterplots between the AOD, PAR and NPP for different solar zenith angles. As seen in Figure 7a, the correlation (R) between AOD and PAR dir ranges from −0.77 and −0.52 when the solar zenith angle increases from 14 • to 61 • . The exponential model (PAR dir = a × exp(b × AOD)) is used to describe the overall relationship between AOD and PAR dir , which is consistent with other studies [36]. The PAR total shows negative correlations with AOD (−0.85 ≤ R ≤ −0.53), which is consistent with PAR dir . The correlation coefficient between PAR dif and AOD ranges from 0.42 to 0.65, and their relationship can be fit by a power model (PAR dif = a × AOD b + c). As shown in Figure 7c, the AOD has a low correlation with the NPP. The NPP can be positively or negatively correlated with the PAR dir , PAR dif , PAR total depending on the different solar zenith angles. The correlations of PAR dir -NPP and PAR dif -NPP are opposite for different solar zenith angles because when the PAR dir increases, the PAR dif decreases accordingly.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 21 correlation (−0.6 < R < −0.2) over the middle area of China. The PARdir and PARtotal PAR generally have positive correlations with NPP while the diffuse components are negatively correlated with NPP. Correlations with p < 0.05 between PAR and NPP (0 < RPARdir_NPP < 0.6, 0 < RPARtotal_NPP < 0.6 and −0.6 < RPARdir_NPP < 0) are observed in most of China except the northeast and Tibetan Plateau.

Vegetation Cover
The variation and correlation analysis were also conducted for AOD, PAR and NPP in different vegetation covers and regions. As shown in Figure 8a, the croplands are mainly located in eastern China, the grasslands dominate the western area and the forests exist in southern China. Remarkable differences in the AOD, PAR and NPP levels are found for different vegetation covers. The crop

Vegetation Cover
The variation and correlation analysis were also conducted for AOD, PAR and NPP in different vegetation covers and regions. As shown in Figure 8a, the croplands are mainly located in eastern China, the grasslands dominate the western area and the forests exist in southern China. Remarkable differences in the AOD, PAR and NPP levels are found for different vegetation covers. The crop production area is heavily polluted with the highest annual AOD, followed by that of forests, and lowest for grasslands (Figure 8b). The PAR dir and PAR total is lowest in croplands, moderate in forests, and highest in grasslands (Figure 8c and Table 1). The PAR dif in forests is highest, with a mean value of 46.867 Wm −2 , which is slightly higher than that in croplands (45.964 Wm −2 ) (Figure 8d and Table 1). Grasslands have much lower PAR dif values than those in other vegetation covers. The AOD in cropland has the largest increase before 2008 and the largest decrease after 2008. As shown in Table 1, the vegetation cover with a larger AOD increase has the larger descending trend of PAR dir and ascending trend of PAR dif during 2001-2008, and vice versa during the second period. Although the PAR total follows a similar rule with PAR dir , the trend magnitude is smaller than PAR dir . This is because the PAR dif offset the change of PAR dir . The correlation coefficients of AOD-PAR dir are larger than those of AOD-PAR dif , indicating that aerosols have a stronger negative impact on PAR dir than the positive impacts on PAR dif . This result suggests that the aerosol declines will lead to a greater increase in PAR dir than the decrease in PAR dif , which is consistent with the larger trend of PAR dir inTable 1, Table 2 and Figure 5. The correlations between the AOD and PAR are highest in croplands because the higher aerosol loadings in crop production areas can lead to greater change in PAR.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 21 production area is heavily polluted with the highest annual AOD, followed by that of forests, and lowest for grasslands (Figure 8b). The PARdir and PARtotal is lowest in croplands, moderate in forests, and highest in grasslands (Figure 8c and Table 1). The PARdif in forests is highest, with a mean value of 46.867 Wm −2 , which is slightly higher than that in croplands (45.964 Wm −2 ) (Figure 8d and Table  1). Grasslands have much lower PARdif values than those in other vegetation covers. The AOD in cropland has the largest increase before 2008 and the largest decrease after 2008. As shown in Table  1, the vegetation cover with a larger AOD increase has the larger descending trend of PARdir and ascending trend of PARdif during 2001-2008, and vice versa during the second period. Although the PARtotal follows a similar rule with PARdir, the trend magnitude is smaller than PARdir. This is because the PARdif offset the change of PARdir. The correlation coefficients of AOD-PARdir are larger than those of AOD-PARdif, indicating that aerosols have a stronger negative impact on PARdir than the positive impacts on PARdif. This result suggests that the aerosol declines will lead to a greater increase in PARdir than the decrease in PARdif, which is consistent with the larger trend of PARdir inTable 1, Table 2 and Table 5. The correlations between the AOD and PAR are highest in croplands because the higher aerosol loadings in crop production areas can lead to greater change in PAR.   Table 1 show that NPP has an insignificant trend in the first period but significantly increases after 2008 in three vegetation cover type areas. The forests have higher mean levels of NPP (1.808 gCm −2 day −1 ) and larger increasing magnitudes (0.074 gCm −2 day −1 ) than those of croplands and grasslands during 2009-2017. In contrast to PAR, NPP has no obvious relationship with AOD at the annual scale. The correlation between AOD and NPP (Table 3) is much lower than that between AOD and PAR for the three vegetation covers, indicating that AOD has a negative impact on NPP but that impact is not great. This may because, in addition to the PAR, the atmospheric CO2 concentration, atmospheric nitrogen deposition, ozone concentrations, land cover, and climate    Figure 8c and Table 1 show that NPP has an insignificant trend in the first period but significantly increases after 2008 in three vegetation cover type areas. The forests have higher mean levels of NPP (1.808 gCm −2 day −1 ) and larger increasing magnitudes (0.074 gCm −2 day −1 ) than those of croplands and grasslands during 2009-2017. In contrast to PAR, NPP has no obvious relationship with AOD at the annual scale. The correlation between AOD and NPP (Table 3) is much lower than that between AOD and PAR for the three vegetation covers, indicating that AOD has a negative impact on NPP but that impact is not great. This may because, in addition to the PAR, the atmospheric CO 2 concentration, atmospheric nitrogen deposition, ozone concentrations, land cover, and climate factors (temperature and precipitation) also drive the NPP change [12]. The PAR has stronger correlations (0.215, −0.218, 0.187 for PAR dir , PAR dif , and PAR total ) with NPP in forests than in croplands and grasslands ( Table 3), implying that NPP is more sensitive to the change in PAR in forests than in other vegetation cover types. This result is because greater PAR tends to penetrate more light into deeper canopies of forests, to increase the photosynthesis in shaded leaves; thus, the NPP of forest increases more than that in croplands and grasslands [37,38].

Typical Regions
As shown in Table 2, the mean AOD in five typical regions from 2001 to 2017 is almost 2 times (0.599-0.684) higher than that in the entire country (0.337), which is attributed to the increased anthropogenic activities caused by the rapid urbanization and expansion of industries. The AOD in all typical regions has a larger increase and decrease than that at the nationwide level in the first and second periods, respectively. The decrease of AOD during 2009-2017 in five typical regions is partially caused by emission reduction including NO 2 and SO 2 due to the air quality policy in China [39]. The decreasing trend of carbonaceous aerosols (e.g., organic carbon, elemental carbon) after 2008 and dust aerosol after 2010 may also contribute to the decline of AOD [40,41]. In addition, AOD variation may be influenced by thin cloud contamination in the satellite retrievals of AOD, meteorological factors, climate change and economic recession [42,43]. The greatest increase of AOD occurs in NP, accompanied with a largest decrease in PAR dir and PAR total and a second largest increase in PAR dif from 2001 to 2008. After 2008, the largest decreasing AOD trend in SCB corresponds with the largest PAR dir increase and PAR dif decrease. The correlation between the AOD and PAR is highest in CC and lowest in SCB ( Table 4). The AOD shows stronger negative correlations (−0.569-−0.342) with PAR dir than the positive correlations (0.237-0.429) with PAR dif in all typical regions, indicating that PAR dir is more sensitive to changes in AOD. This result is consistent with the larger change of PAR dir than that of PAR dif . In the second period, opposite to the downward tendency of AOD, the NPP presents a significant upward trend in all typical regions. GD is the most productive region with the highest mean NPP of 1.833 gCm −2 day −1 and the largest trend of 0.114 gCm −2 day −1 after 2008. The negative correlation of AOD-NPP ranges from −0.234 to −0.177 and is higher in CC, SCB and GD, which correspond to the regions with a high correlation between PAR and AOD.

Conclusions
In this paper, the spatio-temporal distribution, variation and trends of AOD was analyzed and compared with that of PAR and NPP from 2001 to 2017. We also investigated their correlations at the pixel, regional and national scales. Generally, the AOD varies from high in the southeast to low in the northwest depending on the different elevations and populations. The highest AOD levels are found in the North Plain, Yangtze River Delta, Central China and Sichuan Basin, which are areas with dense anthropogenic activities. Generally, the PAR dif spatial distribution shows a similar pattern as that of AOD, while the PAR dir and NPP have opposite distributions as that of AOD at the national scale. The PAR dir in the Tibetan Plateau is highest, and the PAR dif is lowest mainly because of the higher elevation. The NPP is higher in the southeast than in the northwest, with the highest values in the south area and the lowest values in the Tibetan Plateau.
The inter-annual trend investigation reveals that the AOD has an increasing trend from 2001 to 2008 and a decreasing trend during 2009-2017. The PAR dir and PAR dif have similar and opposite trends in two periods with AOD, respectively. The PAR total has a consistent tendency with PAR dir . The annual peaks of PAR dif and troughs of PAR dir and PAR total match exactly the peaks of AOD during the study period. The sulfate and absorbing organic aerosols (e.g., soot) originating from human activities have a great impact on the fluctuation of the PAR [44,45]. The NPP decreases first and then increases significantly with an amplitude of 0.035 gCm −2 day −1 after 2008, which is opposed to the temporal variation of AOD. The area covered by forest has the highest mean NPP and the largest increase rate during 2009-2017.
The correlations between AOD, PAR and NPP vary with solar zenith, location and vegetation cover. It is shown that the PAR dir and PAR dif have negative (−0.482) and positive correlations (0.408) with AOD across China except in the Tibetan Plateau. Overall, the correlation between AOD and PAR is higher in the east than in the west. The results of correlation analysis, together with the spatio-temporal variation comparison between AOD and PAR, demonstrate that the aerosol particles weaken the PAR dir and enhance the PAR dif . Furthermore, the PAR dir has stronger correlations with AOD than does PAR dif in most eastern areas of China and in different typical regions and vegetation covers, indicating that the aerosols have a stronger negative impact on PAR dir than the positive impacts on PAR dif for typical regions and main vegetation covers. The NPP is positively correlated with the PAR dir and PAR total and negatively correlated with the PAR dif and AOD. The correlations of NPP-PAR and NPP-AOD are less considerate than those of PAR-AOD. The correlation between PAR and NPP is highest in forests, implying that the NPP of forests is more sensitive to PAR than is the NPP in grasslands and croplands. The forest has the greatest NPP increase after 2008 because the PAR dif in forest can increase the level of penetration into the canopy, possibly resulting in increased radiation use efficiency.
This study provides a cohesive understanding of the spatio-temporal variation of aerosol loading and its relationships with solar radiation and plant production, which also provides a perspective of the interaction between aerosols and the ecological environment. The mechanism of aerosols affecting NPP is complicated and influenced by other climate factors, which is required for further study.