{'fraction_non_alphanumeric': 0.03036109064112012, 'fraction_numerical': 0.02236960615737329, 'mean_word_length': 4.952432010917243, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Disentangling the Impacts of Anthropogenic Aerosols on Terrestrial Carbon Cycle During 1850–2014

Abstract Aerosols have a dimming and cooling effect and change hydrological regimes, thus affecting carbon fluxes, which are sensitive to climate. Aerosols also scatter sunlight, which increases the fraction of diffuse radiation, increasing photosynthesis. There remains no clear conclusion whether the impact of aerosols on land carbon fluxes is larger through diffuse radiation change than through changes in other climate variables. In this study, we quantified the overall physical impacts of anthropogenic aerosols on land C fluxes and explored the contribution from each factor using a set of factorial simulations driven by climate and aerosol data from the IPSL‐CM6A‐LR experiments during 1850–2014. A newly developed land surface model which distinguishes diffuse and direct radiation in canopy radiation transmission, ORCHIDEE_DF, was used. Specifically, a subgrid scheme was developed to distinguish the cloudy and clear sky conditions. We found that anthropogenic aerosol emissions since 1850 cumulatively enhanced the land C sink by 22.6 PgC. Seventy‐eight percent of this C sink enhancement is contributed by aerosol‐induced increase in the diffuse radiation fraction, much larger than the effect of the aerosol‐induced dimming. The cooling of anthropogenic aerosols has different impacts in different latitudes but overall increases the global land C sink. The dominant role of diffuse radiation changes found in this study implies that future aerosol emissions may have a much stronger impacts on the C cycle through changing radiation quality than through changing climate alone. Earth system models need to consider the diffuse radiation fertilization effect to better evaluate the impacts of climate change mitigation scenarios.

indirect impacts, for example, on the terrestrial C cycle (Mercado et al., 2009;Zhang et al., 2019), which in turn affects the climate system through climate-carbon feedbacks (Friedlingstein et al., 2006).
One of the main mechanisms by which aerosols affect the land C cycle is through changing the quantity and quality of short-wave radiation received by vegetation, with an impact on photosynthesis. On the one hand, the absorption and scattering of solar radiation due to atmospheric aerosols decrease the amount of incoming solar radiation at the top of the vegetation canopy. In situ observations have shown negative relationships between aerosol optical depth (AOD) and incoming solar radiation at the land surface (Kanniah et al., 2010). This decrease may cause a decrease in photosynthesis in ecosystems which are light-limited, such as tropical forests (Nemani et al., 2003). On the other hand, aerosols alter the angular distribution of incoming light. The scattering of radiation by aerosols significantly increases the fraction of diffuse light, which is able to penetrate deeper and distribute more evenly in the canopy (Roderick et al., 2001). Since sunlit leaves are often light-saturated while shaded leaves are not, an increase in diffuse radiation may increase the photosynthesis of the entire canopy. In situ observations at ecosystem level show that diffuse radiation increases the light use efficiency (LUE) of forests, grasslands, and croplands (Alton, 2008;Alton et al., 2007;Choudhury, 2001;Gu et al., 2002;Healey et al., 1998).
The compensation issue of gross primary production (GPP, the photosynthesis rate) gain due to diffuse radiation increase procured by aerosols and clouds and GPP loss due to total radiation decline has been investigated at the site level (Alton et al., 2007;Hollinger et al., 1994;Kanniah et al., 2010) in the past. There is, however, no consensus on the quantitative role of each process. One reason is that the field measurements were performed in various ecosystems under different environmental conditions, highlighting the importance of understanding this question at a larger scale. Another reason is that the precise effects of decreased light and increased diffuse light fraction depend on the amount of aerosols in the atmospheric column, their optical properties, and the relative occurrences of clear and cloudy sky conditions. At larger scales, the aerosol impacts are generally investigated using atmospheric and land surface models (LSM) (Jones et al., 2003;Mahowald et al., 2011;Zhang et al., 2019). These studies rely on simulations with and without anthropogenic aerosols and estimate the aerosol impacts as the difference between such pairs of simulations. However, due to the lack of representation of direct and diffuse radiation in canopy light transmission in LSMs, the aerosol impacts investigated by the above studies did not include the fertilization effect of diffuse radiation. Mercado et al. (2009) addressed this problem by modifying the canopy light transmission scheme in the Joint UK Land Environment Simulator (JULES) model to represent the differential absorption of diffuse and direct radiation. They quantified the aerosol impacts on the land C sink using this new scheme driven by observation-based climate fields. Their results strictly presented the impacts of aerosols through changing radiation, but did not consider the impacts of aerosol-caused cooling and precipitation changes, which were suggested to also have considerable impacts (Jones et al., 2003;Mahowald et al., 2011;Zhang et al., 2019). Recently, Xie et al. (2020) compared the impacts of aerosol-induced diffuse radiation fraction, temperature, and vapor pressure deficit changes on vegetation productions in China, and found a dominant impact from the changes in diffuse light fraction. However, their study omitted the impacts of aerosol-induced dimming and precipitation changes. Also, their study area is restricted to an industrialized region with strong aerosol emissions, which is not very representative of the entire land area. Therefore, the need still remains to investigate and quantify the full impact of anthropogenic aerosols on the land C budget at the global scale.
Currently, there is no standard way to deal with diffuse radiation and to assess its impacts. Existing LSMs either use empirical equations to partition total radiation into its direct and diffuse components (Zhang et al., 2020), or consider the weighted average of radiation in the clear and cloudy fractions of a gridbox (Mercado et al., 2009). In both approaches, gridboxes are considered to be homogeneous and to receive single values of total radiation and diffuse radiation fraction. However, the averaging of diffuse radiation fraction over a LSM gridbox (∼55 km or larger) and time step (half hourly or longer) can potentially bias land C fluxes (Zhang et al., 2021). Furthermore, LSMs with a representation of diffuse light were evaluated using site level observations (Mercado et al., 2009;Zhang et al., 2020), where cloudy and clear conditions do not coexist at the same time step. Therefore, simulating C fluxes in a heterogeneous gridbox is a complex task. We argue that it is necessary to at least distinguish clear-sky and cloudy conditions within a model gridbox.
ZHANG ET AL. 10.1029/2021EF002035 2 of 21 For the above reasons, we present here a new study to quantify the impacts of anthropogenic aerosols on land C fluxes during the historical period (1850-2014) using a LSM that simulates photosynthesis in the clear and cloudy gridbox fractions separately. A set of factorial simulations based on the fields of the IPSL-CM6A-LR climate model is designed to investigate the total impact of anthropogenic aerosols and to separate the impacts of each aerosol-mediated mechanism including changing (a) radiation quality, (b) radiation quantity, (c) temperature, and (d) precipitation. The study period is the historical period from 1850 to 2014, during which anthropogenic emission aerosols have increased strongly and exerted a significant effect on climate change. The objectives of this study are to (a) quantify the impacts of aerosol-induced changes of diffuse radiation and climate change on the land C budget during the historical period; (b) determine whether the fertilization by diffuse radiation globally compensates for the effect of reduction of incoming radiation due to anthropogenic aerosols; (c) understand how aerosol-induced changes in radiation, temperature, and precipitation individually affect land C fluxes and the historical carbon budget.

ORCHIDEE_DF Subgrid Version
In order to investigate the full impacts of aerosols including the diffuse radiation fertilization effect on the land C fluxes, we performed simulations with ORCHIDEE_DF, a new development of ORCHIDEE trunk (v5453) (Krinner et al., 2005;Zhang et al., 2020). ORCHIDEE_DF has a two-stream canopy light transmission model which distinguishes direct and diffuse radiation and considers the differential light absorption of sunlit and shaded leaves. ORCHIDEE_DF has been intensively evaluated over 159 flux sites and proved to be able to capture the observed diffuse radiation fertilization effect (Zhang et al., 2020). A detailed description of ORCHIDEE_DF can be found in Zhang et al. (2020).
As stated in Section 1, the homogenization treatment of clear and cloudy fractions in a LSM gridbox may introduce biases in the C flux simulation. To avoid the issue, we modified ORCHIDEE_DF to calculate photosynthesis of each plant functional type twice in each gridbox, for example, under clear sky and cloudy sky conditions, respectively. The clear sky and cloudy sky photosynthesis rates are then weighted by cloud cover and integrated to grid level GPP (Figure 1). To realize this subgrid simulation, the incoming shortwave radiation and diffuse radiation fraction of each clear and cloudy subgrid are needed. These subgrid variables are usually not available in climate models, while the total incoming shortwave radiation, cloud fraction, and aerosol optical depth per gridbox can be often found, especially in data sets from Earth system models (ESMs). Using these grid level variables, we developed a module to obtain the missing subgrid clear and cloudy sky radiation and diffuse radiation fraction in ORCHIDEE_DF ( Figure 1). How we partition radiation and calculate diffuse radiation fraction is detailed in Appendix A.

Forcing Data and Experimental Design
The climate forcing in this study is derived from two simulations performed with the IPSL-CM6A-LR , an ESM consisting of an atmospheric model LMDZ (v3554), an ocean model NEMO (v9455), a land surface model ORCHIDEE (trunk v5661), and a coupler (OASIS). As a model participating to the Coupled Model Intercomparison Project Phase 6 (CMIP6, Eyring et al., 2016), IPSL-CM6A-LR ran a preindustrial control simulation (piControl) driven by constant external forcing representative of the year 1850, after a long spin up. The historical climate fields are taken from the IPSL-CM6A-LR historical r1i1p1f1 simulation, which is one of the simulations performed for the 1850-2014 period driven by all natural and anthropogenic forcing and initialized from the January 1, 1910 state of the piControl simulation. From the same initial conditions of the piControl simulation, another simulation, hist-piAer (r1i1p1f1, hereafter referred to as histNoAA) was performed by prescribing all the forcing terms to be the same as in historical r1i1p1f1 but keeping the anthropogenic aerosols at pre-industrial level (1850).
The two simulations provide daily climate surface or near-surface variables including 2-m maximum (Tmax) and minimum (Tmin) air temperature, precipitation (Precip), incoming shortwave radiation at land surface (SWdown), incoming longwave radiation at land surface (LWdown), near surface specific humidity (Qair), near surface air pressure (Psurf), 10-m eastward wind component (Wind_E), and 10-m northward wind component (Wind_N). As the modeled climate may be biased from the observations, causing errors ZHANG ET AL.
10.1029/2021EF002035 3 of 21 in C fluxes in simulations, we introduced a bias correction procedure using the observation-based CRU-JRA v1.1 data set as the reference to ensure the climate used in this study has the same baseline as the observation-based data set for the preindustrial period (1850)(1851)(1852)(1853)(1854)(1855)(1856)(1857)(1858)(1859). This method has been used in Zhang et al. (2019) and is described in Appendix B. The aerosol information was obtained from the monthly tropospheric aerosol optical depth at 550 nm (od550aer) and monthly stratospheric aerosol optical depth at 550 nm (od550_STRAT) as used to drive the IPSL-CM6A-LR ESM simulations. Because clouds change on short (subdiurnal) timescales and have a strong impact on diffuse radiation, 3-h total cloud fraction (clt) data from the same IPSL-CM6A-LR simulations is used. Due to the lack of observations of preindustrial aerosol and cloud properties and amounts, we did not perform the same bias-correction on aerosol and cloud fraction variables. All these variables which were originally provided on a 144 × 143 global grid were first re-gridded into a 2° × 2° grid as model input, and further interpolated to 30-min time steps in ORCHIDEE_DF. The climate variables are interpolated using the default method in ORCHIDEE LSM to generate reasonable diurnal cycle, and the aerosol optical depth and cloud fraction are interpolated linearly.
For the LSM simulations, to get equilibrated initial C pools, a spin-up simulation using the 1850-1869 bias-corrected climate from the historical scenario ESM simulation (hereafter all climate refers to the bias-corrected climate if not specified) was first performed for 400 years using the spin-up analytic mode of ORCHIDEE (Lardy et al., 2011), which accelerates the reach of equilibrium state for slow soil carbon pools. After the spin-up, a set of LSM experiments separating the impacts of anthropogenic aerosols were set up. The detailed information about the experimental design can be found in Table 1. ZHANG ET AL.  The hist and the NoAA LSM simulations were respectively driven by the full set of climate variables from the two ESM simulations described above. The difference between the two simulations corresponds to the full impact of anthropogenic aerosols due to both their effects on climate and diffuse light fraction. As radiation, precipitation, and temperature are the three most important factors controlling land carbon fluxes (Nemani et al., 2003), we performed the NoAASW, the NoAAprecip, and the NoAATair LSM simulations, which use the historical ESM climate but replace the SWdown, the precipitation and the Tair fields, respectively, by their counterpart from the histNoAA scenario. These factorial simulations aim to investigate the relative contributions of the aerosol effects from each one of these climate variables. The difference between the hist and the NoAASW simulations represents the impacts of the aerosol-induced full change in solar radiation, including both quality and quantity. Beside these experiments, we also set up the NoAARqlt LSM simulation, which used the same forcing as the hist simulation but recalculated the diffuse light fraction of the clear area of each gridbox using the tropospheric AOD from the histNoAA simulation. This simulation is used to investigate the impacts of the fraction of diffuse radiation (Fdf) change caused by anthropogenic aerosols alone. In the end, a LSM simulation without stratospheric aerosols from volcanoes (NoVA) was performed to investigate the volcanic impacts alone. It should be noted that the difference between the hist and the NoVA simulations includes only the impacts of the Fdf change from volcanoes but not the full climate impact caused by volcanic aerosols because the other climate fields driving the two simulations are the same. All the LSM simulations used the annual land cover and land-use fields from the ESA-LUH2v2 data set.

Aerosol-Induced Climate Change
The variation of bias-corrected global land yearly mean climate from the IPSL-CM6A-LR ESM historical and histNoAA experiments is shown in Figure 2. As the emission of anthropogenic greenhouse gases and aerosols did not significantly increase before the early part of the 20th century, the climate remains relatively stable during this period and the two scenarios show no substantial difference. For both scenarios, the land mean temperature, precipitation, and SWdown generally vary around global averages of 13°C, 800 mm yr −1 , and 196 W m −2 , respectively. The average global cloud cover over land regions is around 51%-52%, and the land-averaged tropospheric AOD average stays around 0.1 before 1940. The global stratospheric AOD is smaller than 0.01 in volcano-free years and can increase by 10-fold during large volcanic eruptions, for    (Figure 2a), both the historical and histNoAA Tair show an increase after the 1950s. However, the Tair under the histNoAA scenario increases faster, resulting in a 0.6°C warmer climate during the last 2 decades compared with the historical scenario. In terms of SWdown, a dimming effect is found along with the cooling from anthropogenic aerosols ( Figure 2c). The historical SWdown decreases by about 3 W m −2 since the 1950s, while the histNoAA SWdown shows no significant decrease during the same period. Compared with Tair and SWdown, no substantial differences are found on precipitation and cloud fraction between the historical and the histNoAA scenarios.
The spatial distributions of the ESM-simulated aerosol-induced climate changes are shown in Figure 3. Considering that the historical and the histNoAA scenarios differ significantly only after 1950s regarding tropospheric AOD (Figure 2e), we focus on the period 1961-2010 to investigate the aerosol impact. The Tair from IPSL-CM6A-LR indicates a widespread cooling effect of anthropogenic aerosols over land regions ( Figure 3a) except western North America, where the mean Tair values are not significantly different between the two scenarios. In contrast to the ubiquitously lower Tair, changes in precipitation and an aridity index (precipitation-potential evapotranspiration [P-PET]) are more heterogeneous (Figures 3b and 3f). A slight to medium drying effect of aerosols is found in central Africa, India, South China, and eastern US, while in Argentina and regions around the Mediterranean, anthropogenic aerosols are found to have a wetting effect. Compared with the precursor version of the ESM (IPSL-CM5A-LR) climate (Zhang et al., 2019), the climate from IPSL-CM6A-LR does not show a strong drying effect due to aerosols in Amazonia, while it neither shows a significant wetting effect of aerosols in this region, as detected by the other two CMIP5 models used in Zhang et al. (2019).
In terms of SWdown changes ( Figure 3c), a dimming effect of anthropogenic aerosols is found in 78% continental regions. The magnitude of the dimming is generally 0-5 W m −2 except over industrialized regions: East Asia, Europe, and East US, where the dimming is over 5 W m −2 . A 0-3 W m −2 brightening is found in ZHANG ET AL.
10.1029/2021EF002035 6 of 21 Alaska, East Russian, Central Africa, and Australia, where aerosols in one region have resulted in dynamical changes leading to fewer clouds and more incoming solar radiation ( Figure 3d).
As expected, emissions of anthropogenic aerosols and their precursors cause an increase of tropospheric aerosol optical depth, with the largest increase found in the most polluted regions (East Asia, India, Europe, and eastern US) ( Figure 3e). It should be noted that the emissions of anthropogenic aerosols do not always increase tropospheric AOD. For instance, a lower tropospheric AOD is found in Central Africa under the historical scenario compared with the histNoAA scenario ( Figure 3e). This might be due to an altered aerosol distribution in response to aerosol-induced atmospheric circulation changes.

Response of Land C Fluxes to Aerosols
The evolution of GPP and total ecosystem respiration (TER, the sum of autotrophic and heterotrophic respiration) from different LSM experiments are shown in Figure 4. All the experiments show a GPP increase from 124-125 PgC yr −1 in 1850 to about 165-168 PgC yr −1 in 2014 ( Figure 4a). This GPP increase is probably, to a large extent, due to the fertilization effect of rising atmospheric CO 2 concentration (Piao et al., 2013). Despite this consistent trend, differences among simulations remain significant. Compared with the NoAA simulation, the hist simulation has lower GPP in most of the years before the 1960s, but higher GPP afterward ( Figure 4c), indicating that the response of GPP to changing anthropogenic aerosols has changed by the 1960s, when aerosol emissions started to grow quickly ( Figure 2). Before 1960, the GPP difference between the hist and the NoAA simulations is very similar to the impacts of Tair alone (hist-NoAATair) ZHANG ET AL.
10.1029/2021EF002035 7 of 21 ( Figure 4c), while changes in SWdown or precipitation alone (NoAASW and NoAAprecip compared with hist) have small impacts on global GPP. This result indicates that the aerosol impacts on GPP is dominated by temperature changes before 1950s. However, after the 1960s, the full aerosol impact (hist-NoAA) on GPP starts to deviate from the temperature-alone induced change (hist-NoAATair) and becomes positive. This change is related to the increasing impact of aerosol-induced SWdown changes, as the impact from aerosol-induced precipitation changes remains small but the GPP difference between the hist and the NoAASW simulations increases strongly from <1 to 3 PgC yr −1 during this period. In addition, the NoAARqlt simulation shows very similar GPP compared to the NoAASW simulation, implying that the GPP change is more attributable to the change in radiation quality (diffuse radiation fraction) rather than quantity (dimming). By comparing GPP from the NoVA and the hist simulations, we are able to better understand the impact of the Fdf change caused by volcanic aerosols alone. For the years with strong volcanic eruptions, GPP increases of 5-10 PgC yr −1 can be found in all simulations except in NoVA (Figure 4a). This difference is generally larger than the slowly varying and weaker impacts of anthropogenic aerosols, however, it lasts for a short time. The stratospheric AOD can increase by over 0.15 after large eruptions, such as Pinatubo in 1991 (Figure 2f). This is larger than the tropospheric AOD increase caused by anthropogenic aerosols (Figure 2e). For years without volcanic eruptions, the hist and the NoVA GPP are generally identical, as expected.
In terms of TER, the global averages show similar increasing trend among the factorial simulations (Figure 4b). From 1850 to 2014, the global TER increases from 119-120 to 150-153 PgC yr −1 . This increase should be mainly due to the increase in GPP, which provides more organic matter for respiration, and the increase in temperature, which accelerates the turnover of soil C pools. The difference in TER between the hist and the other experiments shows very similar patterns as the difference in GPP (Figures 4c and 4d). Before the 1950s, anthropogenic aerosols cause a 0-1 PgC yr −1 decrease in TER mainly through their impact on temperature (cooling, causing a reduced GPP). While after the 1960s, a positive impact of aerosol-caused SWdown changes (quality and quantity) on TER is found and the total TER response to aerosols finally becomes positive after the 2000s. The impacts of volcanic aerosols on TER are smaller but last much longer than GPP. After each eruption, following a sudden increase, a gradual decline of TER can be found during the following 3-4 decades or until the next volcano eruption (Figure 4d).
ZHANG ET AL.
10.1029/2021EF002035 8 of 21 The cumulative net biome production (NBP, the net C flux, calculated as GPP−TER−C release from natural and anthropogenic disturbances) changes caused by aerosols are displayed in Figure 5. As a result of the decreased GPP due to aerosols before 1950s, aerosol-induced NBP change (hist-NoAA) is negative (anthropogenic aerosols decrease land C sink) before the 1950s. After 1950s, the cumulative NBP difference becomes positive and grows quickly. At the end of the study period (2014), the hist simulation shows a 22.6 PgC higher cumulative NBP than the NoAA simulation (Figure 5b), indicating a large positive impact of anthropogenic aerosols on the cumulative land C storage during the historical period.
Besides the full impact of anthropogenic aerosols, we also investigated the impact of each environmental factor on cumulative NBP in Figure 5 by comparing the hist simulation and each factorial simulation. In comparison to hist, the cumulative NBP from the NoAATair simulation decreases before 1880s and stays around −6 PgC until 1950s (Figure 5a). An investigation on the historical and histNoAA climate shows that during the first two decades, the historical climate is much cooler than histNoAA in Eurasia, leading to smaller GPP in the hist simulation ( Figure S1). Given that the radiative forcing from anthropogenic aerosol emissions remains small at global scale during this period, the large difference in temperature between ESM simulations historical and histNoAA should be interpreted as ESM's internal variability (most likely stochastic) triggered by the change in implementation of aerosol emissions in the two ESM simulations. After the 1950s, the impact of aerosol-induced Tair changes reversed and the difference in cumulative NBP increased to a positive value of 3.6 PgC in 2014.
Compared with the impacts of Tair, the impacts of precipitation on cumulative NBP (hist-NoAAprecip) show much smaller impacts, within 0-3 PgC during the entire study period (Figure 5a), indicating limited impacts of aerosol-induced precipitation changes.
In terms of radiation, the cumulative NBP differences between the hist and the NoAARqlt is small before the 1950s, indicating limited impacts of aerosol-induced radiation quality changes during this period. This small impact is not surprising as the anthropogenic aerosol-induced change in light quality is small during this period (Figures 2d-2f). After 1950s, along with the increasing aerosol emissions, strong increases in cumulative NBP difference between the hist and the NoAARqlt simulations are found (Figure 5a). The cumulative NBP difference between the hist and the NoAASW is largely consistent with the difference between the hist and the NoAARqlt, but with slightly smaller magnitude. Until 2014, aerosol-induced total radiation changes (hist-NoAASW) and diffuse light fraction changes alone (hist-NoAARqlt) enhanced the cumulative global land C sink by 14.7 and 17.7 PgC, respectively. As the total radiation changes include changes in light quantity (total SWdown), cloud fraction and the diffuse light fraction, the aerosol-induced changes in light quantity and cloud fraction have inhibited the cumulative land C sink by 3.0 PgC. Considering the full impacts of anthropogenic aerosols of 22.6 PgC, the changes in diffuse light fraction alone contributed to about 78% of the increase in land C sink from anthropogenic aerosols. If we only exclude the first 30 years when the climate internal variability shows a large impact on C fluxes and consider the period ZHANG ET AL.
Unlike anthropogenic aerosols, volcanic aerosols are emitted sporadically during years with large eruptions. In response, the immediate (1-2 years after eruption) land C sink can surge by over 6 PgC depending on the intensity of the eruption (Figure 5a). These anomalies do not disappear right after the periods during which the stratosphere is influenced by the volcanic eruptions. Instead, the cumulative NBP difference between the hist and the NoVA simulation retains a memory and relaxes slowly toward zero during the following decades, at least until the next eruption takes place. In 2014, the diffuse radiation changes caused by volcanic aerosols cumulatively increased the cumulative land C sink by 3.9 PgC, which is less than the persistent changes in diffuse light from anthropogenic aerosols.

Regional C Flux Changes
To further understand how different factors affected the aerosol impacts regionally, we further compared the NBP difference between the hist and the factorial simulations in different latitudes ( Figure 6).
In the northern high latitudes (>60°N), all factorial simulations except the NoAATair show small cumulative NBP difference (<0.6 PgC) compared with the hist simulation (Figure 6a), indicating small impacts of aerosol-induced radiation and precipitation changes in this region. In contrast, the temperature-alone changes from anthropogenic aerosols cause a remarkable decline in cumulative NBP (hist-NoAATair) over those temperature-limited ecosystems during the entire study period. This decline is small before the 1950s but becomes much more pronounced afterward, resulting in a 2.3 PgC lower cumulative C sink in the hist than in the NoAATair simulation by 2014. The full impact of anthropogenic aerosols (hist-NoAA) shows very similar trend and interannual variation compared with the impact of aerosol induced temperature changes only (hist-NoAATair), implying that anthropogenic aerosols are mainly causing a decline in C sink through their cooling effect in this region.
In contrast to the northern high latitudes, the northern mid latitudes (30°-60°N)  (hist-NoAA) is generally in line with impacts exerted by SWdown (hist-NoAASW) or by only diffuse radiation fraction (hist-NoAARqlt), slightly offset by the negative impact of aerosol-induced Tair changes (hist-NoAATair), indicating that the impact of aerosols in this region is mainly through the changes in light quality offsetting the effect of cooling. The difference between the hist and the NoAAprecip simulations remains small in this region. Until the end of the study period in 2014, anthropogenic aerosols cumulatively increased NBP in this region by 8.5 PgC. SWdown, Fdf, and Tair, respectively increased NBP by 9.3, 8.9 PgC and decreased NBP by 2.4 PgC.
Unlike the northern mid and high latitudes, where aerosol induced NBP changes are dominated by a single environmental factor, no single factor is found to explain the full impacts of anthropogenic aerosols at low latitudes (30°S-30°N, Figure 6c). In this region, the cumulative NBP differences between the hist and each of the NoAA, the NoAATair, and the NoAASW simulations decrease in the first half of the study period and increase rapidly after the 1950s. Until 2014, the overall effect of anthropogenic aerosols increases the cumulative C sink in this region by 14.9 PgC, with changes in Tair, Fdf, and SWdown explaining 8.5, 8.1, and 3.9 PgC, respectively. Note that in the tropics, cooling induced by aerosols has a positive effect on NBP as it reduced the evaporative demand and allows more photosynthesis, unlike in the mid and high latitudes where biomes are mainly temperature limited. Compared with other variables, the impacts of aerosol-induced precipitation changes are small for the low latitudes.
The impacts of aerosols on the land C sink in the southern mid latitudes (30°S-60°S) is smaller than in other latitudes because the land area in this region is relatively small (Figure 6d). Until the end of the study period, the anthropogenic aerosols enhanced the C sink by 1.0 PgC (hist-NoAA), which is close to the single impact of aerosol-induced precipitation changes (hist-NoAAprecip, 1.0 PgC) and slightly smaller than the impact of SWdown changes (hist-NoAASW, 1.5 PgC). The changes in Tair cause a small decline of cumulative NBP (0.5 PgC) before 1900 and an increase of C sink after 1960s. At the end of the study period, the difference between hist and NoAATair is only −0.3 PgC in cumulative NBP. It should be noted that unlike the other latitudes where the impact of radiation quality changes is larger than the full impacts of radiation changes, in this region, the cumulative NBP difference between hist and NoAARqlt is much smaller than that between hist and NoAASW. This is because the anthropogenic aerosol emissions are small in this region and have little impact on tropospheric AOD and thus diffuse radiation fraction (Figure 3e). In addition, the SWdown increased in a considerable proportion of the land area in this region in response to anthropogenic aerosol-induced cloud fraction changes (Figures 3c and 3d), this brightening may be responsible to the large difference in radiation impacts between this region and the other latitudes.

Spatial Distribution of the Aerosol Impacts
Since the aerosol-induced flux changes mainly occur after 1950s, to further understand how this happens, we examine the detailed regional spatial patterns of GPP and NBP differences between the hist and each factorial simulation during 1961-2010 in Figures 7 and 8.
The full impacts of anthropogenic aerosols on GPP generally show a latitudinal pattern (Figure 7a), with most regions in low latitudes experiencing a positive impact of anthropogenic aerosols and most mid to high latitude regions a negative impact, except for East Asia, eastern US, and Europe, where the most pronounced tropospheric AOD increases in response to anthropogenic aerosol emissions are found (Figure 3). The strongest GPP decrease in response to aerosols is found in central and western US and northern Russia, with a magnitude of over 50 gC m −2 yr −1 .
The impacts of aerosol-induced Tair changes show a similar latitudinal pattern as the full aerosol impacts (Figure 7b), that is, a decreased GPP at mid to high latitudes and an increased GPP at low latitudes. In contrast to the full impacts of anthropogenic aerosols, the Tair changes alone from aerosols generally decreased GPP in East Asia, eastern US, and Europe. In contrast to Tair, aerosol-caused precipitation changes have relatively small impacts on GPP in most regions, except for semi-arid regions (Figure 7c). In over 85% of the land area, the precipitation changes cause less than 10 gC m −2 yr −1 changes in GPP (generally <1% of annual GPP). Nevertheless, in relatively dry regions like the Sahel and Central India, GPP decreases significantly in response to aerosol-induced precipitation changes (Figure 7c), which partly explains the total aerosol-induced GPP decline in these regions (Figure 7a). Similar positive responses of GPP to aerosol-induced ZHANG ET AL.

10.1029/2021EF002035
11 of 21 precipitation changes are found in the semi-arid regions in southern South America and southern Africa. The impact of radiation (Figures 7d and 7e) changes translates into a spatial pattern strongly related to the anthropogenic aerosol emissions (Figure 3). Both NoAASW and NoAARqlt have over 50 gC m −2 yr −1 higher GPP than hist in the most polluted regions, that is, East Asia, eastern US, and Europe. Meanwhile, the dimming due to aerosols tends to decrease the GPP in these polluted regions (Figure 7f). Generally, the impact of aerosol-caused overall radiation changes on GPP are more dominated by changes in radiation quality (Figures 7d and 7e), except regions in India, Southeast Asia, and Amazon, where GPP is decreased due to the changes in light quantity (mainly dimming) and cloud fraction (Figures 7d and 7f). The simulated light limitation of GPP at low latitudes is in line with analysis of climatic data and satellite observations of vegetation activity by Nemani et al. (2003).
The NBP response to anthropogenic aerosols generally shows a similar spatial pattern as that of GPP, with a positive response in low latitudes and a negative response in high latitudes (Figure 8a). Compared to GPP, the area of NBP decrease in mid and high latitudes is much smaller. Large regions in North America and Russia show negative impacts of anthropogenic aerosols on GPP but positive impacts on NBP (Figures 7a  and 8a). These patterns are in line with the Tair change-induced C flux changes (Figures 7b and 8b), that is, at the mid latitudes, aerosol-induced cooling enhances the regional land C sink by decreasing more respiration than photosynthesis, which is consistent with a previous study using simulations driven by CMIP5 climate data (Zhang et al., 2019). In contrast to temperature, which has strong impacts on heterotrophic respiration, the aerosol-induced radiation changes do not directly alter respiration, but only indirectly through litter input. Not surprisingly, the spatial pattern of radiation change-induced GPP and NBP changes are ZHANG ET AL.

10.1029/2021EF002035
12 of 21 similar (Figures 7d-7f and 8d-f). NBP is found significantly increased in response to the radiation changes in forested regions with high aerosol emissions such as East and South Asia, eastern US, and Europe. The responses of NBP to aerosol-induced precipitation changes are small on the 50-years time scale (Figure 8c).

Aerosol Impacts on GPP for Clear and Cloudy Days
Because under cloudy condition the Fdf is always 1, the impacts of aerosols on C fluxes strongly depends on whether the sky is clear or cloudy. Figure 9 clearly illustrates this difference. In clear days (Figure 9a  anthropogenic aerosols are found to enhance the GPP of most land regions, especially of the most industrialized regions including Europe, India, East Asia, and East US. The maximum GPP difference can be over 200 gC m −2 yr −1 , which is found in East Europe. Some exceptional negative impact on GPP is found in high latitude where ecosystems are limited by temperature. In contrast to the clear days, the impact of aerosols on GPP in cloudy days show a meridional pattern with negative impacts in mid and high latitudes and positive impacts in tropical regions (Figure 9b), which is generally consistent with the Tair-alone caused GPP changes (Figure 7b). Some exceptions are found in India with negative GPP impact of over −150 gC m −2 yr −1 . This negative impact might be attributed to strong decline of SWdown during cloudy days in this region ( Figure S2), possibly due to aerosol-induced circulation changes. Because NBP has a lagged response to GPP, the aerosol-induced NBP changes under cloudy and clear days are much noisier than the GPP changes ( Figure S3).

Global Impacts of Aerosols on C Fluxes
Anthropogenic aerosols are found to decrease the simulated global GPP by 0.22 PgC yr −1 during the entire study period . Along with the decline in GPP, anthropogenic aerosols also induced a decrease of global TER of 0.37 PgC yr −1 , which is stronger than the GPP decrease, leading to a global increase of land C sink by 0.14 PgC yr −1 , or cumulatively 22.6 PgC since 1850. This increase of land C sink is about one half of the cumulative NBP (51 PgC) during this period and can compensate about 2 years of anthropogenic CO 2 emissions during the last decade (Friedlingstein et al., 2014). It should be noted that there is a cumulative NBP difference of about 6 PgC before 1880 due to a particular cold period of 2-3 decades at the beginning of the ESM hist simulation, attributable to ESM's internal variability (Figures 5 and S1). If we account for the anthropogenic aerosol impacts since 1880 and masked out the period before, the cumulative aerosol impacts on NBP would be larger, with the magnitude of 26.5 PgC ( Figure S4).
When considering only the period after the 1950s, we found an increase of GPP in response to anthropogenic aerosols 0.27 PgC yr −1 , which differs significantly from the GPP decline found before the 1950s. This reversed pattern is very different from a recent study using a comparable approach but omitting the changes in light quality omitted (Zhang et al., 2019). They showed that aerosol impacts on GPP did not change much before and after the 1950s. In our model, the changes in diffuse radiation after the 1950s have significantly altered the aerosol impacts.
Among the large overall impact of anthropogenic aerosols, the changes in diffuse radiation fraction accounted for a 17.7 PgC (17.6 PgC after 1880) larger cumulative land storage. This large enhancement of land C storage due to light quality change alone indicates that the diffuse radiation fertilization effect reported at flux sites for short term periods (Alton et al., 2007;Cheng et al., 2015;Gu et al., 2002;Niyogi et al., 2004;Misson et al., 2005;Kanniah et al., 2013;Knohl & Baldocchi, 2008;Oliphant et al., 2011;Williams et al., 2014) also have a long-term historical influence on the terrestrial carbon budget. Mercado et al. (2009) used JU-LES to investigate the impacts of Fdf changes on land C sink and reported an increase in C sink of 0.37 PgC yr −1 during 1960-1999. Compared with their estimation, our simulations show very similar impact from aerosol-caused light quality changes (0.373 PgC yr −1 ).

The Impacts of Radiation Quality Versus Radiation Quantity
One of the key questions on aerosol impacts is whether the effect of increased diffuse radiation offsets the effect of decreasing total radiation (dimming). By investigating C fluxes at different AOD levels, Niyogi et al. (2004) found that higher AOD increases C sink at forest and cropland sites but decreases C sink at grassland sites. This could be explained by the different canopy architecture of different ecosystems. The physical mechanism underlying the diffuse radiation fertilization effect relates to shadows. With diffuse light, leaves that are packed together in a vertical sense can still receive sufficient light to maintain a high photosynthesis rate. Alternatively, under sunny conditions, leaves cast shadows that go a long way below each leaf. In short vegetation with small LAI (e.g., grassland), the shadows are more casted on soil than leaves, leading to a weaker diffuse radiation fertilization effect. In contrast, in tall vegetation with larger LAI (e.g., forests), the shadows under sunny conditions are more casted on leaves in the lower part of the ZHANG ET AL.

10.1029/2021EF002035
14 of 21 canopy, causing a deficit in radiation that cannot maintain a high photosynthesis rate on these leaves. In an extreme case, for a canopy with only one layer of leaves, all leaves should be sunlit and the increase in diffuse radiation would give no benefits to GPP. In this study, we found that substantial C sink increases are only found in regions with mean yearly LAI over 2 ( Figure S5).
Apart from LAI, the overall radiation impact (quality and quantity) is suggested to also depend on whether the diffuse radiation is from the scattering of aerosols or clouds, as the latter media is thought to be less conservative regarding total downward shortwave radiation (i.e., have stronger dimming effect, Alton et al., 2008). In this study, we used a sub-grid scheme to separate the impact from cloud and aerosols. And our results shows that the GPP increase due to aerosols mainly happens in days with fewer clouds (Figure 9), which is consistent with another model-based study (Yue & Unger, 2017) Using the NoAASW and the NoAARqlt simulations we were able to separate the effect of changes in light quantity and light quality. Globally the overall radiation impact is an increase of carbon storage of 14.7 PgC, this large enhancement indicates that the aerosol-induced radiation quality change has much stronger impact than radiation quantity changes. If we omit the difference in cloud cover, which is small at global scale ( Figure 2), and consider the difference between the NoAASW and the NoAARqlt simulations as the impact of radiation quantity due to aerosols, we estimate the dimming-induced cumulative land C sink decrease to be 3.0 PgC during the study period. Regionally, our results also show that the radiation quality dominates the radiation change-caused C flux changes in most regions. Some exceptions are found in a few tropical regions such as Southeast Asia, where the amounts of anthropogenic aerosols are relatively small (Figures 7  and 8) and the forests are more light limited (Nemani et al., 2003), therefore the dimming-caused GPP and C sink decline dominates.

The Main Factors Causing C Flux Changes
Only few studies have attempted to quantify the aerosol impacts on the terrestrial carbon cycle at global scale, and fewer disentangled the relative contributions of different environmental factors. By analyzing model output with a statistical method, Zhang et al. (2019) found that the increase of the land C sink in response to anthropogenic aerosols is mainly due to their cooling effect, which enhances tropical GPP and causes a stronger decline in TER than in GPP in mid latitudes. In this study, using factorial simulations with a process-based land surface model evaluated against many flux tower measurements (Zhang et al., 2020), we found similar latitudinal pattern of cooling-caused C flux changes (Figures 7 and 8) although the increased global land C sink due to this effect is much smaller than in Zhang et al. (2019). This difference in magnitude is mainly due to the internal variability of the climate system before 1880. When considering the period after 1880, the temperature impacts (10.0 PgC) detected in this study is of similar magnitude as in Zhang et al. (2019) (11.6-41.8 PgC) ( Figure S4).
In contrast to temperature, Zhang et al. (2019) suggested that the aerosol-induced dimming has limited direct impacts on land C fluxes as most land ecosystems are temperature or precipitation limited rather than radiation limited (Nemani et al., 2003;Piao et al., 2009). However, this study detected a larger contribution of aerosol-induced radiation changes than temperature to global NBP changes (Figure 5b). This does not contradict the study by Zhang et al. (2019) because the large contribution is attributed to changes in light quality but not quantity (Figures 7 and 8). The light quantity changes induced by aerosols, which are generally equivalent to the radiation changes investigated by Zhang et al. (2019) in a LSM with one-stream canopy light transmission, are found to slightly decrease the GPP and NBP in intensely polluted regions (Figures 7 and 8). However, this dimming-induced decrease is not comparable in magnitude neither to the impacts of light quality changes, nor to those of temperature changes.
Compared with temperature and radiation, the impacts from changing precipitation are generally small in this study. However, this small impact from precipitation changes should be treated with caution because precipitation in current ESMs often has large uncertainties, especially in the tropics (Mehran et al., 2014). As shown by Zhang et al. (2019), if anthropogenic aerosols have a strong drying effect in very productive tropical regions, the regional C budget can be strongly affected. In IPSL-CM6A-LR climate, aerosol-induced drying is mainly found in the Sahel region and in South Asia (Figure 3b). In response, the GPP in these ZHANG ET AL.

10.1029/2021EF002035
15 of 21 regions is slightly smaller in the hist than in the NoAAprecip simulation ( Figure 7c). This regional change has limited impact on the global C fluxes.

The Impacts of Volcanoes
Volcanoes are a very important sporadic source of stratospheric aerosols that affect light globally and they have been reported to significantly enhance GPP at the site level due to their impact on Fdf (Gu et al., 2002). However, at the global scale, there remain no available reliable observational datasets to investigate the impacts of volcanic aerosols. The eruption of Mount Pinatubo is covered by several atmospheric inversion C flux data sets (e.g., Rödenbeck et al., 2003), but it is difficult to separate the impacts of aerosols from those of climate variations, especially given the fact that an El Nino event occurred right after the eruption. Based on models, the eruption of Mount Pinatubo is estimated to have induced an increase in land C sink of 1-2.5 PgC yr −1 in the 1-2 years following the eruption, which is suggested to be strongly contributed by the Fdf changes caused by the volcanic aerosols (Mercado et al., 2009;Roderick et al., 2001). In this study, our simulations show a sudden increase of global land NBP of 3.2 PgC yr −1 in response to the aerosols emitted by Mount Pinatubo and about 1.2-1.3 PgC yr −1 in response to the eruptions of Santa María and El Chichón. These volcanic aerosol impacts partly explain the observed drop of atmospheric CO 2 growth rate after the eruption and are generally in line with the previous estimates (Mercado et al., 2009;Roderick et al., 2001). It should be noted that we only considered the change in diffuse radiation fraction in response to volcanic aerosols, but did not include the cooling and dimming caused by these eruptions. To attribute the full impacts of volcanoes, additional simulations with IPSL-CM6A-LR would be needed to provide the climate fields.
Because volcanic aerosols have a strong cooling effect on the global land surface through reducing the solar radiation, stratospheric aerosol injection technique, which mimics the eruption of volcanoes, has been proposed as a possible solar radiation management technique to mitigate global warming (Boucher et al., 2013). Some strategies consist in adjusting aerosol injections to keep the stratospheric aerosols at a high level (i.e., Tjiputra et al., 2016), in order to keep a low surface temperature. In this study, we found that the impact of such a sudden change in stratospheric aerosols has persistent impacts on land C budget after 3-5 decades ( Figures 5 and 6), implying that the stratospheric aerosol injection technique may have long-lasting atmospheric carbon dioxide removal consequences, but with regionally diverse impacts. This effect needs to be considered when evaluating current solar radiation management strategies.

Limitation, Implication and Future Work
In this study, we set up factorial simulations to investigate the aerosol impact during the historical period using IPSL-CM6A-LR. Although we tried to reduce the uncertainty in the ESM climate through the bias correction, there remain unavoidable uncertainties.
First, the cloud fraction from the current ESMs is still not very accurate, compared with satellite observation. The error in cloud fraction is not likely to alter the impact of diffuse radiation in the clear and cloudy parts of the gridbox alone because the radiation fluxes in the clear part of the gridbox is extracted from the look up table (Appendix A) and the diffuse radiation fertilization effect under cloudy conditions is small ( Figure 9). Nevertheless, the error in cloud fraction may still affect the total impact of aerosols because it is used as a weighting factor to combine the clear and cloudy C fluxes at the subgrid scale. A recent study evaluated the cloud fraction from CMIP6 models and found that current models can generally capture the variation of cloud cover but with a SE of 0.5 compared with recent satellite observations (Vignesh et al., 2020), which is not that huge. In this study, we did not correct the bias in cloud fraction like the other climate variables because there remain no reliable cloud cover reconstructions at the preindustrial level. Even though it is possible to do the correction using present day observations, we still cannot ensure to have correct cloud cover for the early decades of our simulations and under the histNoAA scenario due to the lack of understanding on the evolution of cloud cover. To reduce this uncertainty, ESMs which can better simulate cloud cover are needed.
Second, the look up table for sub-grid radiation partitioning uses a different atmospheric radiation transfer model than the one used in IPSL-CM6A-LR, which may result in a mismatch between the total and diffuse radiation fluxes at the subgrid level. Besides this, the IPSL-CM6A-LR ESM uses a three-dimensional cloud ZHANG ET AL.

10.1029/2021EF002035
16 of 21 distribution and complex overlap assumptions to simulate radiation, the use of two-dimensional cloud cover to partition radiation to clear and cloudy subgrid in this study (Appendix A) will also cause mismatches in the subgrid level radiations. Furthermore, due to the lack of aerosol type information in IPSL-CM6A-LR to match the aerosol types in the look-up table, we have to assume a single type of aerosol in each region, which introduces some uncertainty in the Fdf.
Third, as shown in our results ( Figures 5 and S1), the internal variability of the climate system can also affect the C fluxes. To eliminate this impact, ensemble simulations with IPSL-CM6A-LR starting from slightly different initial state are needed to provide the histNoAA climate with different internal climate variabilities. With limited resources, we are not able to do such simulations now. Nevertheless, the impacts from climate internal variability are not likely to change our main conclusions because the magnitude of the cumulative NBP change due to climate internal variability (about 6 PgC) is much smaller than the aerosol impact we detected after the 1950s (over 20 PgC). In addition, the GPP and cumulative NBP from the IPSL-CM6A-LR historical ensemble simulations (with different climate internal variability) shows high consistency, and have standard deviations less than 2.5 PgC in most of the years ( Figure S6), which is much less than the aerosol impacts detected in this study.
In spite of the bias correction of climate variables, there remain uncertainties in the aerosol-induced climate changes (i.e., difference between historical and histNoAA). As shown by Zhang et al. (2019), the ESM uncertainty in aerosol-induced climate changes, especially precipitation changes, can potentially cause large uncertainty in the C fluxes. However, the histNoAA is not a mandatory experiment in CMIP6, and few models have performed it. Also, not all CMIP6 data have been published currently. It thus remains difficult to quantify the uncertainty in this effect from a multimodel perspective at this moment. One may argue that many CMIP6 ESMs have the hist-aer experiments (i.e., only consider the changes of tropospheric aerosols but keep the other anthropogenic forcings constant) and pairing these experiments with their corresponding CMIP6 piControl experiments (no anthropogenic forcings considered) may also be used to quantify the aerosol impacts. However, the piControl experiments have different climate baseline compared with the historical experiments, thus the aerosol impacts estimated from hist-aer and piControl climate must be biased from the real aerosol impacts.
Besides the forcing field, uncertainties may also arise from the land surface model used. In this study we used only ORCHIDEE_DF to investigate the Fdf and climate effects. Although the capacity of ORCHIDEE_ DF to capture the diffuse light fertilization effect has been carefully evaluated at a large number of flux towers (Zhang et al., 2020), there remain structural uncertainties due to the other processes in ORCHIDEE_DF. To understand this uncertainty, a multi-model study on aerosol impacts would be needed to confirm the considerable impact from aerosol-induced radiation quality changes found in this study. Candidate LSMs for such a multimodel study must have a two-stream canopy radiation transmission to distinguish direct and diffuse light, as it is the case of the JULES (Mercado et al., 2009), YIBs (Yue & Unger, 2015), iTem LSM (Chen & Zhuang, 2014), and CLM (Oleson et al., 2013) models. Besides the uncertainties in current model processes, the lack of other physiological processes may also bias the results. For instance, the ORCHIDEE_ DF model includes no nutrient limitation on vegetation growth, which means the deposition of nutrients along with the deposition of aerosols is not considered in our simulations. It has been suggested that N deposition can be a main factor causing land C sink in some ecosystems (e.g., Magnani et al., 2007;Vuichard et al., 2019). Understanding how this effect interacts with aerosol-induced climate and Fdf requires land surface models with additional complexity.
It should be also noted that the aerosol impacts investigated in this study lack a range of potentially important feedbacks. The increased land C sink will slow the growth of the atmospheric CO 2 concentration, which can decrease the CO 2 fertilization effect on vegetation and decrease land C sink. Meanwhile, the slowed growth rate of atmospheric CO 2 will slow the warming and increase land C sink. Besides these mechanisms, the changes of land C cycle may also alter the land surface albedo and roughness by affecting the vegetation, which will cause biophysical feedbacks to climate. To investigate these feedbacks, coupled simulations with ESMs that consider the diffuse light impacts are needed. Because the climate bias correction technique cannot be used in coupled simulations, the ESMs need to be carefully calibrated.
ZHANG ET AL.

Conclusions
In summary, we performed and reported here a set of factorial simulations using the ORCHIDEE LSM forced with climate fields from the IPSL-CM6A-LR ESM under different configurations. We calculated the differential impact of diffuse and direct radiation in the cloudy and clear gridbox fractions, in order to understand the impact of anthropogenic aerosols on terrestrial C fluxes during the period since 1850. We found that, by the year 2014, anthropogenic aerosols cumulatively enhanced land C sink by 22.6 PgC, out of which 78% is attributed to the impact of aerosol-induced radiation quality changes. This enhancement of radiation quality changes is much stronger than the dimming caused by aerosols. Although the cooling due to anthropogenic aerosols on average only slightly enhances the global land C sink, it has considerable impacts on its temporal evolution and at different latitudes. Compared with radiation and temperature, precipitation changes due to aerosols have smaller impacts on C fluxes in our study.
The diffuse radiation impact on land C sink detected in this study is of considerable magnitude. This large impact highlights the necessity of considering the radiation quality in simulating land C fluxes, especially for the projection of future climate under different aerosol emission scenarios, as we should remind that the climate-carbon feedback is one of the most important ones in the Earth system (Friedlingstein et al., 2006(Friedlingstein et al., , 2014. Unfortunately, the direct and diffuse radiative fluxes are only distinguished and adequately treated in a limited number of LSMs. There are almost no studies evaluating this process in current LSMs at large scale. Such studies are thus urgently needed to better quantify current aerosol impacts and project future C fluxes under different aerosol emission scenarios.

Appendix A: Radiation Partitioning and Diffuse Radiation Calculation at Subgrid Level
At each time step, each gridbox is first divided into a clear and a cloudy part according to the total cloud fraction (clt) of this gridbox (Figure 1). In the clear sky sub-grid area, the incoming direct (SWdown clear,dir ) and diffuse (SWdown clear,dif ) shortwave radiation at the surface are obtained considering tropospheric and stratospheric aerosol optical depth at 550 nm, as well as solar zenith angle from a look up table. The look up table is built with an atmospheric radiative transfer model named Streamer (Key & Schweiger, 1998) as initially described in Mercado et al. (2009). Using SWdown clear,dir and SWdown clear,dif we are able to calculate the sum of their counterparts in the cloudy sub-grid area (SWdown cloudy,dir + SWdown cloudy,dif ) from the average incoming downward shortwave radiation field (SWdown), which is provided in IPSL-CM6A-LR climate: where SWdown cloudy,dir and SWdown cloudy,dif represent the direct and diffuse radiative fluxes of the cloudy subgrid area. Here, we assume SWdown cloudy,dir to be zero following Mercado et al. (2009). This is a good assumption because the optical depth of clouds is generally larger than a few units except for thin high-level clouds. The SW cloudy,dif can then be estimated as: where SWdown and clt come from the forcing files while SWdown clear,dir and SWdown clear,dif come from the look up table. It should be noted that as ORCHIDEE_DF has a time step of 30 min but the forcing fields have coarser temporal resolution, we have to interpolate all forcing data into 30 min in ORCHIDEE_DF. Furthermore, the IPSL-CM6A-LR model, which has generated the forcing fields, and the look up table use different radiative transfer scheme and the IPSL-CM6A-LR radiation scheme relies on the full 3D structure of clouds. This causes unavoidable mismatch between the different terms of Equation A2, which may occasionally lead to unreasonable negative SWdown cloudy,dif or larger SWdown cloudy,dif than SWdownclear,dir + SWdown clear,dif . In the first case, the SWdown cloudy,dif is set as zero and SWdown clear,dir and SWdownclear,dif are linearly scaled to keep the total SWdown conserved. In the second case, SWdown cloudy,dif is capped ZHANG ET AL.