Modeling the Effects of Global and Diffuse Radiation on Terrestrial Gross Primary Productivity in China Based on a Two-Leaf Light Use Efficiency Model

Solar radiation significantly affects terrestrial gross primary productivity (GPP). However, the relationship between GPP and solar radiation is nonlinear because it is affected by diffuse radiation. Solar radiation has undergone a shift from darker to brighter values over the past 30 years in China. However, the effects on GPP of variation in solar radiation because of changes in diffuse radiation are unclear. In this study, national global radiation in conjunction with other meteorological data and remotely sensed data were used as input into a two-leaf light use efficiency model (TL-LUE) that simulated GPP separately for sunlit and shaded leaves for the period from 1981 to 2012. The results showed that the nationwide annual global radiation experienced a significant reduction (2.18 MJ m−2 y−1; p < 0.05) from 1981 to 2012, decreasing by 1.3% over this 32-year interval. However, the nationwide annual diffuse radiation increased significantly (p < 0.05). The reduction in global radiation from 1981 to 2012 decreased the average annual GPP of terrestrial ecosystems in China by 0.09 Pg C y−1, whereas the gain in diffuse radiation from 1981 to 2012 increased the average annual GPP in China by about 50%. Therefore, the increase in canopy light use efficiency under higher diffuse radiation only partially offsets the loss of GPP caused by lower global radiation.


Introduction
Carbon sequestration is extremely important for global change studies [1]. Terrestrial gross primary productivity (GPP), which indicates the ability of vegetation to use light energy through fixation of carbon dioxide during photosynthesis, is a primary component of the terrestrial carbon cycle [2]. Solar radiation substantially affects the ecosystem GPP [3].
Global radiation (R g ) includes both direct and diffuse components. Compared with direct radiation, diffuse radiation (R d ) is more readily distributed evenly in the canopy and shows higher photosynthetic saturation. Therefore, higher R d increases the canopy light use efficiency (LUE) and The Rg at the site scale was interpolated to the national scale by the following steps. First, the Ångström model (Equation (1)) [29] was calibrated using data from 103 meteorological stations recording both Rg and sunshine duration measurements during 1981-2012. Second, to compensate for the lack of global radiation measurements, Rg for the remaining 651 meteorological stations was calculated by using sunshine duration and the calibrated Ångström model in step 1. Third, data from 12 sites with Rg measurements in 2014 were used to validate the calibrated Ångström model (see Supplementary Materials, Figure S2). Finally, using radiation data at the 754 sites, Rg was interpolated to the national scale with an inverse distance weighted method. The R g at the site scale was interpolated to the national scale by the following steps. First, the Ångström model (Equation (1)) [29] was calibrated using data from 103 meteorological stations recording both R g and sunshine duration measurements during 1981-2012. Second, to compensate for the lack of global radiation measurements, R g for the remaining 651 meteorological stations was calculated by using sunshine duration and the calibrated Ångström model in step 1. Third, data from 12 sites with R g measurements in 2014 were used to validate the calibrated Ångström model (see Supplementary Materials Figure S2). Finally, using radiation data at the 754 sites, R g was interpolated to the national scale with an inverse distance weighted method.
Considering the heterogeneity of terrain and climate in China, the country was divided into eight regions: Northeast China (R1), Inner Mongolia (R2), Northwest China (R3), North China (R4), Central China (R5), South China (R6), Southwest China (R7), and the Qinghai-Tibet Plateau (R8). The R g estimation model was fitted for each region separately to convert sunshine duration into R g , which was comparable with the method of Ren et al. [25,30]. The Ångström model [29] was selected to estimate R g because of the advantages of a simple structure, clear physical meaning, stable model coefficients, and good fitting results in addition to its widespread use [31]. This model was as follows: where R g is the global radiation reaching the surface, R 0 is the daily extraterrestrial radiation, n is the actual daily sunshine duration observed at a given site, N is the potential daily sunshine duration, and a and b are coefficients. Both R 0 and N can be calculated based on latitude and date [32]. Root mean square error (RMSE) values of the radiation derived using Ångström models for each of the eight regions of China are shown in Table 1. The R d at the site scale was interpolated to the national scale based on the following steps. First, data from 81 meteorological stations recording both R g and R d during 1981-2012 were used to calibrate the logistic model (Equation (2)) [27,30]. Second, to compensate for the lack of R d measurements, R g estimated with the Ångström model was used to calculate R d for the remaining 673 meteorological stations using the calibrated Equation (2) in step 1. Third, data from 12 sites with R g measurements in 2014 were used to validate the calibrated Ångström model (see Supplementary Materials, Figure S3). Finally, R d was interpolated to the national scale using R d at the 754 sites.
A simple and accurate logistic model [27,30] was used to estimate R d as follows: where R d is the diffuse radiation, and c and d are the nonlinear regression coefficients. Table 2 shows the RMSE values of estimated R d derived from logistic models for the eight regions of China. The coefficient of determination (R 2 ) for R g and R d measurements and estimations ranged from 0.66 to 0.79 (Table 1) and 0.63 to 0.81 (Table 2), which were comparable with the results of Singh [33] and Ren et al. [30]. Validation (see Supplementary Materials, Figures S2 and S3) indicated that the R 2 for R g measurements and estimations was 0.85-0.96 and that for R d was 0.55-0.75 at validation sites. Estimated global and R d in this study showed the similar values to those of Ren et al. [30] (see Supplementary Materials, Figure S4).

LAI Data
The LAI at a 0.0727 • grid for 1981-2012 was generated using the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Advanced Very High Resolution Radiometer (AVHRR) data. The dataset was generated by Liu et al. [34] and has been widely used in studies of the carbon cycle and vegetation response to climate change of terrestrial ecosystems [35][36][37][38][39]. GLOBMAP LAI V3 during 1981 to 2012 was generated by fusing LAI inverted from MODIS reflectance data with Global Inventory Modelling and Mapping Studies (GIMMS) AVHRR Normalized Difference Vegetation Index (NDVI) data. The LAI from 2001 to 2012 was first derived from the MODIS land surface reflectance data (MOD09A1) reflectance based on the GLOBCARBON LAI algorithm. For the fusion of MODIS and AVHRR remotely sensed data, the relationships between GIMMS NDVI and MODIS LAI were established pixel by pixel over a period (2001)(2002)(2003)(2004)(2005)(2006)) that they overlap. The AVHRR LAI from 1981 to 2000 was then generated using these relationships.

Land-cover Data
Given the lack of land-cover data before 2000, we used land-cover data in 2001 for the entire period simulations. The MODIS land-cover data at 500-m spatial resolution (MCD12Q1) for 2001 was resampled to 0.0727 • to drive the TL-LUE model during 1981 and 2012. The International Geosphere-Biosphere Programme classification system was used in this study, including categories of evergreen needle-leaf forest (ENF), deciduous coniferous forest (DNF), evergreen broad-leaved forest (EBF), deciduous broad-leaved forest (DBF), evergreen broad-leaved forest (ENF), deciduous coniferous forest (DNF), evergreen broad-leaved forest (EBF), mixed forest (MF), sparse shrubland (OS), dense shrubland (CS), grassland (GRASS), farmland (CROP), and other surface cover types (NOV).

GPP Data
Flux data from 13 available flux tower sites of the ChinaFLUX network were compared with the simulated GPP. These sites included three broad-leaved forest (BF), two needle-leaved forest (NF), one MF, one shrubland (Shrub), four GRASS, and two CROP sites, providing a total of 100 site years (Table 3).

TL-LUE Model
The TL-LUE model [12] is based on MODIS GPP [40] and the Boreal Ecosystem Productivity Simulator (BEPS) model [41]. Model validation was performed at both the site [42,43] and the regional scale [44]. Additional information on the model was introduced in [12,42,43]. The model parameters used for the TL-LUE models followed Zhou et al. [43].
GPP was calculated with the following Equation (3): where ε msu and ε msh are the maximum LUE of sunlit and shaded leaves, respectively; APAR su and APAR sh are the absorbed PAR of sunlit and shaded leaves, respectively; and f (VPD) and f (T amin ) are stress factors of Vapour Pressure Deficiency (VPD) and air temperature, respectively. APAR su and APAR sh were calculated with the following Equations (4) and (5): where α is albedo; PAR dif_u is diffuse radiation under the canopy; C is multiple scattered radiation; β is the leaf angle, which was considered to be 60 • ; θ is the solar zenith angle; PAR dif and PAR dir are incoming direct and diffuse radiation, respectively (Chen et al., 1999). LAI su and LAIsh are the LAI of sunlit and shaded leaves and was calculated with Equations (6) and (7): LAI sh = LAI−LAI su (7) On the basis of Equations (3)-(5), GPP derived from diffuse radiation can be calculated with the following Equation (8): × LAI msh × f (VPD) × g(Ta min) (8) where GPP dif and GPP dir are GPP contributed by diffuse and direct radiation, respectively.

Modeling Scenarios
To quantitatively assess the impacts of variation in R g and diffuse radiation fractions on GPP, three simulations were performed (see Supplementary Materials, Table S3).
(i) Scenario 1 (Scenario): Both historical meteorological data with temporal variations and LAI were used to simulate GPP for the period 1981-2012.
(ii) Scenario 2 (Scenario_R g ): Temporally constant R g data, other historical meteorological data with temporal variations, and LAI were used to simulate GPP for the period 1981-2012. The differences between simulated GPP values in S1 and S2 reflect the effects of change in R g on GPP simulation.
Based on the assumption that the annual R g remained unchanged from 1981 to 2012, and the variation tendency of daily R g variation to R g averaged values during 1981-2012 was consistent, daily radiation for the period 1982-2012 was calibrated using Equation (9): where i and j are the line and row of any given pixel, and R k (i,j)|m is the daily R g of pixel (i, j) on the kth day of the mth year.
(iii) Scenario 3 (Scenario_R d ): Temporally constant R d , other historical meteorological data with temporal variations, and LAI were used to simulate GPP for the period 1981-2012. The difference between simulated GPP values in S1 and S3 reflect the effects of R d on GPP simulations.
Based on the assumption that the annual diffuse radiation fraction remained constant from 1981 to 2012, and the variation tendency of daily R d variation to R d averaged values during 1981-2012 was consistent, the daily diffuse radiation fraction for the period 1982-2012 was calibrated using Equation (10): where i and j are the row and column of a given pixel, and f k di f (i,j)|m is the diffuse radiation fraction in pixel (i, j) for the kth day of the mth year. The diffuse radiation fraction was calculated using Equation (11): where f dif is the diffuse radiation fraction, R d and R g are diffuse and global radiation, respectively.

Spatiotemporal Distributions of Annual Radiation in China
The spatial distributions of the average annual and standard deviation of R g , R d , and diffuse radiation fraction in China over the period 1981-2012 are shown in Figure 2. The R g (Figure 2a) was unevenly distributed across China, with a range of 3255-7154 MJ m −2 y −1 and a national average of 5258 MJ m −2 y −1 (Figure 2a). The standard deviation of R g ranged from 65 to 359 MJ m −2 y −1 (Figure 2b). Generally, R g values on the western plateau were higher than those on the eastern plain. The highest R g values were recorded on the southern part of the Qinghai-Tibet Plateau (>6500 MJ m −2 y −1 ) because of the thin air density in this region. The northeastern part of the Sichuan Basin and the Yunnan-Guizhou Plateau had the lowest annual R g values (<3800 MJ m −2 y −1 ). The Rg exhibited broad variability during 1981-2012 but generally showed a significant (p < 0.05) negative trend of 2.18 MJ m −2 y −1 , yielding a rate of decrease of 1.3% y −1 (Figure 3a). Specifically, Rg decreased markedly and then increased slightly during the interval 1980-1990. During the interval 1991-1993, Rg declined rapidly, possibly because of the volcanic eruption of Mount Pinatubo in 1991. This volcanic eruption led to the release of about 20 million tonnes of sulfide gas into the stratosphere, producing a reduction in Rg. The Rg values were gradually restored after 1994 [9,25,30]. Values of Rg also decreased slightly after 1999, consistent with the findings of Che et al. [22] and Wild et al. [17].
The simple linear trend in the average annual Rd in China was positive from 1981 to 2012, although this trend was not significant (Figure 3b). The Rd decreased during the interval 1981-1990 but increased and then decreased rapidly during the interval 1991-1994 because of the eruption of Mount Pinatubo. From the mid-1990s, Rd increased slightly and then decreased again. The R g exhibited broad variability during 1981-2012 but generally showed a significant (p < 0.05) negative trend of 2.18 MJ m −2 y −1 , yielding a rate of decrease of 1.3% y −1 (Figure 3a). Specifically, R g decreased markedly and then increased slightly during the interval 1980-1990. During the interval 1991-1993, R g declined rapidly, possibly because of the volcanic eruption of Mount Pinatubo in 1991. This volcanic eruption led to the release of about 20 million tonnes of sulfide gas into the stratosphere, producing a reduction in R g . The R g values were gradually restored after 1994 [9,25,30]. Values of R g also decreased slightly after 1999, consistent with the findings of Che et al. [22] and Wild et al. [17].
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 21 The annual average diffuse radiation fraction in China increased markedly between 1981 and 2012 by 0.0002 y −1 (p < 0.05) (Figure 3c). The spatial distribution of the variation in annual Rg, Rd, and diffuse radiation fraction in China over the period 1981-2012 is shown in Figure 4. Decreases in Rg mainly occurred in North China, Hainan, and on the Qinghai-Tibet Plateau (Figure 4a), which is consistent with the findings of Tang et al. [45] and Lin et al. [46]. Conversely, Rg increased only over a small part of China during the The simple linear trend in the average annual R d in China was positive from 1981 to 2012, although this trend was not significant (Figure 3b). The R d decreased during the interval 1981-1990 but increased and then decreased rapidly during the interval 1991-1994 because of the eruption of Mount Pinatubo. From the mid-1990s, R d increased slightly and then decreased again.
The annual average diffuse radiation fraction in China increased markedly between 1981 and 2012 by 0.0002 y −1 (p < 0.05) (Figure 3c).
The spatial distribution of the variation in annual R g , R d , and diffuse radiation fraction in China over the period 1981-2012 is shown in Figure 4. Decreases in R g mainly occurred in North China, Hainan, and on the Qinghai-Tibet Plateau (Figure 4a), which is consistent with the findings of Tang et al. [45] and Lin et al. [46]. Conversely, R g increased only over a small part of China during the period 1981-2012. Generally, the changes in R g over most parts of China were not significant, predominantly because China experienced "dimming" prior to 1990, whereas R g increased slightly over large parts of China post-1990 [23].
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 21 Increases in Rd predominantly occurred over the Qinghai-Tibet Plateau, Shanghai-Hangzhou-Wuhan region, Guangdong, and Guangxi and its surrounding areas (Figure 4b). Decreases in Rd occurred in the northwestern part of Xinjiang, western parts of Inner Mongolia, and Shenyang. The national average diffuse radiation fraction increased by 0.2% decade −1 during the period 1981-2012 ( Figure 4c). Initially, the diffuse radiation fraction increased and then decreased during the 1980s. In 1993, after the eruption of Mount Pinatubo, the diffuse radiation fraction attained its peak during the study period, and then quickly returned to its average value.

Spatial Distribution of GPP in China over the Period 1981-2012
At the site scale, comparison of the simulated GPP values derived from the TL-LUE model (GPPTL) versus site flux data-derived GPP showed good agreement (R 2 = 0.82; N = 100; p < 0.0001) ( Figure 5). At the national scale, the average annual GPP value in China was estimated to be 7.64 Pg C y −1 , which lies within the range of values reported by previous studies (Table 4). Li et al. [47] estimated the average annual GPP from 2001 to 2010 as 6.04 Pg C y −1 , and the annual mean GPP value based on MODIS for China for the same period was estimated as 5.47 Pg C y −1 . The GPP estimates based on relationships with meteorological data were 6.06 and 7.78 Pg C y −1 during 2001-2010 with two spatial resolutions [48]. Increases in R d predominantly occurred over the Qinghai-Tibet Plateau, Shanghai-Hangzhou-Wuhan region, Guangdong, and Guangxi and its surrounding areas (Figure 4b). Decreases in R d occurred in the northwestern part of Xinjiang, western parts of Inner Mongolia, and Shenyang. The national average diffuse radiation fraction increased by 0.2% decade −1 during the period 1981-2012 (Figure 4c). Initially, the diffuse radiation fraction increased and then decreased during the 1980s. In 1993, after the eruption of Mount Pinatubo, the diffuse radiation fraction attained its peak during the study period, and then quickly returned to its average value.

Spatial Distribution of GPP in China over the Period 1981-2012
At the site scale, comparison of the simulated GPP values derived from the TL-LUE model (GPPTL) versus site flux data-derived GPP showed good agreement (R 2 = 0.82; N = 100; p < 0.0001) (Figure 5). At the national scale, the average annual GPP value in China was estimated to be 7.64 Pg C y −1 , which lies within the range of values reported by previous studies (Table 4). Li et al. [47] estimated the average annual GPP from 2001 to 2010 as 6.04 Pg C y −1 , and the annual mean GPP value based on MODIS for China for the same period was estimated as 5.47 Pg C y −1 . The GPP estimates based on relationships with meteorological data were 6.06 and 7.78 Pg C y −1 during 2001-2010 with two spatial resolutions [48]. The multi-year mean GPP in China for the period 1981-2012 is shown in Fig. 6a. In the southeastern and northeastern parts of China, the GPP was higher than 1000 g C m −2 y −1 . On the southeastern coast, Yunnan-Guizhou Plateau, and the southern part of the Qinghai-Tibet Plateau, the GPP was higher than 2500 g C m −2 y −1 . Low values of GPP occurred over desert areas of China and the remaining parts of the Qinghai-Tibet Plateau. Given the scarcity of vegetation in these areas, their annual GPP values were less than 50 g C m −2 y −1 . The northern Tianshan Mountains showed high GPP values of up to 1000-2000 g C m −2 y −1 . Inner Mongolia and the southeastern part of the Qinghai-Tibet Plateau formed part of a transition zone to lower values, showing annual GPP values of 200-1000 g C m −2 y −1 . Therefore, the TL-LUE model's simulation of the average annual distribution of GPP was reasonable and was largely consistent with previous simulations [44,47]. The standard deviation of GPP during 1981-2012 showed a national average of 75 g C m −2 y −1 with a range of 0-838 g C m −2 y −1 (Figure 6b).  Therefore, the TL-LUE model's simulation of the average annual distribution of GPP was reasonable and was largely consistent with previous simulations [44,47]. The standard deviation of GPP during 1981-2012 showed a national average of 75 g C m −2 y −1 with a range of 0-838 g C m −2 y −1 (Figure 6b). The multi-year mean GPP in China for the period 1981-2012 is shown in Fig. 6a. In the southeastern and northeastern parts of China, the GPP was higher than 1000 g C m −2 y −1 . On the southeastern coast, Yunnan-Guizhou Plateau, and the southern part of the Qinghai-Tibet Plateau, the GPP was higher than 2500 g C m −2 y −1 . Low values of GPP occurred over desert areas of China and the remaining parts of the Qinghai-Tibet Plateau. Given the scarcity of vegetation in these areas, their annual GPP values were less than 50 g C m −2 y −1 . The northern Tianshan Mountains showed high GPP values of up to 1000-2000 g C m −2 y −1 . Inner Mongolia and the southeastern part of the Qinghai-Tibet Plateau formed part of a transition zone to lower values, showing annual GPP values of 200-1000 g C m −2 y −1 . Therefore, the TL-LUE model's simulation of the average annual distribution of GPP was reasonable and was largely consistent with previous simulations [44,47]. The standard deviation of GPP during 1981-2012 showed a national average of 75 g C m −2 y −1 with a range of 0-838 g C m −2 y −1 (Figure 6b).   Table 4. Comparison between gross primary productivity (GPP) of major vegetation types and total annual GPP estimated using TL-LUE and other models (ENF means evergreen needle-leaved forest; EBF means evergreen broad-leaved forest; DNF means deciduous needle-leaved forest; DBF means deciduous broad-leaved forest).  [48] Note: The units for mean GPP and national total GPP are g C m −2 y −1 and Pg C y −1 , respectively. a Only data for net primary productivity (NPP) were reported (NPP/GPP = 0.5)

Spatial Distribution of Diffuse Radiation Contribution to Terrestrial Ecosystem GPP
The spatial distribution of GPP driven by R d was consistent with that of GPP driven by R g (Figure 7a-b). GPP driven by Rd ranged from 0 to 4348 g C m −2 y −1 with a national average of 586 g C m −2 y −1 and a national average standard deviation of 56 g C m −2 y −1 . Figure 7c,d show the proportion of GPP fixed by R d to the total GPP and its standard deviation. The proportion of GPP fixed by R d to the total GPP ranged from 0.43 to 0.88 with a national average of 0.71, with very low standard deviation ranging from 0.0007 to 0.12. The highest proportion was in the northern Sichuan Basin and the northern parts of the Yunnan-Guizhou Plateau, where the diffuse radiation fraction was highest. Conversely, in most parts of northwestern and northeastern China, the diffuse radiation fraction was less than 0.5 (Figure 2e), where the proportion of GPP contributed by R d to GPP derived from R g was greater than 0.5, reflecting the higher LUE of R d .

Spatial Distribution of Diffuse Radiation Contribution to Terrestrial Ecosystem GPP
The spatial distribution of GPP driven by Rd was consistent with that of GPP driven by Rg (Figure  7a-b). GPP driven by Rd ranged from 0 to 4348 g C m −2 y −1 with a national average of 586 g C m −2 y −1 and a national average standard deviation of 56 g C m −2 y −1 . Figure 7c,d show the proportion of GPP fixed by Rd to the total GPP and its standard deviation. The proportion of GPP fixed by Rd to the total GPP ranged from 0.43 to 0.88 with a national average of 0.71, with very low standard deviation ranging from 0.0007 to 0.12. The highest proportion was in the northern Sichuan Basin and the northern parts of the Yunnan-Guizhou Plateau, where the diffuse radiation fraction was highest. Conversely, in most parts of northwestern and northeastern China, the diffuse radiation fraction was less than 0.5 (Figure 2e), where the proportion of GPP contributed by Rd to GPP derived from Rg was greater than 0.5, reflecting the higher LUE of Rd.

Simulated Interannual Variability of Total GPP in China
The total amounts of GPP for China during the period 1981-2012 from three simulations are shown in Figure 8. In S1, the mean annual total GPP for the period 1981-2012 ranged from 7.43 to 8.47 Pg C y −1 . The total amount of GPP was lowest in 2000, likely because of the large-scale severe drought in China at that time. In S2, the multi-year mean annual GPP for the period 1981-2012 was 7.73 Pg C y −1 , which was 0.09 Pg C y −1 higher than the average value for S1. The difference in GPP between S1 and S2 reflects the effects of the 1.3% reduction in Rg, which caused a reduction of 1.7% in GPP. In S3, the multi-year mean national total GPP was 7.58 Pg C y −1 for the period 1981-2012, which was higher than the average value for S1. The difference in GPP between S1 and S3 reflects the

Simulated Interannual Variability of Total GPP in China
The total amounts of GPP for China during the period 1981-2012 from three simulations are shown in Figure 8. In S1, the mean annual total GPP for the period 1981-2012 ranged from 7.43 to 8.47 Pg C y −1 . The total amount of GPP was lowest in 2000, likely because of the large-scale severe drought in China at that time. In S2, the multi-year mean annual GPP for the period 1981-2012 was 7.73 Pg C y −1 , which was 0.09 Pg C y −1 higher than the average value for S1. The difference in GPP between S1 and S2 reflects the effects of the 1.3% reduction in R g , which caused a reduction of 1.7% in GPP. In S3, the multi-year mean national total GPP was 7.58 Pg C y −1 for the period 1981-2012, which was higher than the average value for S1. The difference in GPP between S1 and S3 reflects the increase in the diffuse radiation fraction over the period 1981-2012, which positively affected GPP in China and partially compensated for the loss in GPP because of the reduction in R g .
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 21 increase in the diffuse radiation fraction over the period 1981-2012, which positively affected GPP in China and partially compensated for the loss in GPP because of the reduction in Rg.

Spatial Distribution of Radiation Effects on GPP
The average GPP in the period 1981-2012 simulated in S2 ( Figure 9a) and S3 (Figure 9c) were similar to the distribution for S1 ( Figure 6). The difference (Figure 9b) between the multi-year mean GPP simulated in S1 and S2 showed that any reduction in Rg would result in a reduction in GPP over most of China. In North China, the decrease in Rg was significant, simultaneously reducing the GPP by more than 50 g C m −2 y −1 . The Rg over Central and Southwest China showed no significant decrease, but these areas experienced a decrease in GPP that was higher than 50 g C m −2 y −1 , which was predominantly because of the high LAI in these regions (see Supplementary Materials, Figure S5). Compared with areas with lower LAI, the reduction in Rg substantially reduced GPP for regions with higher LAI. During the period 1981-2012, the areas that showed the highest GPP increase were along the southeastern coast, northern Sichuan, as well as northeastern, northern, and some western parts of Xinjiang.
The difference in GPP values between S1 and S3 are shown in Figure 9d. In most regions of China, the enhanced diffuse radiation fraction led to an increase in GPP, especially in Central China and Yunnan Province. The significant increase in GPP in Central China mainly resulted from the increase in Rd. In Yunnan Province, although the increase in diffuse radiation fraction was not significant, the GPP increased because of the high LAI (see Supplementary Materials, Figure S5). On the Qinghai-Tibet Plateau, although the diffuse radiation fraction significantly increased, the increase in GPP was not as strongly significant as that in Central China because of the scarcity of vegetation in the former region. Similarly, along the southeastern coast, in Xi'an and western parts of Xinjiang, the reduction in diffuse radiation fraction led to a reduction in GPP.

Spatial Distribution of Radiation Effects on GPP
The average GPP in the period 1981-2012 simulated in S2 (Figure 9a) and S3 (Figure 9c) were similar to the distribution for S1 ( Figure 6). The difference (Figure 9b) between the multi-year mean GPP simulated in S1 and S2 showed that any reduction in R g would result in a reduction in GPP over most of China. In North China, the decrease in R g was significant, simultaneously reducing the GPP by more than 50 g C m −2 y −1 . The R g over Central and Southwest China showed no significant decrease, but these areas experienced a decrease in GPP that was higher than 50 g C m −2 y −1 , which was predominantly because of the high LAI in these regions (see Supplementary Materials, Figure S5). Compared with areas with lower LAI, the reduction in R g substantially reduced GPP for regions with higher LAI. During the period 1981-2012, the areas that showed the highest GPP increase were along the southeastern coast, northern Sichuan, as well as northeastern, northern, and some western parts of Xinjiang. Figure 9. Spatial distribution of the mean annual GPP simulated by scenario 2 (a); spatial distribution of the difference between the mean annual GPP in simulations 1 and 2 (b); spatial distribution of the annual mean GPP in simulation 3 (c); and spatial distribution of the difference between the mean annual GPP in simulations 1 and 3 (d) (unit: g C m −2 y −1 ). Simulation 1 assumed that the global radiation and diffuse radiation ratio changed with time; simulation 2 assumed that global radiation remained unchanged from 1981 to 2012; and simulation 3 assumed that the diffuse radiation fraction remained unchanged from 1981 to 2012 with a spatial resolution of 0.0727°.

Discussion
In China, annual mean Rd has increased during the past 30 years [25]. An increase in Rd leads to an increase in canopy LUE because photosynthetic saturation is reduced, thereby enhancing the ecosystem GPP [3,5]. Therefore, investigating the effect of Rd on GPP is important in studies of the carbon cycle [5,8,9]. There are many causes of radiation changes, such as changes in the solar constant, clouds, aerosols, and water vapor [55]. As the world's largest developing country, rapid economic development and rapid population growth in China have led to the emission of increasing quantities of man-made aerosols into the atmosphere. Aerosols can simultaneously change the amounts of direct and diffuse radiation reaching the surface [56]. Li et al. [57] reported that radiation decreased by about 24.1 W m −2 from September 2004 to September 2005 because of the higher aerosol concentration by using aerosol observation data from the Xianghe Station in China. In addition, volcanic eruptions can also strongly affect global and diffuse radiation [58,59]. Ren et al. [25] observed that the eruption of Mount Pinatubo in the Philippines in 1991 caused a rapid increase in the amount of diffuse radiation in China in the early 1990s.
The TL-LUE model, which provides a distinct means to investigate the effects of Rd changes on GPP, simulates GPP for sunlit (diffuse and direct radiation) and shaded leaves (diffuse radiation) separately. In the present study, during the period 1981-2012, the average GPP contributed by Rd in China was 5.69 Pg C y −1 , which accounted for 74.48% of the total annual GPP. This finding shows that the contribution of Rd to GPP in China was greater than that of direct radiation. The higher Figure 9. Spatial distribution of the mean annual GPP simulated by scenario 2 (a); spatial distribution of the difference between the mean annual GPP in simulations 1 and 2 (b); spatial distribution of the annual mean GPP in simulation 3 (c); and spatial distribution of the difference between the mean annual GPP in simulations 1 and 3 (d) (unit: g C m −2 y −1 ). Simulation 1 assumed that the global radiation and diffuse radiation ratio changed with time; simulation 2 assumed that global radiation remained unchanged from 1981 to 2012; and simulation 3 assumed that the diffuse radiation fraction remained unchanged from 1981 to 2012 with a spatial resolution of 0.0727 • .
The difference in GPP values between S1 and S3 are shown in Figure 9d. In most regions of China, the enhanced diffuse radiation fraction led to an increase in GPP, especially in Central China and Yunnan Province. The significant increase in GPP in Central China mainly resulted from the increase in R d . In Yunnan Province, although the increase in diffuse radiation fraction was not significant, the GPP increased because of the high LAI (see Supplementary Materials, Figure S5). On the Qinghai-Tibet Plateau, although the diffuse radiation fraction significantly increased, the increase in GPP was not as strongly significant as that in Central China because of the scarcity of vegetation in the former region. Similarly, along the southeastern coast, in Xi'an and western parts of Xinjiang, the reduction in diffuse radiation fraction led to a reduction in GPP.

Discussion
In China, annual mean R d has increased during the past 30 years [25]. An increase in R d leads to an increase in canopy LUE because photosynthetic saturation is reduced, thereby enhancing the ecosystem GPP [3,5]. Therefore, investigating the effect of R d on GPP is important in studies of the carbon cycle [5,8,9]. There are many causes of radiation changes, such as changes in the solar constant, clouds, aerosols, and water vapor [55]. As the world's largest developing country, rapid economic development and rapid population growth in China have led to the emission of increasing quantities of man-made aerosols into the atmosphere. Aerosols can simultaneously change the amounts of direct and diffuse radiation reaching the surface [56]. Li et al. [57] reported that radiation decreased by about 24.1 W m −2 from September 2004 to September 2005 because of the higher aerosol concentration by using aerosol observation data from the Xianghe Station in China. In addition, volcanic eruptions can also strongly affect global and diffuse radiation [58,59]. Ren et al. [25] observed that the eruption of Mount Pinatubo in the Philippines in 1991 caused a rapid increase in the amount of diffuse radiation in China in the early 1990s.
The TL-LUE model, which provides a distinct means to investigate the effects of R d changes on GPP, simulates GPP for sunlit (diffuse and direct radiation) and shaded leaves (diffuse radiation) separately. In the present study, during the period 1981-2012, the average GPP contributed by R d in China was 5.69 Pg C y −1 , which accounted for 74.48% of the total annual GPP. This finding shows that the contribution of R d to GPP in China was greater than that of direct radiation. The higher contribution of R d to GPP was largely because of the high fraction of R d and was also associated with the higher LUE of R d . The present results showed that the increase in R d ratio could compensate for the reduction in GPP because of the lower photosynthetically effective radiation. These results will help to quantify the impacts of R d on GPP and enhance knowledge of the possible causes of carbon budget variation.
Although the mean annual GPP and the GPP of each vegetation type are comparable to those reported by many previous studies (Table 4), some marked differences are apparent. The possible reasons are, first, that model structures are different, not only between process models and LUE models, but also among various LUE models. Second, the modeling area and period are set differently. For example, in the Carnegie-Ames-Stanford approach (CASA) ( Table 4), Piao et al. [50] estimated that the total amount of GPP for the period 1982-2009 only varied by 2.66-3.16 Pg C y −1 , whereas the average annual GPP reported by Zhu et al. [51] for 1993 was 6.24 Pg C y −1 , in which maximum LUE was retrieved based on observations at 690 sites. Third, the types of input data used are inconsistent, covering different time ranges and differing in spatial and temporal resolutions.
Several uncertainties in the present study are acknowledged. First, the R g of 1981 as a baseline was used to investigate the effect of R g and diffuse radiation fraction on GPP in terrestrial ecosystems. The advantages of this method were that the annual R g remained unchanged from 1981 to 2012 and the variation tendency of daily R g variation to R g values in 1981 was consistent. However, this method would result in a degree of uncertainty regarding daily R g . For example, for a year that included many cloudy days, the R g of a sunny day would be overestimated.
Second, the TL-LUE model has been validated worldwide [42][43][44], showing that the model performance for GPP simulations is satisfactory. In the present study, we obtained a high R 2 (0.82) value for comparisons between the model-simulated GPP and observed GPP values. However, an extremely limited number of flux sites in China were used to optimize LUE sun and LUE shaded [43], which would cause their underestimation for some ecosystems, especially for crops. Compared with C 3 plants, C 4 plants photosynthesize by the combination of vascular bundle sheath cells and mesophyll cells, which have a higher CO 2 assimilation rate, resulting in higher LUE [60]. Therefore, the inclusion of a few sites dominated by C 4 plants to parameterize the model would result in underestimation of LUE sun and LUE shaded used in the model and the subsequent simulation would underestimate GPP values at some sites in China ( Figure 5). Therefore, additional flux tower sites in China should be used to validate the model in the future.
In addition, underestimation of LAI values would result in underestimation of simulated GPP values. For example, the LAI values at the GRASS sites used in the present study were significantly lower than those of the MOD15A2 LAI product (Figure 10), which resulted in underestimation of GPP values at these sites. The GPP uncertainties caused by LAI uncertainties have been investigated and the sensitivity of GPP estimates derived from TL-LUE was lower than that of large-leaf LUE models at the site scale [43]. In the present study, the uncertainties of GPP simulated by the TL-LUE model caused by LAI uncertainties at the national scale were investigated ( Figure 10). Overestimation of LAI by 10% would result in GPP overestimation by about 1.2-10%, whereas in most areas with higher LAI the overestimation of GPP was no more than 8%; similarly, underestimation of LAI by 10% would result in GPP overestimation by a similar magnitude to that mentioned above but in the inverse direction ( Figure 11). duration, omitting many of the drivers of dimming and Rd, which would result in Rd and Rg uncertainties, especially in heterogeneous areas. Figure 10. Comparison of leaf area index (LAI) data used in this study and the MODIS LAI product (MOD15A2) for grassland sites, as part of the validation process for estimation of GPP (details for each site were described in Table 3).

Haibei Neimeng
Duolun Dangxion Figure 10. Comparison of leaf area index (LAI) data used in this study and the MODIS LAI product (MOD15A2) for grassland sites, as part of the validation process for estimation of GPP (details for each site were described in Table 3).
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 21 duration, omitting many of the drivers of dimming and Rd, which would result in Rd and Rg uncertainties, especially in heterogeneous areas. Figure 10. Comparison of leaf area index (LAI) data used in this study and the MODIS LAI product (MOD15A2) for grassland sites, as part of the validation process for estimation of GPP (details for each site were described in Table 3).

Haibei Neimeng
Duolun Dangxion Furthermore, radiation measurements in China are extremely limited, comprising 103 and 81 sites for global radiation and diffuse radiation, respectively; therefore, an interpolation method was used to infer the national-scale R g and R d . The Ångström model is commonly used for radiation interpolation [29,30]. Interpolated radiation data are widely used for ecosystem-level simulations of water and carbon dynamics [44,61]. However, because of the limited available measurements, national R g and R d were interpolated to 500 m based on daily extraterrestrial radiation and sunshine duration, omitting many of the drivers of dimming and R d , which would result in R d and R g uncertainties, especially in heterogeneous areas.

Conclusions
In this study, we modeled the effects on GPP of solar radiation variation caused by changes in diffuse radiation in China over the period 1981-2012 using a two-leaf light use efficiency model (TL-LUE), which was driven by national global radiation in conjunction with other meteorological data and remotely sensed data. The main findings are as follows: (i) The annual mean R g in China significantly decreased over the period 1981-2012 by 2.18 MJ m −2 y −1 , whereas the diffuse radiation increased. A significant reduction in R g and increases in the diffuse radiation fraction were observed in the Beijing-Tianjin-Hebei region and on the Qinghai-Tibet Plateau. The annual mean R d predominantly showed an upward trend over the study period, especially on the southeastern coast and the Qinghai-Tibet Plateau.
(ii) During the period 1981-2012, the average annual total GPP of China's terrestrial ecosystems was 7.64 Pg C y −1 , comprising a contribution from R d of 5.69 Pg C y −1 . The total GPP decreased by 0.09 Pg C y −1 , although this decrease was partially compensated for by about 50% because of the increase in R d . Clearly, the increase in LUE of canopy light associated with R d partly compensated for the reduction in total GPP caused by the reduced R g in China.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/20/3355/s1, Figure S1: Distribution of sunshine duration, diffuse and global radiation station in China. Figure S2: Observed and interpolated daily global radiation at 12 sites in 2014; Figure S3: Observed and interpolated daily R d at 12 sites in 2014; Figure S4: Comparisons of estimated R g and R d between this study and that of Ren et al. [30]; Figure S5: Spatial distribution of mean annual LAI in China from 1981 to 2012; Table S1: Proportions of sites distributed in different aerosol optical depth (AOD) levels. Table S2: Proportions of sites distributed in different AOD levels. Table S3: Scenario simulations of the impact of radiation change on GPP.

Conflicts of Interest:
The authors declare no conflict of interest.