Effects of Aerosols on Gross Primary Production from Ecosystems to the Globe

: Aerosols affect the gross primary productivity (GPP) of plants by absorbing and scattering solar radiation. However, it is still an open question whether and to what extent the effects of aerosol on the diffuse fraction (D f ) can enhance GPP globally. We quantiﬁed the aerosol diffuse fertilization effect (DFE) and incorporated it into a light use efﬁciency (LUE) model, EC-LUE. The new model is driven by aerosol optical depth (AOD) data and is referred to as AOD-LUE. The eddy correlation variance (EC) of the FLUXNET2015 dataset was used to calibrate and validate the model. The results showed that the newly developed AOD-LUE model improved the performance in simulating GPP across all ecosystem types (R 2 from 0.6 to 0.68), with the highest performance for mixed forest (average R 2 from 0.71 to 0.77) and evergreen broadleaf forest (average R 2 from 0.34 to 0.45). The maximum LUE of diffuse photosynthetic active radiation (PAR) (3.61 g C m − 2 MJ − 1 ) was larger than that of direct PAR (1.68 g C m − 2 MJ − 1 ) through parameter optimization, indicating that the aerosol DFE seriously affects the estimation of GPP, and the separation of diffuse PAR and direct PAR in the GPP model is necessary. In addition, we used AOD-LUE to quantify the impact of aerosol on GPP. Speciﬁcally, aerosols impaired GPP in closed shrub (CSH) by 6.45% but enhanced the GPP of grassland (GRA) and deciduous broadleaf forest (DBF) by 3.19% and 2.63%, respectively. Our study stresses the importance of understanding aerosol-radiation interactions and incorporating aerosol effects into regional and global GPP models. The results show that the AOD-LUE model has higher potential in some poorly simulated vegetation types.


Introduction
Accurate and reliable estimates of gross primary productivity (GPP) are important for understanding the terrestrial carbon cycle and predicting plant production status [1,2]. Gross primary productivity is the major indicator to measure the material production (e.g., food and fiber) capacity and carbon uptake rate of an ecosystem [3,4]. Aerosols can affect the productivity of plants by increasing the diffuse radiation reaching the surface, thereby affecting the carbon cycle [5][6][7][8]. However, it is difficult to determine how changes in aerosol-induced radiation affect GPP changes. Thus, it is necessary to incorporate the influences of aerosols into GPP models to quantify how aerosols affect the terrestrial carbon cycle [6,9].
The enhancement or inhibition effect of aerosols on GPP has been an issue of constant debate [10,11], and the effective quantification of GPP by aerosols has become a scientific problem that needs to be solved urgently. For two reasons, the estimate of the aerosol effects on GPP is controversial. First, moderate aerosol loading regions can increase the diffuse fraction (D f ) of solar radiation. Diffuse light can penetrate canopies more effectively and increase plant photosynthesis via the aerosol diffuse fertilization effect (DFE) [12]. On the contrary, in the dense aerosol loading regions, aerosols have a negative impact on GPP because the diffuse radiation is easily saturated and the direct light is strongly attenuated [13,14]. In addition, the negative effect of high-load aerosols on GPP is also manifested in that aerosol particles are deposited on leaves, thereby blocking stomata, directly affecting plant photosynthesis and related protein expression [10]. For example, a modeling study showed that for a part of the northern coniferous forest, a relatively high aerosol loading (AOD > 2) resulted in a decrease in total radiation and thereby the GPP [4]. Second, the effect of aerosols on GPP may also be limited to the height and leaf area index (LAI) and variations in different ecosystems. Increased diffuse radiation can enhance GPP for forests (i.e., broadleaf forest) with complex canopy structure and high LAI [13,15,16] but not for forests with low canopy (i.e., wetlands) and LAI (i.e., needleleaf forest) [17,18]. However, few studies have currently considered both the regional differences in aerosol loading distribution and the differences in vegetation structure characteristics [19][20][21]. There is an urgent need to systematically study the specific impact of aerosols on GPP from a global perspective.
Existing studies have proved that incorporating the impact of aerosols on GPP in the GPP models is an important and effective way to reduce the uncertainty of GPP estimates [22][23][24]. The transmission of direct light and diffuse light in the canopy can differ due to aerosols [25], and therefore some models need to divide the canopy into sunlit and shaded leaves [26,27]. For example, some current two-leaf models and land surface process models [28,29] consider the influence of aerosols and improve the performance in simulating GPP by calculating photosynthetically active radiation and GPP for shaded and sunlit leaves separately [30]. However, most of the current light use efficient (LUE) models treat the entire canopy as a leaf [31], ignoring the differences in photosynthetic active radiation (PAR) and LUE caused by aerosols among different leaves within the canopy [16,32,33]. Regional and global GPP estimates have substantial uncertainty partly because the influence of aerosol on radiation is not considered [34][35][36]. Yuan et al. [31] found that ignoring the promotion of aerosols on diffuse radiation (increasing the diffuse radiation by 43%) led to an underestimation of the European GPP. Therefore, we emphasize the need to quantify the relationship between aerosols and diffuse radiation, which is a factor that cannot be ignored in future GPP simulation improvements.
In this study, we improved a widely used LUE model, EC-LUE, by including aerosol-radiation effects, and examined the potential impact of aerosols on GPP at a global scale. The EC-LUE is a widely used tool to estimate large-scale real GPPs because of its reasonable performance and few input parameters [37]. Data from the EC-sites and satellites were used to monitor GPP at the global scale and evaluate the impact of aerosols on it. The objectives of this study were: (1) to quantify the relationship between aerosols and diffuse radiation; (2) to determine whether aerosols and diffuse radiation can be included in the GPP model; and (3) to quantify the aerosol impact on GPP at the global scale.

Eddy Covariance Flux Tower Data
Data from the FLUXNET2015 dataset (Available online: http://www.fluxdata.org, accessed on 11 March 2020) explain ecosystem-scale processes such as the exchange of carbon dioxide, water, and energy between the biosphere and the atmosphere. The data processing and site description details can be found on the FLUXNET2015 dataset website (Available online: https://fluxnet.org/data/fluxnet2015-dataset/, accessed on 9 June 2020). In this paper, a total of 70 eddy covariance (EC) sites were selected (Table S1), and each of these sites had more than five years of data that were able to capture interannual variations and long-term changes ( Figure 1). According to the classification scheme developed by the International Geosphere-Biosphere Programme, the land cover classes are divided into 10 types (Table S1)

Eddy Covariance Flux Tower Data
Data from the FLUXNET2015 dataset (Available online: http://www.fluxdat cessed on 11 March 2020) explain ecosystem-scale processes such as the exchang bon dioxide, water, and energy between the biosphere and the atmosphere. The d cessing and site description details can be found on the FLUXNET2015 dataset (Available online: https://fluxnet.org/data/fluxnet2015-dataset/, accessed on 9 Jun In this paper, a total of 70 eddy covariance (EC) sites were selected (Table S1), and these sites had more than five years of data that were able to capture interannu tions and long-term changes ( Figure 1). According to the classification scheme de by the International Geosphere-Biosphere Programme, the land cover classes are into 10 types (Table S1)  The daily GPP variable (GPP_NT_VUT_REF) was estimated using net ecosy change (NEE) methods [38]. The air temperature (Ta), sensible heat (H), and lat (LE) variables were quantified for each site. Using marginal distribution sampling the micrometeorological variables have been gap-filled with reanalysis datasets, w carbon flux measurement data were gap-filled with LE and H. The daily averag were aggregated to the 8-day time step to match the time step of other data.

Remote Sensing Data
Aerosol optical depths (AOD) for 550-nm wavelengths are available from th of Global Modeling and Assimilation (GMAO) of NASA (Available https://disc.gsfc.nasa.gov/datasets/, accessed on 14 July 2021). The data is par NASA atmospheric reanalysis dataset (MERRA-2) for the satellite era [39]. In the project, historical climate analyses are conducted at a range of time scales c weather and climate, and the NASA EOS observations are placed in a climate con The SYN1deg Version 4.1 Cloud and Earth Radiant Energy System (CERES) [40] provided the direct PAR and diffuse PAR data for clear skies for the period fr The daily GPP variable (GPP_NT_VUT_REF) was estimated using net ecosystem exchange (NEE) methods [38]. The air temperature (Ta), sensible heat (H), and latent heat (LE) variables were quantified for each site. Using marginal distribution sampling (MDS), the micrometeorological variables have been gap-filled with reanalysis datasets, while the carbon flux measurement data were gap-filled with LE and H. The daily average values were aggregated to the 8-day time step to match the time step of other data.

Remote Sensing Data
Aerosol optical depths (AOD) for 550-nm wavelengths are available from the Office of Global Modeling and Assimilation (GMAO) of NASA (Available online: https://disc.gsfc. nasa.gov/datasets/, accessed on 14 July 2021). The data is part of the NASA atmospheric reanalysis dataset (MERRA-2) for the satellite era [39]. In the MERRA project, historical climate analyses are conducted at a range of time scales covering weather and climate, and the NASA EOS observations are placed in a climate context. The SYN1deg Version 4.1 Cloud and Earth Radiant Energy System (CERES) product [40] provided the direct PAR and diffuse PAR data for clear skies for the period from 2000 to 2014. The PAR total was calculated as the sum of PAR dir and PAR dif . The SYN1deg Ed4A product is designed to provide top of the atmosphere (TOA) and surface flux data with the highest temporal resolution by incorporating hourly GEO imager data. The CERES PAR data are available online (Available online: https://ceres-tool.larc.nasa.gov/ord-tool/ jsp/SYN1degEd41Selection.jsp, accessed on 3 May 2021).

Description of the AOD-LUE Model
We developed a new LUE model, AOD-LUE, by improving the EC-LUE model [37] by incorporating the relative effects of diffuse and direct radiation on GPP. A total of 35 EC sites encompassing ten vegetation types were selected for calibration, and the remaining 35 EC sites were used for model verification (Figure 1). The EC-LUE model is expressed as: where fPAR is the fraction of incident daily photosynthetically active radiation absorbed by the vegetation canopy (g C m −2 MJ −1 APAR) (%), PAR is the incident daily photosynthetically active radiation (MJ m −2 ), ε max is the maximum light use efficiency without environmental stress (APAR) [2], min is the minimum function of two scalars varying from 0 to 1 and represents the reduction of potential LUE under limiting environmental conditions, and T s and W s , are temperature and water downward regulation scalars, respectively.
To explicitly consider the Df effect of aerosols, we used the D f fitting to divide PAR into PAR dir and PAR dif in the EC-LUE model. At the same time, the maximum LUE parameter (ε) will also be parameterized for PAR dir and PAR dif separately. We use a nonlinear regression equation (NLS) to fit the LUE values of the two types of radiation respectively. The AOD-LUE model is as follows: where LUE dir is the maximum LUE (g C m −2 MJ −1 APAR) affected by PAR dir , LUE dif is the maximum LUE affected by PAR dif . In the overall model representation, through the NLS equation, LUE dir and LUE dif were estimated to be 1.68 (g C m −2 MJ −1 ) and 3.61 (g C m −2 MJ −1 ) respectively.
In the model, fPAR is approximated as linear function of NDVI: where a and b are empirical constants. According to Sims et al. [41], a and b are set to 1.24 and −0.168, and NDVI is obtained directly from 1 km MODIS data. The influence of aerosol on solar radiation is mainly manifested in two aspects, the enhancement of PAR dif and the weakening of PAR dir . In order to consider the influence of aerosol on the diffusion part of global irradiance, the D f was calculated as: where PAR total is the total photosynthetically active radiation, and PAR total is the total PAR. Aerosol and D f were fitted by the Michaelis-Menten equation: where V max is the maximum rate achieved by the ecosystem when the substrate concentration is saturated. The Michaelis constant K M is the concentration of substrate at which a reaction rate is halved of V max . The concentrations at which V max and K M were optimized in this study were 0.6 and 0.14 respectively [11,42].
Estimated T s based on the equation developed for terrestrial ecosystem model (TEM): where T min (0 • C), T max (40 • C), and T opt (23.5 • C) correspond to the minimum, maximum, and optimum daily air temperatures ( • C) for photosynthesis. T s = 0 if T min or T max fall below or above the minimum air temperature [37]. The evaporative fraction (W s ) formula is as follows: where LE is EC-measured latent heat flux (W m −2 ), and H is sensible heat flux (W m −2 ). The determination coefficient (R 2 ), RMSE, and SD of the predicted value and the observation value are calculated to evaluate the performance of the model at the site level [43,44].

Data Analysis and Prediction of GPP
This study used ArcGIS 10.4 to calculate the annual average values of AOD, PAR dir , PAR dif and PAR total from 2000 to 2014, and plot the global average. A linear trend analysis method was used to study the annual trend of each variable, and then to test the significance of the regression slope.
In this research, an orthogonal regression method was used to compare each spatial data product with the corresponding site-level measurement value. Orthogonal regression can consider the errors of independent variables and dependent variables at the same time [45]. The correction coefficient obtained by orthogonal regression was used to correct the spatial data set for the prediction of GPP globally. The correction coefficients were derived by regressing a flux tower variable value against a spatially derived variable value based on slope and intercept. For the global GPP simulation, the AOD-LUE model was driven by the MERRA-2 reanalysis dataset. The 1km NDVI data we obtained are aggregated these data to 0.1 degree resolution and monthly time step.
In addition, we used the AOD-LUE model to estimate GPP with the 0.1 • spatial resolution and monthly temporal resolution for the globe during 2000-2019 ( Figure S1a). Using linear regression analysis in Python3, we determined the long-term trend in GPP from 2000 to 2019 ( Figure S1b).

Aerosols Contribution to Changes in GPP
In order to represent different biomes on the global spatial scale, the global potential natural vegetation (PNV) map at 1 km spatial resolution (Available online: https://doi. org/10.7910/DVN/QQHCIK, accessed on 18 February 2021) was used for the prediction of GPP. The PNV map includes 10 vegetation types. We used PVN vegetation classification because it assumes an undisturbed natural vegetation state and can well reflect the zonal characteristics of vegetation [46]. The aerosol-radiation effect has a higher latitudinal direction in the influence of GPP, and the use of PVN can better show the influence of aerosol on GPP in different latitudes vegetation types.
To isolate the contribution of aerosols to GPP, we also simulated GPP with a global spatial resolution of 0.1 • and a monthly temporal resolution from 2000 to 2019 using the EC-LUE model. The simulation based on the AOD-LUE model considered aerosol radiation effects and was referred to as S0 hereafter, while the simulation based on the EC-LUE model ignored aerosol radiation effects and was referred to as S1 hereafter. Two metrics: ∆GPP (Equation (8)) and ∆GPP (%) (Equation (9)) were used to measure the differences in GPP between the S0 and S1 simulations:

Spatio-Temporal Variations of AOD, PAR dif , PAR dir and PAR total
The global mean values of AOD, PAR dir , PAR dif , and PAR total had large spatial differences from 2000 to 2014 ( Figure 2). The high values of AOD and PAR dif (>0.8 and 50W m −2 , respectively) were found in regions with rapidly developing economies such as eastern China and northern India, mainly due to man-made pollution (Figure 2a,c); the moderate values were found in desert areas (e.g., the Arabian Peninsula and North Africa) likely due to frequent sandstorms; the low values (<0.2 and 20W m −2 , respectively) were found in the Qinghai-Tibet Plateau and Andes Mountains. The PAR dir distribution patterns (Figure 2b) showed the opposite pattern to those of AOD and PAR dif on the Plateau, in the Himalayas, and in the Andes. In addition, PAR total and PAR dir values had similar spatial patterns (Figure 2b,d), and the influence of PAR dir on PAR total was greater than that of PAR dif .

Spatio-Temporal Variations of AOD, PARdif, PARdir and PARtotal
The global mean values of AOD, PARdir, PARdif, and PARtotal had large spatial diff ences from 2000 to 2014 ( Figure 2). The high values of AOD and PARdif (>0.8 and 50W m respectively) were found in regions with rapidly developing economies such as easte China and northern India, mainly due to man-made pollution (Figure 2a,c); the moder values were found in desert areas (e.g., the Arabian Peninsula and North Africa) lik due to frequent sandstorms; the low values (<0.2 and 20W m −2 , respectively) were fou in the Qinghai-Tibet Plateau and Andes Mountains. The PARdir distribution patterns (F ure 2b) showed the opposite pattern to those of AOD and PARdif on the Plateau, in Himalayas, and in the Andes. In addition, PARtotal and PARdir values had similar spat patterns (Figure 2b,d), and the influence of PARdir on PARtotal was greater than that PARdif.  (Figure 3b,d). The value for PARdif show a significant growth trend, while that for AOD dropped noticeably in tropical regio (e.g., Indonesia and Congo Basin) (Figure 3a,c). Values of PARdif, PARdir and PARtotal had increasing trends in most parts of the Southern Hemisphere (Figure 3b,d). At The trends of AOD, PAR dir , PAR dif , and PAR total mean values over the period 2000-2014 varied across regions ( Figure 3). For example, in India and southeastern China, values for AOD and PAR dif exhibited significant growth trends (Figure 3a,c), while those for PAR dir and PAR total showed significant declines (Figure 3b,d). The value for PAR dif showed a significant growth trend, while that for AOD dropped noticeably in tropical regions (e.g., Indonesia and Congo Basin) (Figure 3a,c). Values of PAR dif , PAR dir and PAR total all had increasing trends in most parts of the Southern Hemisphere (Figure 3b,d). At the global scale, the annual averages of AOD and PAR dif showed upward trends (0.03 Wm −2 and 0.1 Wm −2 ) (Figure 3e,g), while the annual averages of PAR dir (−0.3 Wm −2 ) and PAR total displayed downward trends (Figure 3f

High Correlations between AOD and Df
The AOD had a high correlation with the diffusion fraction (Df) (Figure 4; R 2 = p < 0.01), and the influence of AOD on Df can be expressed by the distribution of the M aelis-Menten equation (Equation (5): Df = Vmax × AOD/(Km + AOD), where AOD is range from 0.1 to 1). Among them, Vmax represents the threshold of Df that could b fected by AOD (Figure 4), that is to say, AOD will not have an impact on Df once D ceeds this value (0.6). In addition, Km represents the AOD value (0.14) when Df is hal (i.e., 0.3) (Figure 4). According to this equation, we can know the increase velocity o fusivity under different aerosol loading. When the AOD loading was low (from 0 to Df increases rapidly (from 0 to 0.3), which was a first-order reaction. As the value of increases (from 0.14 to 1), the growth rate of Df starts to slow down (from 0.3 to 0.6), w was a second-order reaction. Overall, most of the AOD values were low, indicating a scatter at low AOD.  -d) show the significance level of the trends of AOD, PAR dif, PAR dir , and PAR total , respectively, with the red area indicating significant trends (p-value < 0.05) and the black area indicating insignificant trends. The temporal variation of annual AOD, PAR dir , PAR dif and PAR total that were spatially averaged over the globe are illustrated in (e-h), respectively; the red lines are the fitting lines.

High Correlations between AOD and D f
The AOD had a high correlation with the diffusion fraction (D f ) (Figure 4; R 2 = 0.57, p < 0.01), and the influence of AOD on D f can be expressed by the distribution of the Michaelis-Menten equation (Equation (5): D f = V max × AOD/(K m + AOD), where AOD is set to range from 0.1 to 1). Among them, V max represents the threshold of D f that could be affected by AOD (Figure 4), that is to say, AOD will not have an impact on D f once D f exceeds this value (0.6). In addition, K m represents the AOD value (0.14) when D f is half V max (i.e., 0.3) (Figure 4). According to this equation, we can know the increase velocity of diffusivity under different aerosol loading. When the AOD loading was low (from 0 to 0.14), D f increases rapidly (from 0 to 0.3), which was a first-order reaction. As the value of AOD increases (from 0.14 to 1), the growth rate of D f starts to slow down (from 0.3 to 0.6),

Improvement of GPP Estimation by Considering the Aerosol-Radiation Effect
Compared with the EC-LUE model, the improved GPP model that considered ae sol-radiation effects (i.e., the AOD-LUE model) can more accurately estimate GPP (Figu 5a). The AOD-LUE model had higher R² (0.68) and slope (0.93) than the EC-LUE mo (R 2 = 0.6 and slope = 0.89). The data points of the AOD-LUE model values were mo concentrated in the vicinity of the 1:1 line than those of the EC-LUE model, which in cated that the GPP values estimated by the AOD-LUE model were closer to the true valu (i.e., flux tower GPP). Additionally, the simulated value of the AOD-LUE model d creased when the tower's GPP was less than 7.5, and vice versa. The degree of impro ment of the AOD-LUE model differed across EC sites (Figure 5b). To be more specific, AOD-LUE model was better optimized at the EC sites where the R 2 for the EC-LUE mo was lower (Figure 5b). The AOD-LUE model can be greatly improved (is about 40-120 when the R 2 value of the EC-LUE model is in the range of 0-0.4., but when the R 2 value the EC-LUE model was high (i.e., R 2 greater than 0.6), the R 2 value of the AOD-LUE mo was not greatly improved (less than 10%).

Improvement of GPP Estimation by Considering the Aerosol-Radiation Effect
Compared with the EC-LUE model, the improved GPP model that considered aerosolradiation effects (i.e., the AOD-LUE model) can more accurately estimate GPP (Figure 5a). The AOD-LUE model had higher R 2 (0.68) and slope (0.93) than the EC-LUE model (R 2 = 0.6 and slope = 0.89). The data points of the AOD-LUE model values were more concentrated in the vicinity of the 1:1 line than those of the EC-LUE model, which indicated that the GPP values estimated by the AOD-LUE model were closer to the true values (i.e., flux tower GPP). Additionally, the simulated value of the AOD-LUE model decreased when the tower's GPP was less than 7.5, and vice versa. The degree of improvement of the AOD-LUE model differed across EC sites (Figure 5b). To be more specific, the AOD-LUE model was better optimized at the EC sites where the R 2 for the EC-LUE model was lower (Figure 5b). The AOD-LUE model can be greatly improved (is about 40-120%) when the R 2 value of the EC-LUE model is in the range of 0-0.4., but when the R 2 value of the EC-LUE model was high (i.e., R 2 greater than 0.6), the R 2 value of the AOD-LUE model was not greatly improved (less than 10%).
For each vegetation type, the simulation performance of the AOD-LUE model outperformed the EC-LUE model ( Table 1). As compared with the EC-LUE model, the AOD-LUE model produced higher R 2 and less SD across all vegetation types. The improvement was more obvious in EBF and CSH, in which the EC-LUE model has low R 2 (0. 35 The results show that the AOD-LUE model has higher potential in some poorly simulated vegetation types. (i.e., flux tower GPP). Additionally, the simulated value of the AOD-LUE mo creased when the tower's GPP was less than 7.5, and vice versa. The degree of im ment of the AOD-LUE model differed across EC sites (Figure 5b). To be more spe AOD-LUE model was better optimized at the EC sites where the R 2 for the EC-LU was lower (Figure 5b). The AOD-LUE model can be greatly improved (is about 4 when the R 2 value of the EC-LUE model is in the range of 0-0.4., but when the R 2 the EC-LUE model was high (i.e., R 2 greater than 0.6), the R 2 value of the AOD-LU was not greatly improved (less than 10%).

Impact of Aerosols on GPP
The impact of aerosols on GPP (∆GPP = ±400 g C m −2 yr −1 ) showed spatial heterogeneity on a global scale (Figure 6a). Aerosols had positive effects on GPP in many areas across the globe, such as northern China, Central Asia, India, Europe, tropical and southern Africa, and Mexico, and the positive impact of ∆GPP could be up to 300−400 g C m −2 yr −1 (6-8%) (Figure 6a,c). The negative effects of aerosols on GPP were mainly concentrated in tropical Asia, western Europe, Amazon, and coastal areas in some regions (Figure 6a). From the perspective of latitude distribution, aerosols had a large weakening effect on ∆GPP (200 g C m −2 yr −1 ) in tropical and high latitude (220 g C m −2 yr −1 ) regions but had a large enhancement effect for most subtropical (110 g C m −2 yr −1 ) and temperate (80 g C m −2 yr −1 ) regions (Figure 6b,d). From the percentage change in the spatial distribution of GPP caused by the aerosol, it was found that there was a large substantial increase in GPP in Central Asia, the Sahara Desert, and northwestern China. The GPP values in the Amazon rainforest, Southeast Asia, and boreal forest areas declined significantly due to the impact of the aerosol (Figure 6c). cluding tropical vegetation and alpine vegetation (Figure 7). Specifically, there was a ne ative correlation between ΔGPP and aerosols in Evergreen Broadleaf Forest (EB (−6.45%), Evergreen Broadleaf Forest (EBF) (−3.97%), Evergreen Needleleaf Forest (EN (−3.36%) and Wetlands (WET) (−4.67%). Aerosol positively influenced ΔGPP in the main Grassland (GRA) (3.19%) and Deciduous Broadleaf Forest (DBF) (2.63%), graminoid an Open Shrub (OSH) (2.38%), Crop (CRO) (2.19%) and Woody Savanna (WSA) (1.49%).

Model Improvement by Incorporating the Effects of Diffuse Radiation
The traditional ecosystem models typically do not consider the impact of diffuse PAR [47], and have a certain deviation in simulated GPP under clear and cloudy conditions [48]. However, the changes in diffuse radiation under natural conditions are often caused by changes in aerosols, and our results found that aerosols can significantly increase diffuse PAR. Gu et al. [30] hypothesized that aerosols would increase NPP due to the increased quantum yield of diffuse light. However, the effects of diffuse light on GPP have not been quantified globally. In addition, we also found that when there was a higher D f , aerosols weakened plant photosynthesis by reduced direct PAR. Therefore, compared with previous modeling studies that only considered the total incident PAR [49,50], our research results revealed that the relative effects of diffuse PAR and direct PAR on plant growth should be considered in GPP models.
Most LUE models define GPP as the product of PAR and LUE absorbed by the vegetation canopy [2]. When PAR is affected by aerosol and changes, maximum LUE will also change [51]. In the estimation of GPP based on remote sensing, the GPP estimation caused by the structural changes of LUE can also be considered as a key component of the total error budget [21,52,53]. The LUE models like EC-LUE typically use a constant LUE independent of the biological community [48], and assume that LUE is independent of the diffuse PAR [54], which leads to the underestimation of GPP on cloudy days with more diffuse radiation.
In this research, through the AOD-LUE model, we quantified the maximum LUE parameter of PAR dif and PAR dir (3.61 and 1.68, respectively). The maximum LUE of PAR dif is much higher than that of PAR dir mainly because canopy photosynthesis tends to use light more effectively under diffuse light than direct light [55]. Most of the modeling and observational studies also proved that LUE increases when the D f increases [56,57]. The regression relationship between the GPP estimated by the EC tower data and the GPP simulated by the AOD-LUE model has a slope of 0.68, which is higher than the slope for the original EC-LUE model (Figure 5a). This means that the model with the effects of aerosols added further shows the spatial heterogeneity of global GPP [27,58].

Differences in the Impact of Aerosol Diffuse Fertilization Effect on GPP in Different Regions
According to our modeling, aerosol-induced changes in PAR have a strong impact on GPP depending on aerosol loading and cloud thickness. Moderate aerosol increases the photosynthesis of plants by increasing the D f of total PAR, while the further increase in aerosol loading can have the opposite effect due to the strong attenuation of total PAR. Moderate aerosol loading would increase plant photosynthesis and the total partition of PAR to D f , while further increases in aerosol loading could have an opposite effect due to a strong attenuation of total PAR. The key mechanism of aerosols affecting GPP is its effect on total PAR and partitioning of total PAR into direct and diffuse fractions. Diffuse PAR is more effective than direct PAR in penetrating vegetation canopy [11,30,42] and therefore has a higher LUE coefficient (3.61 g C m −2 MJ −1 ) than direct PAR (1.68 g C m −2 MJ −1 ). In addition, the aerosol diffuse fertilization effect (DFE) was very closely related to cloud cover [12,59], which caused the D f of the area to easily reach saturation and reduce the availability of light [60]. Therefore, the overall impact of aerosol DFE on GPP may be positive in areas with low cloud cover and aerosol loading, and may be negative in areas with higher cloud cover and aerosol loading.
The negative impact of aerosols on GPP is the greatest in the tropics (∆GPP < −4%). This is mainly because tropical regions with strong evapotranspiration have limited potential for aerosol DFE due to thick cloud cover [12,61], which cannot compensate for the reduction of total PAR and direct PAR (Figure 3b,d). However, our results showed that ∆GPP > 2% is common in subtropical areas with low cloud cover, especially in some arid and semi-arid regions (∆GPP > 3%) (Figure 6d). The cloudiness in these areas is generally low and does not inhibit the aerosol DFE [25,61] likely because the vegetation photosynthesis in these areas does not saturate [8,62] and the change of aerosol to diffuse PAR increases the potential LUE. The mean value of diffuse PAR in these areas was relatively high and exhibited an increasing trend.
However, compared with subtropical regions with moderate aerosols, southeastern China did not exhibit increases in GPP due to aerosol DFE despite a heavy aerosol loading (AOD > 0.5). Previously, Yue and Unger [63] considered the impact of aerosols and ozone pollutants on net primary productivity (NPP) and found that the aerosol DFE would increase NPP by 0.2 Pg C (5%) in eastern China. The reason why our research results are different from them is likely because we have considered that the increase of D f is limited by aerosol loading and the corresponding offsetting effect of direct PAR attenuation [25].
Although aerosols have positive and negative effects on GPP, they are even offset by positive and negative effects in the global effect. What cannot be ignored, is the change of scattered radiation caused by aerosols. Although the relationship between scattered radiation and aerosol coincidence leads to a small change in GPP, it can be compensated by the increase of LUE [10].

Differences in the Impact of Aerosol DFE on Different Vegetation Types
Our results show that aerosol DFE increased GPP for most broad-leaved forests. However, in coniferous forests, the aerosol greatly reduced the total PAR and reduced the ∆GPP by 3.36%. One possible explanation is that the aerosol DEF of coniferous forests is not as obvious as that of broad-leaved forests [64], which may be caused by the sparse canopy and lower leaf area index (LAI) of coniferous forests [12,65]. Another possible reason is that the negative effect by the reduction of total PAR more than offset the DEF. As a result of the increase in the proportion of leaves exposed to moderate light levels, terrestrial ecosystems with high LAI were more sensitive to diffuse radiation [66]. A modeling study by Alton et al. [67] showed that under the condition of an increased diffuse fraction, the actual LUE increased by 33% for a broad-leaved forest with an LAI of 5.05 but by only 6% for a Scottish coniferous forest with an LAI of 2. The simulations by Knohl and Baldocchi [15] using a multilayer canopy model also indicated that the aerosol DFE decreased with a decrease in LAI.
In addition, previous studies have suggested that GPP in some ecosystems with fewer leaves or open canopies, such as grasslands and wetlands, is also less sensitive to diffuse radiation [17,50]. However, in our research, we found that in some drought-tolerant grasslands (such as WSA) with low LAI, the aerosol DFE caused a significant increase in GPP (Figure 7). Jing et al. [29] found that CO 2 absorption increased significantly due to aerosols and thick clouds in grasslands with extremely low LAI (0.37). There is a significant relationship between canopy structure and the aerosol DFE for the multilayer arctic shrub system with low LAI (1.5) [68]. Some previous studies also have suggested that the LAI of temperate coniferous forests and temperate deciduous forests was similar, and LAI alone cannot explain the differences in the response of GPP between plant functional types (PFTs) to the aerosol DFE [15,69]. Therefore, our results support Park et al.'s [11] view that canopy structure may be more important in determining the aerosol DFE than LAI.

Limitations and Future Needs
There are still some uncertainties in our quantitative study of the influence of aerosol on GPP. First, the current flux tower sites measure less diffuse PAR when conducting carbon and water flux observations [29,36]. Although most process-based models can now differentiate between sunlit and shaded leaves [21,36,70], field measurements are still necessary in order to parameterize the effects of diffuse PAR on LUE and photosynthesis. However, most of the diffuse radiation data used nowadays are mainly derived from experience or numerical models [67,71]. Therefore, the lack of field-based diffuse PAR measurements can lead to uncertainty in simulated effects of scattered radiation on GPP by numerical models.
Second, the indirect effects of aerosols on meteorological changes (i.e., precipitation) can also affect GPP. Yue et al. [63] found that aerosol-induced drought strongly reduced NPP, which may be due to aerosol-induced changes in evaporation and precipitation [22,72]. In addition, the indirect effects of aerosols can also change the size and distribution of cloud droplets and cloud albedo [32,73] and increase the depth and number of clouds. Clouds can also disturb the scattered radiation, surface temperature and precipitation and thus have complex effects on the terrestrial carbon cycle. Because the relationship between aerosols and clouds is so complicated, current research cannot completely separate them [74,75]. To resolve these meteorological feedbacks and accompanying mechanisms, future efforts are needed to fully couple the terrestrial biosphere, atmosphere chemistry and climate in earth system models [63].
Third, the input data also affect the simulation of the impact of aerosols on GPP. The effects of spatial resolution and uncertainty of remote sensing data on carbon fluxes cannot be ignored [76,77]. Some uncertainties in the simulation of the AOD-LUE model may be caused by the scale mismatch between the input data sets. The spatial resolution of the AOD data used in this article was 0.5 • × 0.625 • , while the spatial resolution of other simulated flux tower datasets was generally less than 0.03 • × 0.03 • . This kind of uncertainty in GPP simulation caused by data scale mismatch is inevitable [1,53]. For example, researchers Wang et al. [35] found that spatial PAR data explained 57% of the PAR variation detected at EC tower sites [78,79]. To achieve a global level of spatial resolution for GPP modeling on a large scale, the GPP model research center still must improve the quality of remote sensing data.

Conclusions
The relationship between aerosols and diffuse radiation was quantified and a new LUE model developed in this study to explicitly account for the diffuse fertilization effect of aerosols on GPP. By changing the structure of the widely used EC-LUE model, the maximum LUE of diffuse PAR and direct PAR were incorporated and parameterized. In comparison with the original model (EC-LUE), the new model (AOD-LUE) shows fairly good performance across different biomes (R 2 = 0.68, p < 0.01). The model was then used to simulate GPP at the global level from 2000-2014. We then used the AOD-LUE model to simulate GPP at the global scale for the period 2000-2014. With the global-scale GPP estimates, we quantified the effects of aerosol radiation on GPP regionally and globally and for different vegetation types. Although the total effect is low, it still needs attention, because the total effect is offset by both positive and negative aspects, but the positive and negative effects of different regions are very large. Aerosols mainly have a relatively large negative impact on Closed Shrub (CSH) (∆GPP = −6.45%), and a positive impact on arid and semi-arid regions (GRA) (∆GPP = 3.19%) in the subtropical zone. Our results showed that the effects of aerosol radiation have severely affected the global GPP. Separating PAR to diffuse and direct components and incorporating their relative effects to LUE models are effective for improving the accuracy of the GPP simulation at regional to global scales.