Pathway dependence of ecosystem responses in China to 1.5 °C global warming

China is currently the world's largest emitter of both CO₂ and short-lived air pollutants. Ecosystems in China help mitigate a part of the country's carbon emissions, but they are subject to perturbations in CO₂, climate, and air pollution. Here, we use a dynamic vegetation model and data from three model inter-comparison projects to examine ecosystem responses in China under different emission pathways towards the 1.5 ^∘C warming target set by the Paris Agreement. At 1.5 ^∘C warming, gross primary productivity (GPP) increases by 15.5±5.4 % in a stabilized pathway and 11.9±4.4 % in a transient pathway. CO₂ fertilization is the dominant driver of GPP enhancement and climate change is the main source of uncertainties. However, differences in ozone and aerosols explain the GPP differences between pathways at 1.5 ^∘C warming. Although the land carbon sink is weakened by 17.4±19.6 % in the stabilized pathway, the ecosystems mitigate 10.6±1.4 % of national emissions in the stabilized pathway, more efficient than the fraction of 6.3±0.8 % in the transient pathway. To achieve the 1.5 ^∘C warming target, our analysis suggests a higher allowable carbon budget for China under a stabilized pathway with reduced emissions in both CO₂ and air pollutants.

1 Introduction

The past decade has seen record-breaking warming largely related to anthropogenic greenhouse gas emissions (Mann et al., 2017). This warming trend presents a challenge to achieving the temperature control target of 1.5 ^∘C above the pre-industrial (PI) level set by the 2015 Paris Agreement. Many studies have shown that a conservative level of warming such as 1.5 ^∘C is necessary to limit climatic extremes (Nangombe et al., 2018), avoid heat-related mortality (Mitchell et al., 2018), reduce economic loss (Burke et al., 2018), and alleviate ecosystem risks (Warszawski et al., 2013) compared to stronger anthropogenic warming. To achieve this target, each country must aim to control its greenhouse gas emissions. A full understanding of regional ecosystem response to the changing climate and environmental stress is essential to reduce uncertainties in allowable carbon budget estimates at 1.5 ^∘C (Mengis et al., 2018). China is covered with a wide range of terrestrial biomes (Fang et al., 2012). While China's ecosystem response to possible future climate has been explored (Wu et al., 2009; He et al., 2017; Dai et al., 2016), impacts on the regional carbon budget of differing pathways to the 1.5 ^∘C target are not known.

There are two distinct pathways to the 1.5 ^∘C global warming. One is a fast process in which global temperature passes 1.5 ^∘C and continues to increase (scenarios assuming high CO₂ emissions and no climate mitigation) while the other is a stabilized process with an equilibrium warming right below 1.5 ^∘C that lasts for decades before the end of 21st century (scenarios including climate mitigation; James et al., 2017). The stabilized pathway is the one proposed by the 2015 Paris Agreement. However, the unprecedented warming in 2016 results in an increase of global average temperature by 1.1 ^∘C above PI (https://public.wmo.int/, last access: 27 February 2019), suggesting that the 1.5 ^∘C limit may be broken in the near future under a transient pathway (Henley and King, 2017). A few studies have compared allowable carbon budgets between these two pathways (Collins et al., 2018; Millar et al., 2017), but none has estimated the mitigation potential of regional ecosystems with joint impacts of changes in climate, CO₂, and air pollution under different pathways.

Here, we apply the Yale Interactive terrestrial Biosphere (YIBs) model (Yue and Unger, 2015, 2018) to investigate the response of terrestrial ecosystem productivity in China to both stabilized and transient global warming of 1.5 ^∘C relative to the PI period. We focus on the changes of gross primary productivity (GPP) and net ecosystem exchange (NEE). GPP represents total canopy photosynthesis through gross carbon assimilation. NEE is the residue after subtraction of GPP from ecosystem (plant + soil) respiration (Reco − GPP), indicating the net carbon sink from land to atmosphere. The larger the GPP values, the stronger the carbon assimilation by ecosystems. In contrast, the more negative the NEE, the stronger the carbon sink of the land. The YIBs model is driven with meteorology from an ensemble of climate models in Climate Model Intercomparison Project Phase 5 (CMIP5). The stabilized global warming pathway is represented by the RCP2.6 low-emissions scenario that yields an equilibrium change in global mean temperature (ΔGMT) of 1.49 ^∘C by 2050–2070 with selected climate models (Fig. S1 in the Supplement). The transient pathway is represented by RCP8.5 high-emission scenario in which ΔGMT grows rapidly and realizes a transient 1.5 ^∘C around the years 2021–2041. We select the present-day period of 1995–2015 as a reference.

2 Methods

2.1 Datasets

2.1.1 CMIP5 data

We use both daily and monthly meteorology predicted by CMIP5 models (https://esgf-node.llnl.gov/, last access: 1 February 2020). The daily data are used as input for YIBs. In total, we select 15 climate models (Table S1 in the Supplement) with all available daily meteorology, including surface air temperature, precipitation, specific humidity, surface downward shortwave radiation, surface pressure, and surface wind speed, for a historical scenario and two future scenarios (RCP2.6 and RCP8.5). These two scenarios assume distinct emission pathways of both CO₂ and air pollutants, with the RCP2.6 scenario projecting much lower CO₂ and pollution concentrations than RCP8.5. Simulated annual GMT is smoothed with a 21-year window to remove decadal variations. The ensemble changes of GMT relative to PI period (1861–1900) from two scenarios are examined (Fig. S1a). The low-emission scenario RCP2.6 yields an equilibrium ΔGMT of 1.85 ^∘C by 2100. We remove eight climate models predicting stabilized ΔGMT higher than 1.85 ^∘C by the end of century. The seven remaining models yield an ensemble warming close to 1.5 ^∘C (1.49 ^∘C for 2050–2070, Fig. S1b). Meanwhile, ΔGMT in the high-emission scenario RCP8.5 grows fast and realizes a transient 1.5 ^∘C warming around the years 2021–2041. Daily meteorology from seven selected models (Table S1) are then interpolated to the uniform 1^∘ × 1^∘ resolution and used to drive YIBs to simulate terrestrial carbon fluxes in China for 1850–2100. Due to the large data storage, we retain only the domain of (15–60^∘ N, 60–150^∘ E) covering Chinese territory. We bias correct modeled meteorology with WFDEI (WATCH Forcing Data methodology applied to ERA-Interim reanalysis) data (Weedon et al., 2014):

$$\begin{array}{}\text{(1)}& V_{\mathrm{d}}^{\mathrm{s}}=V_{\mathrm{d}}\times S_{\mathrm{w}}/S_{\mathrm{m}}.\end{array}$$

Here V_d is the original daily variables and $V_{\mathrm{d}}^{\mathrm{s}}$ is the scaled value. S_w is the two-dimensional WFDEI value averaged for 1980–2004 and S_m is the modeled values averaged for the same period. In this case, the average climate from each individual model matches observations at present day; meanwhile, climate variability from models are retained to estimate uncertainties in carbon fluxes.

2.1.2 TRENDY-v6 data

We acquire the global GPP and NEE datasets from 1901 to 2016 simulated by 14 dynamic global vegetation models (DGVMs) participating in the TRENDY project (Table S2). All DGVMs are implemented following the same simulation protocol and are driven by consistent input datasets, including CRU-NCEP climate data, and atmospheric CO₂ concentrations, but fixed present-day land use (Le Quéré et al., 2018).

2.1.3 ACCMIP O₃ data

We use monthly output of surface O₃ concentrations from 12 models joining the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP, Lamarque et al., 2013; Table S3). The ACCMIP models have a wide range of horizontal and vertical resolutions, natural emissions, chemistry schemes, and interaction with radiation and clouds. However, these models apply the same anthropogenic and biomass burning emissions specified for CMIP5 RCP scenarios (e.g., RCP2.6 or RCP8.5), though different models perform simulations at different time slices. Here, we use surface O₃ and interpolate the original output to 1^∘ × 1^∘ resolution. We fill the temporal gaps between two adjacent time slices using a linear fitting approach. In this way, we derive the monthly O₃ from 1850 to 2100 for each model and their ensemble average at each grid point.

2.1.4 Diffuse radiation data

The original CMIP5 archive does not provide a diffuse component of shortwave radiation. Here, we use empirical relations between total and diffuse radiation from 11 studies to calculate hourly diffuse radiation (Table S4). The diffuse fraction k_d in all equations depends on clearness index k_t, which is defined as the ratio between global solar radiation I_t and extra-terrestrial solar radiation I₀ (Ghosh et al., 2017):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{t}}=I_{\mathrm{t}}/I_{\mathrm{0}},\text{(3)}& {\displaystyle }& {\displaystyle }{I}_{\mathrm{0}}=I_{\mathrm{sc}}\left[\mathrm{1}+\mathrm{0.033}\cos \left({\displaystyle \frac{\mathrm{360}N}{\mathrm{365}}}\right)\right]\cos \mathit{\phi }.\end{array}$$

Here I_sc=1367 W m⁻² is solar constant, N is Julian day of the year, and φ is solar zenith. The empirical equations are evaluated using hourly total and diffuse radiation from Modern-Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al., 2011) during 2008–2012. For each grid in China, we calculate hourly diffuse radiation (D_c) using MERRA total radiation and compare it with the standard output (D_m). Statistical metrics including correlation, normalized mean bias (NMB), and normalized root mean square error (NRMSE) are used to evaluate the performance of empirical equations:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{NMB}=\left(\stackrel{\mathrm{‾}}{D_{\mathrm{c}}}-\stackrel{\mathrm{‾}}{D_{\mathrm{m}}}\right)/\stackrel{\mathrm{‾}}{D_{\mathrm{m}}},\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{NRMSE}=\sqrt{\sum {\displaystyle \frac{{\left(D_{\mathrm{c}}-D_{\mathrm{m}}\right)}^{\mathrm{2}}}{n}}}/\stackrel{\mathrm{‾}}{D_{\mathrm{m}}}.\end{array}$$

Here $\stackrel{\mathrm{‾}}{D_{\mathrm{c}}}$ and $\stackrel{\mathrm{‾}}{D_{\mathrm{m}}}$ are mean values of calculated and MERRA diffuse radiation, respectively. The evaluation is performed month by month for 2008–2012 and n is the number of daytime samples (grids with total radiation > 5 W m⁻²). The value of n varies from month to month with a minimum of 540 000 in December 2010. Evaluation shows that the empirical model M01 (Lam and Li, 1996) yields the highest correlation and the lowest NRMSE (Fig. S2). As a result, we use the M01 model to derive diffuse radiation from CMIP5 models.

2.2 Model

We apply the YIBs model (Yue and Unger, 2015; Yue et al., 2017) to simulate historical and future (1850–2100) ecosystem productivity. The YIBs model dynamically calculates leaf area index (LAI) and tree height based on carbon assimilation and allocation. Leaf-level photosynthesis is calculated hourly using the well-established Farquhar et al. (1980) scheme and is upscaled to canopy level by the separation of sunlit and shaded leaves (Spitters, 1986). Sunlit leaves can receive both direct and diffuse radiation, while shaded leaves receive only the diffuse component (Yue and Unger, 2017). The assimilated carbon is in part used for maintenance and growth respiration, and the rest is allocated among leaf, stem, and root for plant growth (Clark et al., 2011). Soil respiration is calculated as the loss of carbon flows among 12 soil carbon pools (Schaefer et al., 2008). The YIBs model considers nine plant functional types (PFTs) including evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), evergreen broadleaf forest (EBF), shrubland, tundra, C₃ grassland, C₄ grassland, C₃ cropland, and C₄ cropland. The land cover is prescribed based on satellite retrievals from the Moderate Resolution Imaging Spectroradiometer (MODIS; Hansen et al., 2003) and the Advanced Very High Resolution Radiometer (AVHRR; Defries et al., 2000). For this study, we fix the land cover to isolate impacts of CO₂ and climatic changes. Other studies also show only moderate changes in vegetation fraction and composition at a low warming level (Warszawski et al., 2013). The YIBs model can be applied at the site, regional, and global scales. The site-level model has been evaluated with measured carbon fluxes from 145 FLUXNET sites (Yue and Unger, 2015). For this study, all simulations are performed at the 1^∘ × 1^∘ resolution over China. During the period of 1982–2011, YIBs predicts an average GPP of 7.17 Pg C yr⁻¹ in China (Fig. S3), close to the 7.25 Pg C yr⁻¹ estimated in the benchmark product (Jung et al., 2009).

YIBs calculates O₃ damage to plant photosynthesis using a flux-based parameterization (Sitch et al., 2007). The inhibition rate of GPP is dependent on both ambient O₃ concentrations and stomatal conductance. Compared to meta-analyses of hundreds of data sources from China (Table S5) and the world (Yue and Unger, 2018), the scheme shows good performance in estimating GPP responses to O₃ for DBF, EBF, C₃, and C₄ herbs (Fig. S4). The predicted O₃ damaging effects on ENF might be underestimated. The YIBs model separates the effects of diffuse and direct light on plant photosynthesis (Spitters, 1986). Simulated GPP responses to direct and diffuse radiation show good agreement with observations at 24 global flux tower sites from FLUXNET network (Yue and Unger, 2018). In general, diffuse radiation is more efficient at enhancing canopy photosynthesis compared to the same level of direct radiation.

2.3 Simulations

We perform two main groups of simulations, one for RCP2.6 and the other for RCP8.5. For each group, seven sub-groups are designed with varied climatic or CO₂ forcings (Table S6). In each sub-group, separate runs are conducted for the YIBs model driven with climate variables from seven selected CMIP5 models (Table S1), making a total of 98 runs. A baseline group simulation (HIST_2000) is performed with fixed meteorology and CO₂ after the year 2000. Another four sub-group simulations are performed to quantify O₃ effects on photosynthesis (Table S7). These simulations are driven with both CMIP5 meteorology and monthly O₃ concentrations from an ensemble of 12 ACCMIP models. The runs are distinguished with different O₃ damaging sensitivity (high or low) and scenario projections (RCP2.6 or RCP8.5). Monthly O₃ concentrations are downscaled to an hourly time step using the diurnal cycle simulated by a chemistry–climate model NASA ModelE2 (Schmidt et al., 2014). The O₃-affected GPP or NEE are calculated as the average of simulations with low and high sensitivities.

For each run, a 251-year simulation is performed with the historical climate for 1850–2000 and the future climate for 2001–2100. For simulations driven with meteorology from the same climate model, all sensitivity tests apply the same climate forcing during the historical period but utilize varied forcings after the year 2000. For example, RCP26_CO2 is identical to RCP26_MET for the period of 1850–2000. However, after the year 2000, the former runs fixed climatic conditions for the year 2000 but allows changes in CO₂ concentrations year by year for 2001–2100 following the pathway projection, while the latter fixes the CO₂ level for the year 2000 but continues to use day-to-day meteorology after 2000. For all simulations, we initialize vegetation and soil carbon pools in the YIBs model with a 200-year spin-up by recycling meteorology for the year 1850. Contributions of individual factors are calculated as the differences between sensitivity and baseline group (e.g., RCP26_CO2–HIST_2000 for CO₂ fertilization in RCP2.6 scenario).

The main focus of this study is to quantify how the differences of anthropogenic emissions, including both CO₂ and air pollution, which are usually associated, will cause different responses in the land carbon budget to the same global warming target. We highlight that the role of air pollution on the land carbon cycle has always been ignored. The assumptions of land use can be quite uncertain between future pathways (Stehfest et al., 2019), and these assumptions are not necessarily associated with CO₂ and air pollution emissions. As a result, for this study, we consider fixed land cover in all simulations. We do not consider ecosystem–climate interactions as in earlier studies (e.g., Yue et al., 2017), so as to take advantage of climate projections from multiple models.

3 Results

3.1 Changes of atmospheric compositions and radiation

The ensemble concentrations of ACCMIP O₃ show good agreement with ground-based observations from 1580 sites in China (Fig. 1). The spatial correlation is R=0.80 (p<0.01) between observations and the ensemble O₃ concentrations ([O₃]), though the latter is higher by 25 % (Fig. 1a–c). Such overestimation can likely be attributed to the high [O₃] at night in the models, because the evaluation of maximum daily 8 h average (MDA8) [O₃], which mainly occurs in the daytime, shows more reasonable predictions with a lower bias of 10 % (Fig. 1d–f). Since O₃ vegetation damage in general occurs in the daytime, when both plant photosynthesis and [O₃] are at high levels, the ACCMIP [O₃] can be used as input for YIBs to derive long-term O₃ inhibition effects on ecosystem productivity.

Figure 1 Evaluation of surface O₃ with site-level observations. Simulations are ensemble (a) mean and (d) daily maximum 8 h average (MDA8) O₃ for the period of 2005–2015 from 12 ACCMIP models. Observations (b, e) are the average during 2015–2018 from 1580 sites operated by Ministry of Ecology and Environment, China. The correlation coefficients (r), relative biases (b), and number of sites (n, excluding sites lacking data) are shown in the scatter plots (c, 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.

The ensemble radiation from CMIP5 models matches observations at 106 sites in China (Fig. 2). For total shortwave radiation, the model prediction shows high values in the west and low values in the southeast, consistent with observations for a correlation coefficient of R=0.79 (p<0.01) and a mean bias of 8.9 %. The derived diffuse radiation is highest in the southeast, where the total radiation is lowest. Observed diffuse radiation is available only at 17 sites. Compared to these sites, predictions show reasonable spatial distribution with a correlation of R=0.65 (p<0.01) and a low bias of 7.1 %. Both the total radiation and derived diffuse radiation are used as input for YIBs to estimate GPP responses to joint changes in direct and diffuse radiation caused by aerosol removal.

Figure 2 Evaluation of radiation fluxes with site-level observations. Simulations are surface (a) total shortwave radiation (W m⁻²) and (d) diffuse radiation derived with method M01 (Table S4) for the period of 2005–2015 from an ensemble of seven CMIP5 climate models. Observations (b, e) are the average during 2009–2011 from 106 sites operated by the Climate Data Center, Chinese Meteorological Administration. The correlation coefficients (r), relative biases (b), and number of sites (n, excluding sites lacking data) are shown in the scatter plots (c, 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.

Atmospheric compositions and radiation show varied changes in different scenarios. The GMT changes mainly follow those in CO₂ concentrations, which show fast growth in RCP8.5 but slow changes in RCP2.6 (Fig. 3a). The latter assumes a large reduction of carbon emissions globally after the year 2020 (Meinshausen et al., 2011). Global CO₂ levels reduce slightly after the year 2030 in RCP2.6, while GMT continues growing until 2050 due to air–sea interactions (Solomon et al., 2009). As a low-emission scenario, RCP2.6 experiences a slow growth in nitrogen oxide (NO_x) emissions and a continuous reduction after the year 2020 (Fig. S5), resulting in a decline of 6.4 ppb (15.2 %) in surface O₃ over eastern China under 1.5 ^∘C warming at 2060 (Fig. 3b). In contrast, RCP8.5 assumes fast growth of NO_x emissions with delayed controls after the year 2030, leading to surface O₃ enhancements of 6.6 ppb (15.7 %) under 1.5 ^∘C warming at 2030. The lower emissions in RCP2.6 also result in smaller aerosol optical depth (AOD) than RCP8.5 (Fig. S6), leading to higher surface total radiation (Fig. 3c), with lower diffuse radiation (Fig. 3d) due to reducing light extinction (Yu et al., 2006).

Figure 3 Changes in atmospheric compositions and radiation. Results shown are projected future (a) global CO₂ concentrations, and (b) surface O₃ concentrations, (c) total photosynthetically active radiation (PAR), and (d) diffuse PAR at growth season in China. The average (a) CO₂ concentrations at the global warming of 1.5 ^∘C are 442 ppm for RCP2.6 scenario (blue, 2050–2070) and 454 ppm for RCP8.5 scenario (red, 2021–2041). The (b) O₃ concentrations are averaged over east of 110^∘ E in China from 12 ACCMIP models for RCP2.6 (blue) and RCP8.5 (red) scenarios. Each dot represents the value averaged for May to September from a chemistry model. The (c, d) PAR values are averaged over China from seven CMIP5 models for RCP2.6 (blue) and RCP8.5 (red) scenarios. Diffuse PAR is calculated using hourly total PAR and solar zenith angle based on the parameterization M01. Each dot represents the value averaged for May to September from a climate model. For each selected year in (b–d), a period of 11 years (5 years before and 5 years after) is used to derive the decadal mean values.

3.2 Historical ecosystem productivity in China

The ensemble simulations show an increasing trend in GPP in China of 0.011 Pg C yr⁻² over the historical period 1901–2016 (Fig. 4a). A stronger trend of 0.022 Pg C yr⁻² is found after 1960. Such a change is much faster than the trend of 0.013 Pg C yr⁻² estimated by a benchmark product (Jung et al., 2009) for 1982–2011 but close to a recent estimate of 0.02 Pg C yr⁻² from combined machine learning algorithms and eddy flux measurements from 40 sites in China (Yao et al., 2018). The simulated trend is also consistent with the TRENDY ensemble, which predicts trends of 0.013±0.006 Pg C yr⁻² (ensemble ± inter-model uncertainty) for 1901–2016 and 0.022±0.01 Pg C yr⁻² for 1961–2016. The YIBs simulations show variabilities of 0.41±0.23 Pg C yr⁻¹ (6.2±3.9 %, blue shading in Fig. 4a) due to uncertainties in climate from CMIP5 models, much smaller than the value of 1.33±0.16 Pg C yr⁻¹ (19.2±2.6 %, red shading in Fig. 4a) caused by structural uncertainties across different vegetation models.

Figure 4 Historical carbon fluxes in China. Results shown are simulated (a) gross primary productivity (GPP) and (b) net ecosystem exchange (NEE) during the historical period (1901–2016) using YIBs (blue), and the comparison with predictions of 14 terrestrial models from TRENDY project (red). The bold lines are ensemble means with red shadings for inter-vegetation-model uncertainties and blue shadings for inter-climate-model uncertainties. All YIBs simulations are driven with daily meteorology from CMIP5 models. All TRENDY simulations are driven with CRU-NCEP meteorology. The black line in (a) represents benchmark results of 1980–2011 from Jung et al. (2009). The black point with an error bar in (b) represents the synthesis of ground- and model-based estimates of NEE in China by Piao et al. (2009).

NEE in China is negative, suggesting a regional land carbon sink (Fig. 4b). This sink is −94.7 Tg C yr⁻¹ with a trend of −1.7 Tg C yr⁻² during 1901–2016. Such change matches TRENDY simulations, which predict a multi-model mean carbon sink of $-\mathrm{74.1}\pm \mathrm{30.8}$ Tg C yr⁻¹ (uncertainties due to inter-model variations) and a trend of $-\mathrm{1.3}\pm \mathrm{0.7}$ Tg C yr⁻² for the same period. During 1980–1989, the ground-based estimate (Piao et al., 2009) suggests a sink of 177±73 Tg C yr⁻¹ in China, consistent with the sink intensity of 149±20 Tg C yr⁻¹ from the YIBs ensemble prediction. For the recent period of 1980–2000, YIBs estimates a strengthened sink of 154±30 Tg C yr⁻¹ in China, weaker than the estimate of 198±114 Tg C yr⁻¹ with the Dynamic Land Ecosystem Model (DLEM; Tian et al., 2011) but within the estimates of 137–177 Tg C yr⁻¹ based on both ground and satellite data (Fang et al., 2007). The interannual variability in YIBs simulations is much weaker than the estimates in other studies, because the ensemble approach largely dampen variations among different runs. Similar to GPP, the NEE simulations exhibit smaller variability of 62±50 Tg C yr⁻¹ among different YIBs runs than that of 122±57 Tg C yr⁻¹ among different TRENDY models.

Figure 5 Changes in carbon fluxes by global warming of 1.5 ^∘C. Results shown are simulated (a–c) GPP and (d–f) NEE over China between the period of global warming of 1.5 ^∘C and present day (1995–2015) under (a, d) RCP2.6 scenario, (b, e) RCP8.5 scenario, and (c, f) their differences. The period of global warming of 1.5 ^∘C is set to 2050–2070 for RCP2.6 and 2021–2041 for RCP8.5. Simulations are performed using the YIBs vegetation model driven with daily meteorology from seven CMIP5 models. The O₃ damaging effect is included with predicted ensemble O₃ concentrations from 12 ACCMIP models. For each grid, significant changes at p<0.05 are marked with dots. The total changes (Pg C yr⁻¹) over China are shown in each panel.

3.3 Future changes of carbon fluxes

Projected GPP continues to increase in both RCP2.6 and RCP8.5 scenarios after the year 2016 (Fig. S7a). With global warming of 1.5 ^∘C, GPP increases significantly in China, especially over eastern and northeastern parts (Fig. 5). Compared to the present day, GPP with O₃ effects increases by 1.07±0.38 Pg C yr⁻¹ (15.5±5.4 %) in the RCP2.6 scenario (Fig. 5a) and 0.82±0.30 Pg C yr⁻¹ (11.9±5.4 %) in RCP8.5 (Fig. 5b). The spatial pattern of the GPP changes is similar in the two pathways (correlation coefficient R=0.93), except that ΔGPP in RCP2.6 is higher than in RCP8.5 by 30 %, with a positive center over eastern China (Fig. 5c). Projected NEE continues to be more negative in the RCP8.5 scenario after the year 2016 (Fig. S7b). Meanwhile, future NEE reaches the minimum value (or the maximum sink strength) around the year 2025 and then reverses to be less negative in the RCP2.6 scenario (Fig. S7b). By the time of 1.5 ^∘C global warming, NEE changes in China show opposite tendencies between the two pathways. Compared to the present day, NEE increases by 0.03±0.03 Pg C yr⁻¹ ($-\mathrm{17.4}\pm \mathrm{19.6}$ %) in RCP2.6 (Fig. 5d) but decreases by 0.14±0.04 Pg C yr⁻¹ (94.4±24.9 %) in RCP8.5 (Fig. 5e), suggesting that the land carbon sink is slightly weakened in the former but strengthened in the latter. Their differences exhibit widespread positive values in China with high centers in the east (Fig. 5f).

The changes in carbon fluxes follow the variations in atmospheric composition and climate (Figs. 6 and S8–S11). With global warming of 1.5 ^∘C, a dominant fraction of GPP enhancement in China is attributed to CO₂ fertilization (Fig. 6a). For the RCP2.6 scenario, CO₂ alone contributes 0.83 Pg C yr⁻¹ (77 %) to ΔGPP, with the highest enhancement of 0.8 g C m⁻² d⁻¹ over the southeast coast (Fig. S8a). For RCP8.5, CO₂ fertilization increases GPP by 0.95 Pg C yr⁻¹, even higher than the total ΔGPP of 0.82 Pg C yr⁻¹. The larger CO₂-induced ΔGPP in RCP8.5 is due to the higher CO₂ concentrations (454 ppm) than RCP2.6 (442 ppm) at the same 1.5 ^∘C warming (Fig. 3a). The 12 ppm differences in CO₂ concentrations lead to a change of 0.12 Pg C yr⁻¹ (1.7 %) in GPP. This sensitivity of GPP to CO₂, 0.14 % ppm⁻¹, falls within the range of 0.05–0.21 % ppm⁻¹ as predicted by 10 terrestrial models (Piao et al., 2013) and that of 0.01 %–0.32 % ppm⁻¹ as observed from multiple free-air CO₂ enrichment (FACE) sites (Ainsworth and Long, 2005). The higher ΔGPP in RCP2.6 instead yields a weakened NEE (more positive) due to the CO₂ effects (Fig. 6b). The stabilization of CO₂ concentrations in this scenario (Fig. 3a) results in a stabilized GPP after the year 2040 (Fig. S7a). Meanwhile, the 55-year (from 2005 to 2060) carbon accumulation enhances soil carbon storage by 10.5±1.3 Pg C and promotes soil respiration to 0.71±0.19 Pg C yr⁻¹. The stabilized GPP with enhanced soil respiration (NEE = Reco − GPP; Reco includes both soil and plant respiration) together lead to a weakened carbon sink (less negative NEE) in the 1.5 ^∘C warming period (Fig. 7b). In contrast, soil carbon storage increases only 5.2±0.5 Pg C in RCP8.5 due to relatively short time period (from 2005 to 2031) for carbon accumulation, leading to lower soil respiration of 0.41±0.15 Pg C yr⁻¹ in the fast warming pathway. The continuous increase of GPP and lower soil respiration jointly strengthen the land carbon sink (more negative NEE) in China by 0.1 Pg C yr⁻¹ under the RCP8.5 scenario (Fig. 6a).

Figure 6 Attribution of changes in GPP and NEE to individual driving factors. Results shown are the predicted GPP changes in China between the period of global warming of 1.5 ^∘C and present day (1995–2015) caused by all (ALL) or individual driving factors, including CO₂ fertilization, O₃ damage, and meteorological changes (MET). The perturbations from meteorology is a combination of those of temperature (T), radiation (RAD), specific humidity (Q), and soil moisture (SOILM). The contrast is shown between the scenarios of RCP2.6 (blue, 2050–2070) and RCP8.5 (red, 2021–2041). The error bars indicate uncertainties of YIBs simulations using different future meteorology from seven CMIP5 models.

Figure 7 Damaging effects of O₃ to photosynthesis and plant growth. Results shown are ensemble mean changes in (a–c) offline GPP, (d–f) online GPP, and (g–i) leaf area index (LAI) caused by O₃ in the (a, d, g) present day (1995–2015) and 1.5 ^∘C warming under (b, e, h) RCP.6 (2050–2070) and (c, f, i) RCP8.5 (2021–2041) scenarios. The simulations are performed with the YIBs vegetation model driven with meteorology from 7 CMIP5 models and hourly ozone derived from 12 ACCMIP models. The damaging effect is averaged for high and low O₃ sensitivities. For each grid, significant changes at p<0.05 are marked with dots. The mean changes over China are shown in each panel.

Ozone (O₃) damages plant photosynthesis and the land carbon sink (Sitch et al., 2007; Yue and Unger, 2018). In the present day, O₃ decreases GPP by 6.7±2.6 % (uncertainties ranging from low to high damaging sensitivities) in China (Fig. 7d), because of the direct inhibition of photosynthesis by 6±2.4 % (Fig. 7a) and the consequent reduction of 1.8±0.8 % in leaf area index (LAI, Fig. 7g). For 1.5 ^∘C global warming, this weakening effect shows opposite tendencies in the two RCP scenarios, with a reduced GPP loss of 4.7±2.0 % in RCP2.6 (Fig. 7e) but an increased loss of 7.9±3.0 % in RCP8.5 (Fig. 7f). These impacts are predominantly driven by the variations of surface O₃ concentrations in the two scenarios, as predicted O₃ at 1.5 ^∘C warming decreases by 15.2 % in the low-emission pathway but increases by 15.7 % in the high-emission pathway (Fig. 3b). Consequently, changes in O₃ help increase GPP by 0.1±0.03 Pg C yr⁻¹ in RCP2.6 but decrease GPP by 0.14±0.04 Pg C yr⁻¹ in RCP8.5 for the same 1.5 ^∘C warming. Following the benefits to GPP, the lower O₃ decreases NEE (strengthens the sink) by 0.06±0.02 Pg C yr⁻¹ in RCP2.6, offsetting more than half of the negative effect (weakening the sink) from CO₂ (Fig. 6b). For RCP8.5, O₃ impacts make limited contributions to ΔNEE.

Changes in meteorology account for the rest of the perturbations in the carbon fluxes. With global warming of 1.5 ^∘C, the temperature in China increases by 0.90 ^∘C for RCP2.6 and 0.91 ^∘C for RCP8.5 (Fig. S12a and b) compared to the present-day climate. The spatial pattern of these changes is very similar without significant differences (Fig. S12c), leading to almost identical GPP responses (Figs. S8d and S9d). Generally, higher temperature is not beneficial for plant photosynthesis at low latitudes (Piao et al., 2013), where regional summer climate is already warmer than the optimal temperature threshold for leaf photosynthesis (Corlett, 2011). As a result, warming leads to negative changes in GPP over the east. Surface specific humidity exhibits widespread enhancement in eastern China (Fig. S13a and b). Air humidity may rise in a warmer climate because the corresponding enhancement of saturation pressure allows the atmosphere to hold more water vapor. On average, surface specific humidity increases by 0.34 g kg⁻¹ in RCP2.6 and 0.31 g kg⁻¹ in RCP8.5, leading to a promotion of GPP by 0.14 Pg C yr⁻¹ in the former and a similar value of 0.12 Pg C yr⁻¹ in the latter (Figs. S8e and S9e). Precipitation increases by 0.14 mm d⁻¹ (4.6 %) over eastern China in RCP2.6 but decreases by 0.03 mm d⁻¹ (1.2 %) in RCP8.5 (Fig. S12d and e), leading to higher soil moisture in eastern China for RCP2.6 (Fig. S13d and e). Nevertheless, most vegetation in eastern China is not water stressed, leaving moderate GPP responses to soil moisture changes in both RCP scenarios (Figs. S8f and S9f).

For the RCP2.6 scenario, the net effect of climate change causes an increase of 0.15 Pg C yr⁻¹ in GPP with a range from −0.54 to 0.62 Pg C yr⁻¹ (Fig. 6a). Such large variability in ΔGPP is related to the uncertainties in meteorology from different climate models. For RCP8.5, climate-induced GPP change is only 0.04 Pg C yr⁻¹ with a range from −0.6 to 0.26 Pg C yr⁻¹. The discrepancy of ΔGPP for the two pathways is mainly caused by the different radiation impacts, which enhance GPP by 0.2 Pg C yr⁻¹ in RCP2.6 but only 0.11 Pg C yr⁻¹ in RCP8.5 (Fig. 6a). Photosynthetically active radiation (PAR) is higher by 2.8 W m⁻² in RCP2.6 than in RCP8.5 (Fig. 3c). The distinct changes in radiation are related to aerosol radiative effects, because global analyses also show radiation enhancement in regions (e.g., USA and Europe) with aerosol removal (Fig. S14). The lower AOD in RCP2.6 helps increase solar insolation at the surface by reducing light extinction (Yu et al., 2006), and promote precipitation with weaker aerosol semi-direct and indirect effects (Lohmann and Feichter, 2005). Although lower aerosols in RCP2.6 slightly decrease diffuse radiation (Fig. 3d), which is more efficient in increasing photosynthesis (Mercado et al., 2009; Yue and Unger, 2018), the overall enhancement in total radiation helps boost GPP. Climate-induced ΔNEE is −0.02 Pg C yr⁻¹ (strengthened sink) for both pathways (Fig. 6b), resulting from comparable responses of NEE to changes in radiation (R=0.82), temperature (R=0.71), air humidity (R=0.91), and soil moisture (R=0.73) between the two pathways (Figs. S10 and S11).

3.4 Impacts on allowable carbon budget

For a warming target of 1.5 ^∘C, our analyses suggest that a simultaneous reduction of CO₂ and air pollution emissions enhances the efficiency of land carbon uptake compared to a pathway without air pollution emission control. The increased light availability from aerosol removal and decreased surface O₃ jointly promote GPP in China by 0.3 Pg C yr⁻¹, equivalent to 36 % of the CO₂ fertilization. In contrast, air pollution results in a net GPP inhibition of 0.03 Pg C yr⁻¹ under the high-emission pathway, suggesting a detrimental environment for plant health. Compared to RCP8.5, the timing of 1.5 ^∘C warming is delayed by 30 years in RCP2.6, leading to weaker carbon sink in the latter. However, even with the longer period of accumulation, the total carbon loss by O₃ damage is smaller by 3 %–16 % in RCP2.6 relative to RCP8.5 at the same warming level (Fig. 8a).

Figure 8 Accumulated carbon budget in China under 1.5 ^∘C global warming. (a) shows the total carbon loss of ecosystems caused by O₃ damaging effects at different warming thresholds for two emission pathways. The bottom panel shows the accumulated net carbon sink by ecosystems in China at the 1.5 ^∘C global warming. The two solid lines represent emissions of CO₂ and pollutants from the same scenario, either RCP2.6 (blue) or RCP8.5 (red). The dashed lines represent sensitivity experiments with inconsistent CO₂ and pollutants, with the blue (red) line driven with CO₂ from RCP2.6 (RCP8.5) but air pollution from RCP8.5 (RCP2.6). The warming of 1.0 ^∘C is the year 2010 for both RCP2.6 and RCP8.5 scenarios.

The slow warming increases the allowable cumulative anthropogenic carbon emissions. Assuming China's carbon emission fraction of 27 % of the world (the level at year 2017; Le Quéré et al., 2018), the total national emissions allowed are 80.4 Pg C in RCP2.6 and 71.9 Pg C in RCP8.5 from the year 2010 to the 1.5 ^∘C warming, following the global emission rates defined for these scenarios. The ensemble simulations show that ecosystems in China help mitigate 8.5±1.1 Pg C in RCP2.6 and 4.5±0.6 Pg C in RCP8.5 (Fig. 8b). Sensitivity experiments with either reduced CO₂ (but retaining high pollution) or reduced pollution (but retaining high CO₂) reveal land carbon uptakes of 7.3±0.9 and 5.0±0.6 Pg C, respectively. These values are both lower than that in RCP2.6, suggesting that simultaneous control of carbon and air pollution emissions can maximize the mitigation potential of ecosystems. The higher ecosystem assimilation rate in a low-emission pathway (10.6±1.4 % in RCP2.6 vs. 6.3±0.8 % in RCP8.5) over China, which is not considered in CMIP5 models, further buffers the pace of warming to 1.5 ^∘C.

4 Discussion and conclusions

Projection of future ecosystem productivity is subject to uncertainties in climate forcing and biophysical responses. The multi-model ensemble is a good approach to reduce the uncertainty in climate (Flato et al., 2013). In this study, we employ daily meteorology from seven CMIP5 models. A comparison with more CMIP5 models is performed (not shown) and confirms that the changes in meteorology from the seven selected climate models are robust and representative of future projections. As for ecosystem responses, future projections generally showed increasing GPP in China (Mu et al., 2008; Ji et al., 2008; Ju et al., 2007); however, climate change alone usually reduces productivity by inducing hot and drought conditions. In contrast, the YIBs simulations reveal a net positive effect of climate change on GPP though with large uncertainties (Fig. 6a). Such discrepancies are related to structural uncertainties across different vegetation models. Evaluations suggest that biophysical responses to environmental forcings in the YIBs model are generally reasonable as compared to the TRENDY ensemble (Fig. 4).

The YIBs simulations do not consider the nitrogen cycle and its limitation on carbon uptake. Inter-model comparisons show that models without nutrient constraints tend to overestimate GPP responses to CO₂ fertilization (Smith et al., 2016). As a result, the difference of CO₂ contributions in RCP scenarios would be smaller than predicted (Fig. 6a), suggesting that GPP enhancement in RCP2.6 might be even higher than RCP8.5 if the nitrogen cycle is included. In contrast, nitrogen deposition in RCP2.6 would be much smaller than that in RCP8.5 due to emission controls (Fig. S5), leading to lower nitrogen supply for the ecosystem in the former scenario. Consequently, plant photosynthesis is confronted with a stronger nutrient limit in RCP2.6 than that in RCP8.5, resulting in lower CO₂ fertilization efficiency in the former scenario. The net effect of the nitrogen cycle on the land carbon cycle is very uncertain (Zaehle et al., 2014; Huntzinger et al., 2017; Xiao et al., 2015).

For a warming target of 1.5 ^∘C, our analyses suggest that an associated reduction of CO₂ and pollution emissions brings greater benefits to ecosystems in China than a pathway without emission controls. The slow changes of temperature and other environmental variables due to slow growth of CO₂ are helpful for plant adaptation and limit biome shift (Warszawski et al., 2013), and the lower O₃ and higher solar radiation from aerosol removal increase plant photosynthesis. Consequently, China's ecosystems mitigate 10.6±1.4 % of national emissions in the stabilized pathway, more efficient than the fraction of 6.3±0.8 % in the transient pathway, allowing for a larger carbon budget for economic development.