Aerosol optical depth thresholds as a tool to assess diffuse radiation fertilization of the land carbon uptake in China

China suffers from frequent haze pollution episodes that alter the surface solar radiation and influence regional carbon uptake by the land biosphere. Here, we apply combined vegetation and radiation modeling and multiple observational datasets to assess the radiative effects of aerosol pollution in China on the regional land carbon uptake for the 2009–2011 period. First, we assess the inherent sensitivity of China’s land biosphere to aerosol pollution by defining and calculating two thresholds of aerosol optical depth (AOD) at 550 nm, (i) AODt1, resulting in the maximum net primary productivity (NPP), and (ii) AODt2, such that if local AOD < AODt2, the aerosol diffuse fertilization effect (DFE) always promotes local NPP compared with aerosol-free conditions. Then, we apply the thresholds, satellite data, and interactive vegetation modeling to estimate current impacts of aerosol pollution on land ecosystems. In the northeast, observed AOD is 55 % lower than AODt1, indicating a strong aerosol DFE on local NPP. In the southeastern coastal regions, observed AOD is close to AODt1, suggesting that regional NPP is promoted by the current level of aerosol loading, but that further increases in AOD in this region will weaken the fertilization effects. The North China Plain experiences limited enhancement of NPP by aerosols because observed AOD is 77 % higher than AODt1 but 14 % lower than AODt2. Aerosols always inhibit regional NPP in the southwest because of the persistent high cloud coverage that already substantially reduces the total light availability there. Under clear-sky conditions, simulated NPP shows widespread increases of 20–60 % (35.0± 0.9 % on average) by aerosols. Under all-sky conditions, aerosol pollution has spatially contrasting opposite sign effects on NPP from−3 % to +6 % (1.6± 0.5 % on average), depending on the local AOD relative to the regional thresholds. Stringent aerosol pollution reductions motivated by public health concerns, especially in the North China Plain and the southwest, will help protect land ecosystem functioning in China and mitigate long-term global warming.


Introduction
Atmospheric aerosols scatter and absorb solar radiation, while plants rely on sunlight for photosynthesis.Thus, aerosols affect land carbon uptake through radiative perturbations.In particular, observations have demonstrated that aerosols can enhance canopy photosynthesis and light-use efficiency (LUE = GPP/PAR; GPP is gross primary productivity and PAR is photosynthetically active radiation) by increasing diffuse radiation in the lower canopy (Gu et al., 2003;Rap et al., 2015;Strada et al., 2015).This aerosol diffuse fertilization effect (DFE) is subject to the aerosol loading and sky conditions (Cohan et al., 2002;Oliphant et al., 2011) because the potential benefit of increased diffuse radiation in the lower canopy can be offset or even reversed by the concomitant reductions in direct sunlight.China is the world's largest anthropogenic emitter of carbon dioxide, reaching 2.5 Pg C yr −1 (Liu et al., 2015), while the land ecosystems mitigate only around 0.2 Pg C yr −1 with large uncertainties (Piao et al., 2009).At the same time, China experiences frequent haze pollution events due to large emissions of anthropogenic aerosols and precursors (Guo et al., 2014;Wang and Chen, 2016).It is critically important to understand how this haze pollution affects the land carbon sink in China.a Plant functional types (PFTs) include evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), trees (mixture of ENF/DBF/EBF and shrub), herbs (grass and crop), and savanna.b Carbon metrics include gross primary productivity (GPP), net primary productivity (NPP), and net ecosystem productivity (NEP).c Clearness index (CI) is converted to diffuse fraction (DF) with DF = 1.45-1.81× CI (Alton, 2008).
Leaf photosynthesis increases with the solar irradiance before reaching saturation (Farquhar et al., 1980).For a canopy with complex composition and vertical distribution, only sunlit leaves receive both direct and diffuse sunlight.Typically, irradiance is abundant for these leaves and photosynthesis is light-saturated.In contrast, shaded leaves receive only diffuse radiation and their photosynthesis is usually light-limited (He et al., 2013).Existence of aerosol pollution and/or a cloud layer simultaneously increases diffuse radiation but decreases direct radiation.The enhancement of diffuse radiation helps increase photosynthesis by the shaded leaves, but the response of the sunlit leaves depends on the level of aerosol/cloud loading.Moderate reductions in direct sunlight will not decrease photosynthesis of sunlit leaves because the light availability is still saturated.Consequently, the GPP of the whole canopy (sunlit plus shaded portions) increases due to the improved LUE (Knohl and Baldocchi, 2008).Large reductions in direct sunlight may convert the light-saturated regime into a light-limited regime for the sun-lit leaves, leading to reduced LUE and dampened canopy GPP (Alton, 2008).Photosynthetic response to diffuse light is also dependent on plant functional type (PFT).C4 plants are less sensitive to the enhanced diffuse radiation compared to C3 plants because C4 plants do not become light-saturated under high irradiance.As a result C4 plants are more sensitive to the reductions in direct light than C3 plants ( Kanniah et al., 2012).Thus, the net effect of aerosol pollution on canopy carbon uptake depends on the aerosol loading, cloud amount, geographic location, and plant species.
Previous observation-based studies of cloud and aerosol DFE are summarized in Table 1.Most observational studies have detected DFE using long-term ground-based measurements (Niyogi et al., 2004;Cirino et al., 2014).Currently, direct measurements of DFE on photosynthesis in China are limited (Li et al., 2014).Observations suggest that both clouds and aerosols exert similar impacts on land carbon uptake (Kanniah et al., 2012).Many observational studies have found that canopy GPP of trees maximizes with a diffuse fraction (DF) of 0.4-0.7 (Rocha et al., 2004;Alton, 2008).For grass and savanna, the optimal DF is as low as 0.2-0.3,above which the carbon uptake will decrease (Alton, 2008;Kanniah et al., 2013).The appearance of thin cloud may enhance the net ecosystem productivity (NEP) of trees by 7-11 % (Monson et al., 2002;Misson et al., 2005), while thick cloud reduces carbon uptake due to large irradiance attenuation (Rocha et al., 2004;Cheng et al., 2016).Aerosol light scattering in clear sky may increase the NEP of trees by 8-29 % (Misson et al., 2005;Cirino et al., 2014), but decreases the carbon uptake of grassland (Niyogi et al., 2004).Previous vegetation modeling results are generally consistent with observations (Table 1).For example, Knohl and Baldocchi (2008) predicted maximum GPP with a DF of 0.45 for a deciduous forest.
In this study, we apply the Yale Interactive terrestrial Biosphere (YIBs) model (Yue and Unger, 2015) combined with the Column Radiation Model (CRM) to analyze the impacts of aerosol DFE on the net primary productivity (NPP) in China in the present-day world.The main objective is to explore the responses of NPP to current aerosol pollution with and without the appearance of clouds.First, we perform multiple sensitivity experiments to derive the NPP sensitivity to aerosol optical depth (AOD) at 550 nm and compare it with available observations (Table 1).Second, we calculate the aerosol-induced DFE "tolerance" of China's land biosphere by defining and computing two thresholds of AOD, (i) AOD t1 , resulting in the maximum NPP, and (ii) AOD t2 , such that if local AOD < AOD t2 , the aerosol DFE always promotes local NPP compared with aerosol-free conditions.Third, we estimate changes in NPP between simulations with and without aerosol DFE, and relate these changes to the derived AOD thresholds so as to understand the causes of NPP responses to aerosol radiative effects.Our model approach offers a large regional-scale assessment, and is not limited by spatiotemporal and species-specific sampling issues (Table 1).We consider aerosol-induced perturbations to diffuse, direct, and total PAR under both clear-and cloudysky conditions, but ignore the meteorological and hydrological feedbacks from those perturbations.Section 2 describes the measurement data, vegetation and radiation models, and full group of simulations.Section 3 presents the model evaluation, the sensitivity of GPP to aerosol pollution in China, the derived AOD thresholds, and the simulated NPP responses to current levels of aerosol pollution.Section 4 summarizes and discusses the main results.

The Yale Interactive Terrestrial Biosphere Model (YIBs)
The YIBs is a process-based vegetation model that simulates the global terrestrial carbon cycle with dynamic pre-dictions of leaf area index (LAI) and tree growth (Yue and Unger, 2015).Plant photosynthesis is simulated using the well-established Farquhar scheme (Farquhar et al., 1980) and is coupled to stomatal conductance with the Ball-Berry scheme (Ball et al., 1987).The canopy radiative transfer scheme separates diffuse and direct PAR for sunlit and shaded leaves (Spitters, 1986), depending on solar zenith angle and canopy LAI (Sect.2.3).Autotrophic respiration (Ra) is split into maintenance and growth components, and is dependent on leaf temperature and nitrogen content (Clark et al., 2011).The model includes nine PFTs, including evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), shrubland, tundra, C3/C4 grassland, and C3/C4 cropland (Fig. S1 in the Supplement).This land cover is derived based on retrievals from both MODIS (Hansen et al., 2003) and the Advanced Very High Resolution Radiometer (AVHRR) (Defries et al., 2000).The fraction of C4 cropland is derived based on the total crop fraction and the ratio of C4 species (Monfreda et al., 2008).
The YIBs can be used in three different configurations: site level, global/regional offline, and online within a climate model.For this study, we use the regional offline version driven with hourly 1 • × 1 • meteorological forcings from the NASA Modern Era Retrospective analysis for Research and Applications (MERRA).On the global scale, simulated LAI, tree height, phenology, GPP, and NPP show reasonable spatial distribution and long-term trends compared with both in situ measurements and satellite retrievals (Yue and Unger, 2015;Yue et al., 2015a, b).Other carbon fluxes such as NEP, autotrophic respiration, and heterotrophic respiration are also reasonably simulated relative to multiple model ensembles (Yue et al., 2015b).

The Column Radiation Model (CRM)
The CRM is the standalone version of the radiation model used by the NCAR Community Climate Model, which has been updated to the Community Earth System Model (http://www.cesm.ucar.edu/models/).Using temporally varying aerosol profiles (types and concentrations) and meteorological reanalyses, the CRM calculates reflectivity and transmission of atmospheric layers, emissivity and absorptivity of greenhouse gases (GHGs), and Mie scattering and absorption of aerosols.Aerosol optical parameters associated with each aerosol species, including specific extinction coefficients, single scattering albedo, and asymmetric parameters, are adopted from Yue et al. (2010) for mineral dust, Yue and Liao (2012) for sea salt, and the RegCM4 model for other species (Giorgi et al., 2012).These parameters vary with changes in both wavelength and relative humidity (Fig. S2 in the Supplement).Sulfate and nitrate aerosols share the same parameters.For carbonaceous aerosols (black carbon, BC, and organic carbon, OC), half are considered hydrophobic and half are hydrophilic.
The CRM is driven with hourly 1 • × 1 • fields of temperature, humidity, and [O 3 ] at 20 sigma levels interpolated from the MERRA data.Vertical profiles of cloud cover and liquid water path are adopted from the CERES SYN1deg (http: //ceres.larc.nasa.gov),which are determined using remote sensing from MODIS and the Visible and Infrared Sounder (VIRS).Surface albedo, temperature, and pressure are also adopted from MERRA.The CRM utilizes aerosol profiles of all species simulated by the ModelE2-YIBs, a fully coupled chemistry-carbon-climate model (Schmidt et al., 2014).Aerosol components include sulfate, nitrate, BC, OC, dust (clay and silt), and sea salt (accumulation and coarse modes).Concentrations of these pollution species are predicted based on emission inventories of the year 2010 from the Greenhouse Gas and Air Pollution Interactions and Synergies (GAINS) integrated assessment model (Amann et al., 2011).We compare the GAINS v4a inventory with the HTAP inventory adopted from the Emissions Database for Global Atmospheric Research (EDGAR, http://edgar.jrc.ec.europa.eu)and the RCP8.5 inventory from the Intergovernmental Panel on Climate Change (IPCC, http://www.ipcc.ch/)(Fig. S3 in the Supplement).The intercomparison shows that GAINS has similar magnitudes (differences within ±10 %) of emissions for major pollutants over China to HTAP and RCP8.5, except for ammonia, which is higher by 50-80 % in GAINS.Simulated summertime surface PM 2.5 concentrations show high correlations (R = 0.76) and low biases (NMB = 1.6 %) with in situ measurements at 188 sites (not shown).

Canopy radiative transfer and carbon uptake
We use the multilayer canopy radiative transfer scheme proposed by Spitters (1986) to separate diffuse and direct PAR for sunlit and shaded leaves.The canopy is divided into an adaptive number of layers (typically 2-16) for light stratification.Light intensity decreases exponentially with LAI when penetrating into the canopy: where I t is the total PAR at the top of the canopy, L is the LAI from the top of the canopy to layer n, I is the total PAR available for absorption at depth L, and k is the extinction coefficient.Here, ρ is the reflection coefficient calculated as follows: where α is the solar zenith and σ = 0.2 is the scattering coefficient of single leaves.Light absorption at depth L is estimated as follows: where I a is the flux absorbed per unit leaf area.The total PAR at the top of the canopy consists of diffuse (I tf ) and direct (I tr ) We amplify or diminish base AOD by a certain ratio for each simulation.The base AOD (Fig. 3a) is derived with the aerosol profiles from the ModelE2-YIBs climate model and optical parameters from multiple data sources (Fig. S2 in the Supplement).
b For clear-sky simulations, cloud cover and liquid water path are set to zero.c For all-sky simulations, cloud cover and liquid water path are adopted from the Clouds and the Earth's Radiant Energy System (CERES) SYN1deg product.
components, both of which are simulated with the CRM: According to Eq. ( 3), absorption of the diffuse flux is calculated as where k f = 0.8(1 − σ ) 1/2 is the extinction coefficient of the diffuse flux.For the direct flux, it is segregated into diffuse and direct components on its way through the canopy.The total absorption of the direct flux is calculated as where k r = 0.5/ cos α is the extinction coefficient of the direct component of the direct flux.Here the function (1−σ ) 1/2 is applied to account for the scattering effects of leaves for direct light.The absorption of the direct component of the direct flux is calculated as We distinguish between light absorption for shaded and sunlit leaves.Shaded leaves absorb diffuse flux and the diffuse component of direct flux: Sunlit leaves absorb both diffuse and direct radiation: Photosynthesis is calculated separately for shaded and sunlit leaves based on the different light absorption.Canopy photosynthesis (µmol C s −1 per unit leaf area) is the sum of these two parts of leaves: where A sl and A sh are the photosynthesis of sunlit and shaded leaves, respectively.The fraction of sunlit leaf area f sl is calculated as Finally, the total carbon uptake GPP (µmol C s −1 per unit ground area) is estimated as follows:

Simulations
We conduct simulations combining the offline YIBs vegetation model and the CRM.Diffuse and direct PAR at the top of the canopy is predicted with the CRM.These radiative fluxes are then used as input for the YIBs model, which further separates diffuse and direct fluxes absorbed by sunlit and shaded leaves using the Spitters (1986) canopy scheme.We perform two groups of YIBs-CRM sensitivity simulations, 30 for clear-sky conditions and 30 for all-sky conditions, to derive the AOD thresholds for aerosol DFE (Table 2).The YIBs model is driven with meteorological forcings from MERRA, except for direct and diffuse PAR, which is predicted with the CRM.We set up baseline simulations (CLR010 and ALL010) with the aerosol profile from ModelE2-YIBs and validated optical parameters.In other simulations, specific scaling factors ranging from 0.0 to 30 are applied to aerosol concentrations to represent variations in AOD.Due to the disk storage limit for hourly meteorological profiles, we perform CRM simulations for 2009-2011.The simulated PAR is alternately applied as input for the YIBs model, which uses additional meteorological forcings from MERRA for 2000-2011.In this case, the predicted PAR at 2009-2011 is recycled as input for the periods of 2000-2002, 2003-2005, 2006-2008, and 2009-2011 in the YIBs simulations.The last 3 years of YIBs output, including GPP and NPP, are used for analyses.We focus our analyses for the summer (June-July-August) season, when both AOD and carbon fluxes are high.

Benchmark and evaluation observational datasets
To evaluate GPP simulated with the YIBs model, we use global benchmark product upscaled from the FLUXNET eddy covariance data for 2009-2011 (Jung et al., 2009).
For NPP, we use the satellite product of 2009-2011 retrieved by the Moderate Resolution Imaging Spectroradiometer (MODIS) (Zhao et al., 2005).The MODIS dataset provides indirect estimates of global NPP using an empirical light-use function between GPP and meteorology, as well as the modeled plant respiration.The actual values of MODIS NPP may exhibit certain biases at the regional scale compared with site-level observations (Pan et al., 2006).However, the spatial pattern of the product is in general reasonable and has been widely used for model evaluations (e.g., Collins et al., 2011;Pavlick et al., 2013).To evaluate surface radiative fluxes simulated with the CRM, we use groundbased radiation data for 2008-2012 from 106 pyranometer sites in China, provided by the Climate Data Center, Chinese Meteorological Administration (Xia, 2010).For each site and each month, we derive the monthly mean radiative fluxes based on daily data if < 30 % days are missing.We then select sites where all months are available for 2008-2012, leaving a total of 95 sites for the evaluation of total shortwave radiation.Diffuse radiation is not observed at all sites, and only a total of 17 sites meet the criteria for the continuous measurements.For aerosol radiative effects, we use assimilation data of surface radiative fluxes adopted from the SYN1deg product of NASA Clouds and the Earth's Radiant Energy System (CERES) (Wielicki et al., 1996;Rutan et al., 2015).Aerosol effect in CERES is calculated with the Langley Fu-Liou radiative transfer model (Fu and Liou, 1993), using aerosol profiles simulated by the Model for Atmospheric Transport and Chemistry (Rasch et al., 1997).We also use satellite-based AOD data of MOD08_M3 for 2008-2012 retrieved by MODIS onboard the Terra platform (Remer et al., 2005).All gridded data are interpolated onto 1 • × 1 • grids, matching the resolution of both CRM and YIBs models.

Model evaluation
The YIBs model predicts reasonable spatial distribution of carbon fluxes compared with other data products (Fig. 1).For the summer, high GPP and NPP are simulated in the northeast, the southwest, and the southeastern coastal regions, where DBF and ENF trees dominate (Fig. S1 in the Supplement).The correlation coefficients between simulations and proxy data are as high as 0.8 for both GPP and NPP.Predicted GPP is −24 % lower on average than the benchmark product, mainly because the model shows smaller values over the North China Plain.The normalized mean bias (NMB) of NPP is close to 0, because of the regional off-  setting.On the annual mean basis, simulations show higher correlations (R = 0.84 for GPP and 0.86 for NPP) and similar NMB (−21 % for GPP and −16 % for NPP) compared with data products.
The CRM predicts opposite spatial distribution for total solar radiation and the corresponding DF (Fig. 2).For radiation, high values are found in the north and west, and low values in the east and southwest.Such a pattern is related to cloud cover, which is lower in the north but higher in the south, especially the southwest, where cloud cover is usually higher than 80 % (Fig. S4 in the Supplement).Due to the cloud scattering, those regions with low insolation have high DF.Compared with in situ measurements, the simulated total radiation exhibits reasonable spatial characteristics and the correlation coefficient is as high as 0.88.The evaluation of DF shows certain biases but is reasonable over the east (blue points of Fig. 2f), which is the major domain for this study.
The CRM also predicts reasonable AOD and aerosol radiative effects (Fig. 3).Using aerosol concentrations from the ModelE2-YIBs climate model, the CRM simulates high regional AOD centered in the North China Plain.Previous sensitivity tests with ModelE2-YIBs show that such high loading of aerosols is mainly (> 80 %) contributed by anthropogenic emissions.The AOD is generally high over the vast domain of eastern China but low in the western part.Compared with AOD from CERES, which is derived using aerosol concentrations from a chemical transport model, the AOD in CRM presents reasonable spatial distribution with high correlation coefficient of 0.7 and low NMB of −1 %.However, compared with MODIS AOD, which is derived based on satellite retrievals, the predicted AOD is on average overestimated by 30 % in eastern China.For the following analyses, we use the predicted AOD as the benchmark but discuss the influence of its overestimation on the predicted aerosol DFE.We further assess aerosol radiative efficiency (ARE), defined as radiative forcing per unit AOD, over China.A higher magnitude of ARE is predicted in the north and west, because lower cloud coverage there allows larger radiative perturbations by the same level of aerosols.shown are for clear-sky (red empty points) and all-sky (blue solid points) conditions.Separately for clear-sky and all-sky conditions, we first collect all grid cells in eastern China (box region in Fig. 2a) for all 30 sensitivity simulations.We then aggregate all grid cells with non-zero fraction of a specific PFT into 30 AOD bins ranging from 0 to 6 at an interval of 0.2.In each bin, we calculate average diffuse fraction and corresponding GPP, with an error bar indicating one standard deviation.
The comparison of ARE between CRM and CERES shows a high correlation coefficient of 0.87.However, ARE in CRM is smaller by 21 % in magnitude than that in CERES.
X. Yue and N. Unger: Aerosol optical depth thresholds

Sensitivity of GPP to DF and AOD in China
Appearance of cloud and/or aerosols increases DF but decreases total insolation.We examine PFT-specific GPP responses to DF for clear and all-sky conditions (Fig. 4).Under clear-sky conditions, at DF < 0.55, all PFTs show increased GPP in response to increasing DF.At low DF (DF < 0.55), shaded leaves experience low light availability because diffuse radiation is limited.Meanwhile, photosynthesis of sunlit leaves is light-saturated because direct radiation is abundant.Introduction of aerosol pollution, which increases DF, to this system redistributes sunlight to the shaded light-limited leaves (without compromising the total light availability to sunlit leaves), thus increasing the LUE and GPP of the whole canopy.The derived GPP-diffuse sensitivity ( GPP/ DF, units: g C m −2 day −1 per unit change in DF with ±95 % confidence interval) is 5.4 ± 1.2 for ENF, 6.8 ± 2.1 for DBF, 5.6 ± 1.7 for shrub (tundra plus arid shrub), 2.8 ± 0.7 for C3 herbs (C3 grassland and cropland), and 2.2 ± 2.3 for C4 herbs (C4 grassland and cropland) when DF < 0.55.However, at high DF (DF > 0.55), light is no longer saturated for sunlit leaves because of the large attenuation of direct light.Further introduction of aerosols decreases photosynthesis of sunlit leaves, which may offset the carbon gains from enhanced diffuse light by shaded leaves, resulting in a net carbon loss for the whole canopy.Although large variations in GPP, denoted as error bars, exist within each bin, they do not affect the significance of the GPP-diffuse responses.We select the GPP-diffuse sensitivity under clear-sky conditions averaged for all PFTs as an example.The largest binto-bin difference in GPP is calculated as 1.9 g C m −2 day −1 between DF = 0.58 and DF = 0.12.Such GPP difference is significant at p < 0.001 level based on a Student's T test, suggesting that GPP varies significantly when the DF change is pronounced.
Under all-sky conditions, which includes both clear and cloudy skies, DF is always higher than 0.55, because existing cloud cover has already increased the diffuse light fraction in the system.Any further increase in DF by aerosol scattering has limited or even detrimental impacts on whole canopy GPP because the associated aerosol-induced reductions in direct radiation impact photosynthesis by the sunlit leaves.GPP decreases almost linearly for all PFTs in response to increasing DF > 0.55, with GPP-diffuse sensitivity of −6.9 ± 1.4 for ENF, −10.8 ± 2.3 for DBF, −3.8 ± 1.7 for shrubland, −3.8 ± 0.8 for C3 herbs, and −10.1 ± 1.6 for C4 herbs.The GPP response to increases in DF > 0.55 is almost identical under clear-sky and all-sky conditions.
The PFT-specific GPP responses to idealized perturbations in AOD depend strongly on the sky conditions (Fig. 5).Under clear-sky conditions, aerosol promotes GPP for most PFTs if AOD < 2. The maximum possible enhancement of GPP is ∼ 40 % for DBF, ENF, and C3 herbs, and can be as high as 60 % for shrub.Most shrub species (especially tundra) are located in the southwest (Fig. S1 in the Supplement).
Over that region, solar irradiance is abundant at clear days (not shown), allowing more efficient scattering for a given AOD.For a given AOD, the DFE of C4 herbs is the weakest amongst PFTs with a maximum possible GPP enhancement of only ∼ 10 %.C4 plants usually have lower LUE than C3 plants, and as a result more easily become light-limited when aerosol attenuates total irradiance (Still et al., 2009;Kanniah et al., 2012).A clear turning point is found for all species.For C3 plants, aerosol scattering weakens GPP when AOD > 2, because photosynthesis of the sunlit leaves starts to become light-limited due to reductions in direct insolation.For C4 plants, this turning point appears earlier when AOD > 1.Under all-sky conditions, atmospheric aerosol shows neutral effects on GPP when AOD < 1 and detrimental negative effects when AOD > 1 for all PFTs except C4 (Fig. 4).For C4 plants, any addition of aerosol to the atmospheric column decreases GPP.
Our estimates of GPP sensitivity to DF and AOD are reasonable compared with previous studies (Table 1) as summarized below: (1) the maximum enhancement of GPP is 40 % under clear-sky conditions for most PFTs (Gu et al., 2003); (2) the GPP enhancement of C4 plants is the least due to the lowest LUE compared with other PFTs (Still et al., 2009;Kanniah et al., 2012); (3) the aerosol DFE is much stronger under clear-sky than under all-sky conditions (Cohan et al., 2002); (4) both clouds and aerosols exert similar DFE on land carbon uptake (Kanniah et al., 2012;Cirino et al., 2014); and (5) the maximum GPP enhancement appears when DF = 0.4-0.8(Rocha et al., 2004;Alton, 2008;Zhang et al., 2010).Most results listed in Table 1 are based on NEP; however, we evaluate the sensitivity of GPP because it is the direct carbon metric affected by DFE.

Aerosol pollution DFE in China 2009-2011
We apply the idealized GPP responses in Sect.3.2 to estimate the current magnitude of aerosol pollution DFE under realistic background conditions by defining 2 summertime AOD thresholds across China (Fig. 6).The AOD thresholds are derived based on NPP (= GPP − Ra), assuming no impacts of aerosol DFE on the autotrophic respiration Ra (Knohl and Baldocchi, 2008).Now, we consider NPP responses, instead of GPP, because the former represents the net carbon uptake by plants after subtracting autotrophic respiration for maintenance and growth.The first threshold, AOD t1 is defined as the AOD value with the maximum NPP (the uppermost point of the response curve in Fig. 5), representing the saturation of DFE.The second threshold, AOD t2 is the cross point of the response curve at zero NPP.For each model grid cell, if AOD is close or equal to AOD t1 , the peak NPP is expected.If AOD < AOD t2 , aerosol DFE always promotes local NPP relative to aerosol-free conditions (AOD = 0).
We find relatively high AOD t1 in the northeast and low values in the southwest (23-35 • N, 100-105 • E, box d in Fig. 6a), where average cloud cover is 80 % in the summer shown are different for clear-sky (red empty points) and all-sky (blue solid points) conditions in summer (June-August).Separately for clear-sky and all-sky conditions, we first collect all grid cells in eastern China (box region in Fig. 2a) for all 30 sensitivity simulations (Table 2).We then aggregate all grid cells with non-zero fraction of the specific PFT into 30 AOD bins ranging from 0 to 6 at an interval of 0.2.In each bin, we calculate average GPP change relative to aerosol-free conditions, with an error bar indicating one standard deviation.
(Fig. S4 in the Supplement).The pattern of AOD t2 is similar to that of AOD t1 (Fig. 6b), except for high values in the north (42-48 • N, 105-118 • E), where insolation is high, while average cloud fraction is less than 50 %.The values of the AOD thresholds are much higher for clear days; for example, the average clear-sky AOD t1 is 6 times the all-sky values over the east (Fig. S5 in the Supplement).Observed summer AOD in the North China Plain (32-40 • N, 113-120 • E) exceeds AOD t1 by 77 % but is 14 % lower than AOD t2 (Fig. 6c), suggesting limited aerosol DFE in this region.A reduction of 44 % in local AOD (so that AOD = AOD t1 ) leads to the largest benefit for regional carbon uptake.Both AOD t1 and AOD t2 in the southwest are close to 0 (Fig. 6d), indicating that appearance of aerosol always inhibits regional carbon uptake there.Observed AOD on the southeastern coast (23-30 • N, • E) is approximately equal to AOD t1 (Fig. 6e) and that in the northeast (40-47 • N, 122-132 • E) is 55 % lower than AOD t1 (Fig. 6f), indicating amplified carbon uptake by aerosol enhancements there.
To understand the relationship between the current AOD level and DFE, we calculate differences between MODIS AOD and the thresholds (Fig. 7).The pattern is quite similar for annual and summer differences, except that the former is less positive than the latter.The stronger dampening effect by aerosols in summer is related to the higher seasonal AOD (summer mean of 0.43 versus annual mean of 0.38) and cloud amount (summer mean of 65 % versus annual mean of 58 %).Over the southwest, the current AOD level is larger than both AOD t1 and AOD t2 , suggesting that aerosols there on average inhibit NPP and further increases in AOD lead to stronger inhibition.In the North China Plain, AOD exceeds AOD t1 and is close to AOD t2 , indicating that aerosol pollution there has almost neutralized impacts on local NPP.For the southeast, current AOD is lower than AOD t1 and AOD t2 at confined regions along the coast, but is higher than AOD t1 in the inner domain.On average, observed AOD of the box domain e (Fig. 6a) is lower than AOD t2 but close to AOD t1 (Fig. 6e), suggesting that aerosols there generally promote NPP relative to aerosol-free conditions.However, the potential for stronger DFE is limited as further increases in AOD will dampen NPP.For the north and northeast, the current AOD is far below AOD t1 , suggesting that the appearance of aerosols there boosts carbon uptake, and that increases in AOD continue to increase NPP.
We calculate perturbations in NPP for different sky conditions (Fig. 8), resulting from the aerosol-induced perturba-  tions in both direct and diffuse radiation (Fig. S6 in the Supplement).Under clear-sky conditions, aerosol DFE causes widespread enhancement of NPP (Fig. 8a), ranging from 20 to 60 % over the east (Fig. S7 in the Supplement).The absolute changes in NPP are small over the North China Plain (Fig. 8a), where high AOD is observed (Fig. 3).In contrast, fractional changes in NPP exhibit a high center of > 60 % in the North China Plain, consistent with the conclusion of sensitivity tests that aerosols usually promote carbon uptake in clear sky (Fig. 5).Such discrepancy originates from the low background (aerosol-free) NPP in the North China Plain, because the YIBs model applies satellite-based land cover (Fig. S1 in the Supplement), which shows a high fraction of C3 cropland but almost zero tree (including ENF and DBF) coverage in the North China Plain.On average, aerosols increase NPP by 1.14 ± 0.01 Pg C yr −1 (35.0 ± 0.9 %) for the whole China domain under the clear-sky conditions (Fig. S7 in the Supplement).
Under the all-sky conditions, aerosols drive weak patchy NPP responses (Fig. 8b), mainly because of the high DF from existing cloud cover.The spatial pattern of percentage NPP changes (Fig. S7b in the Supplement) highly resembles the AOD differences shown in Fig. 7d but with opposite signs.Over the southwest and part of the North China Plain, the current high level of AOD exceeds AOD t2 and as a result inhibits local NPP by 1-2 %.In the southeastern coastal regions, aerosol DFE is limited, though regional AOD is below AOD t2 .The largest NPP enhancement is predicted over the northeast, where current AOD is far smaller than AOD t1 (Fig. 7c) and cloud amount is moderate (Fig. S4 in the Supplement).Regional changes in NPP range from −3 to 6 % and the total change is 0.07 ± 0.02 Pg C yr −1 (1.6 ± 0.5 %) over China.

Discussion and conclusions
We provide the first assessment of aerosol pollution radiative effects on land carbon uptake in China today based on regional simulations that combine the CRM radiation model and the YIBs vegetation model.The confidence level of our estimate is dependent on the capability of these models to reproduce observed radiative and carbon fluxes, aerosolinduced radiative perturbations, and GPP responses to these perturbations.For the first two aspects, we evaluate CRM and YIBs models using in situ, satellite, and assimilation data (Figs.1-3).The simulated GPP-AOD relationship (Figs.4-5) is reasonable compared with available measurements and modeling results in the literature (Table 1).
Our estimate is subject to limitations and uncertainties.First, calculated aerosol DFE is sensitive to the canopy radiative transfer scheme.We apply the widely used Spitters (1986) scheme to separate diffuse and direct light for sunlit and shaded leaves (two-leaf mode).This scheme invokes Beer's law that assumes light decays exponentially from the top to bottom of canopy.The predicted maximum GPP enhancement of 40 % is at the high end of the range of previous modeling estimates (Table 1).A similar magnitude of GPP change was predicted by Gu et al. (2003) using different parameterizations of light partitioning and leaf photosynthesis (one-leaf mode).Mercado et al. (2009) predicted maximum GPP enhancement of only 18 % for deciduous trees, considering Beer's law for light extinction and a two-leaf model for light partitioning.Cohan et al. (2002) showed a maximum NPP enhancement of 17-30 % depending on the choice of canopy scheme.These discrepancies reveal large uncertainties due to differences in the treatment of the canopy geometry (sphere or non-sphere), canopy properties (e.g., leaf clumping factor, leaf inclination angle, and leaf optical properties), light partitioning (diffuse or non-diffuse), and the upscaling of leaf photosynthesis (one-or two-leaf).More observations are required to evaluate the different parameterizations their use in large-scale vegetation models.
Second, uncertainties in simulated AOD and aerosol radiative effects may affect the derived aerosol DFE.For this study, simulated AOD is calculated based on the threedimensional aerosol concentrations from ModelE2-YIBs climate model.Simulation shows similar spatial pattern and magnitude compared with CERES product (Fig. 3a-c), which also uses aerosol profile from a chemical transport model.However, evaluations with MODIS data show that the predicted AOD is on average overestimated by 30 % in eastern China (not shown).Considering that aerosol radiative efficiency from the CRM is 21 % lower than that from the Fu-Liou model (Fig. 3d-f), our estimate of aerosol radiative perturbations in China might be reasonable due to the offsetting biases in AOD and aerosol radiative efficiency.On the other hand, even if we ignore the uncertainties from the radiative transfer scheme, the 30 % overestimation in AOD does not cause large differences in the derived aerosol DFE.As a check, we calculate NPP in sensitivity experiments CLR007 and ALL007 (Table 2), which employ an AOD level lower by 30 % than CLR010 and ALL010.We find that clear-sky NPP by aerosols is 0.91 Pg C yr −1 (27.7 %) in CLR007, lower than the enhancement of 35.0 % in CLR010 (Fig. 8a).The all-sky NPP by aerosols is 0.07 Pg C yr −1 (1.6 %) in ALL007, equal to the values derived from ALL010.The reason for such similarity between ALL007 and ALL010 can be explained by the GPP-AOD relationship, which shows that all-sky GPP is almost invariant when AOD < 1 (Fig. 5).The results suggest that cloud plays a dominant role in regulating diffuse radiation in China, and the aerosol DFE might be secondary compared with cloud DFE.
Third, calculated aerosol DFE does not account for biotic feedbacks.Photosynthesis is connected to plant physiological processes, such as stomatal conductance and respiration.Observations have shown that aerosol DFE may increase water use efficiency (Rocha et al., 2004;Knohl and Baldocchi, 2008), promoting plant growth and LAI that may further increase canopy photosynthesis.We have assumed no responses in autotrophic respiration so that the derived GPP-AOD relationship can be applied directly to NPP-AOD.However, observations suggest that plant respiration decreases due to the aerosol-induced cooling (Alton, 2008).Ignoring these biotic feedbacks indicates that NPP sensitivity to AOD employed in our estimate may be underestimated.
Fourth, we omit the associated climatic responses to aerosol radiative effects.The aerosol-induced radiative perturbations decrease surface temperature but increase relative humidity (because of decreased saturation water vapor pressure) (Jing et al., 2010;Cirino et al., 2014).Increases in air humidity will help enhance plant water use efficiency, leading to increased photosynthesis.The impact of cooling is uncertain, depending on the relationship between the local background temperature and the optimal temperature of photosynthesis, which is about 25 • C for most species (Farquhar et al., 1980).If leaf temperature is > 25 • C, aerosol-induced cooling is beneficial for photosynthesis.By contrast, if leaf temperature is < 25 • C, the cooling will act to inhibit carbon uptake.Furthermore, cloud modification, caused by both aerosol direct and indirect effects, may exert complex influences on land ecosystems through perturbations in diffuse radiation, surface temperature, and precipitation.Resolving these concomitant biotic, meteorological, and hydrological feedbacks requires earth system models that fully couple the land biosphere, atmospheric chemistry, radiation, and climate.
Finally, the biogeochemical response to aerosol pollutionrelated reactive nitrogen (N) deposition is not included.Simulations with the ModelE2-YIBs climate model show that anthropogenic emissions contribute 90 % of the total reactive inorganic and organic N deposition in China, indicating potentially large impacts of anthropogenic aerosols on regional carbon uptake through the coupled C-N cycle.Using a terrestrial ecosystem model, Tian et al. (2011) proposed that anthropogenic N deposition and fertilizer applications together account for 61 % of the net carbon storage in China for the 1961-2005 period.The carbon sequestered per gram of deposited nitrogen decreases gradually after the year 1985, suggesting that most areas have reached nitrogen saturation.Similarly, based on satellite retrievals, Xiao et al. (2015) revealed that anthropogenic N deposition makes no significant contributions to the increases in vegetation productivity during 1982-2006.Therefore, additional fertilization from aerosol N deposition may be limited because many areas have been N saturated for decades such that our estimate of NPP due to aerosol effects based on radiative changes (Fig. 8b) may be realistic.
Despite these uncertainties, our study reveals strong impacts of aerosol DFE on land carbon uptake in China.Although aerosol DFE widely promotes NPP during clear days, NPP shows strongly spatially contrasting responses under all-sky conditions.Aerosol pollution increases NPP by 2-6 % in the northeast, where both cloud coverage and particle loading are moderate.Aerosol decreases NPP by 2-4 % in both the North China Plain and the southwest.For the North China Plain, the NPP inhibition is related to the high local pollution level that is above the AOD t2 threshold for carbon uptake.Our estimates show that a 44 % reduction in local aerosol AOD would help achieve the maximum benefits for plant productivity in this region.In the southwest, existing cloud cover is already dense, and pollution aerosols inhibit NPP in this region.Reductions in pollution aerosol loading will increase carbon uptake in this region.

Data availability
Data used in this study can be provided upon request to the corresponding author Xu Yue (xuyueseas@gmail.com).
The Supplement related to this article is available online at doi:10.5194/acp-17-1329-2017-supplement.

Figure 1 .
Figure 1.Evaluation of summertime carbon fluxes simulated with the YIBs model.Simulations are (a) GPP and (d) NPP from ALL010.Derived carbon fluxes are averaged for 2009-2011, with (b) GPP from the benchmark upscaling of flux tower data and (e) NPP from MODIS.The aerosol DFE is included in the simulation.The correlation coefficients (r), relative biases (b), and number of 1 • × 1 • grid cells (n) for the comparisons are listed in the scatter plots (c and f).The dashed line represents the 1 : 1 ratio.The red line is the linear regression between simulations (predictand) and data products (predictor).The unit of GPP and NPP is g C m −2 day −1 .

Figure 2 .
Figure 2. Evaluation of summertime radiation fluxes simulated with the CRM.Simulations are (a) surface total shortwave radiation (W m −2 ) and (d) diffuse radiation fraction from ALL010 with 1 • × 1 • resolution.Observations (b and e) are the average during 2008-2012 from 106 sites operated by the Climate Data Center, Chinese Meteorological Administration.The correlation coefficients (r), relative biases (b), and number of sites (n) are shown in the scatter plots (c and f).The blue points in the scatter plots represent sites located within the box regions in eastern China as shown in (a).The dashed line represents the 1 : 1 ratio.The red line is the linear regression between simulations and observations.

Figure 3 .
Figure 3. Evaluation of aerosol radiative effects simulated with the CRM.Panels show the simulated (a) AOD at 550 nm and (d) all-sky aerosol radiative efficiency (ARE, W m −2 per unit AOD) with that from the CERES SYN1deg product (b, e).The correlation coefficients (r), relative biases (b), and number of 1 • × 1 • grid cells (n) for the comparisons over the box domain in panel (a) are listed in the scatter plots (c, f).The dashed line represents the 1 : 1 ratio.The red line is the linear regression between simulations and data products.

Figure 4 .
Figure 4. Mean summer GPP for different diffuse fractions.Resultsshown are for clear-sky (red empty points) and all-sky (blue solid points) conditions.Separately for clear-sky and all-sky conditions, we first collect all grid cells in eastern China (box region in Fig.2a) for all 30 sensitivity simulations.We then aggregate all grid cells with non-zero fraction of a specific PFT into 30 AOD bins ranging from 0 to 6 at an interval of 0.2.In each bin, we calculate average diffuse fraction and corresponding GPP, with an error bar indicating one standard deviation.

Figure 5 .
Figure 5. Sensitivity of summer GPP to changes in AOD.Resultsshown are different for clear-sky (red empty points) and all-sky (blue solid points) conditions in summer (June-August).Separately for clear-sky and all-sky conditions, we first collect all grid cells in eastern China (box region in Fig.2a) for all 30 sensitivity simulations (Table2).We then aggregate all grid cells with non-zero fraction of the specific PFT into 30 AOD bins ranging from 0 to 6 at an interval of 0.2.In each bin, we calculate average GPP change relative to aerosol-free conditions, with an error bar indicating one standard deviation.

Figure 6 .
Figure 6.Thresholds of AOD at 550 nm for aerosol DFE.(a) AOD t1 threshold leading to maximum NPP in summer (June-August).(b) AOD t2 threshold leading to enhanced summer NPP.Aerosol indirect effects and climatic feedback are not included for these thresholds.The average AOD-NPP response curves at four box domains in (a) are shown in (c)-(f), with red symbols indicating two AOD thresholds and blue symbols indicating observed AOD from MODIS.Chinese regions with low aerosol-free NPP (< 0.05 g C m −2 day −1 ) are blanked in (a) and (b).The color scales for (a) and (b) are different.

Figure 7 .
Figure 7. Differences between MODIS AOD and two derived AOD thresholds.Left panel shows AOD between MODIS and AOD t1 for the average of (a) whole year and (c) summer months.Right panel shows AOD between MODIS and AOD t2 for the average of (b) whole year and (d) summer months.For the left panel, regions with positive values indicate that increase (decrease) in local AOD leads to reductions (enhancement) in NPP, while regions with negative values denote that decrease (increase) in local AOD results in reductions (enhancement) in NPP.For the right panel, regions with negative (positive) values indicate current level of AOD always promotes (inhibits) NPP, compared with the aerosol-free conditions.The color scales among panels are different.

Figure 8 .
Figure 8. Changes in summer NPP caused by aerosol DFE in China for (a) clear-sky (CLR010-CLR000) and (b) all-sky (ALL010-ALL000) conditions.The percentage changes in NPP are shown in Fig. S7 in the Supplement.The color scales between panels are different.Units: g C m −2 day −1 .

Table 1 .
Summary of studies about diffuse fertilization effect (DFE).