Absolute Contribution of the Non-Uniform Spatial Distribution of Atmospheric CO₂ to Net Primary Production through CO₂-Radiative Forcing

Abstract

Atmospheric concentrations of CO₂ are the most important driver of the Earth’s climate and ecosystems through CO₂-radiative forcing, fueling the surface temperature and latent heat flux on half-century timescales. We used FGOALS-s2 coupled with AVIM2 to estimate the response of net primary production (NPP) to spatial variations in CO₂ during the time period 1956–2005. We investigated how the induced variations in surface temperature and soil moisture influence NPP and the feedback of the oceans and sea ice on changes in NPP. The spatial variations in the concentrations of CO₂ resulted in a decrease in NPP from 1956 to 2005 when we included ocean and sea ice dynamics, but a slight increase in NPP without ocean and sea ice dynamics. One of the reasons is that the positive feedback of sea temperature to the surface temperature leads to a significant decrease in tropical NPP. Globally, the non-uniform spatial distribution of CO₂ absolutely contributed about 14.3% ± 2.2% to the terrestrial NPP when we included ocean and sea ice dynamics or about 11.5% ± 1.1% without ocean and sea ice dynamics. Our findings suggest that more attention should be paid to the response of NPP to spatial variations in atmospheric CO₂ through CO₂-radiative forcing, particularly at low latitudes, to better constrain the predicted carbon flux under current and future conditions. We also highlight the fundamental importance of changes in soil moisture in determining the pattern, response and magnitude of NPP to the non-uniform spatial distribution of CO₂ under a warming climate.

1. Introduction

Atmospheric CO₂ concentrations are one of the most important drivers of climate in earth system models [1,2]. It can affect the global average surface temperature, the global water cycle, global sea-levels, the loss of Arctic sea ice and the atmospheric vapor pressure deficit through radiative forcing [3,4,5]. Changes in CO₂ concentrations indirectly affect the growth of land plants and the productivity of terrestrial vegetation (gross primary production) [6]—for example, the gross primary production of terrestrial vegetation decreased after the late 1990s as a result of the increased vapor pressure deficit [5].

There is large uncertainty in the magnitude and spatial distribution of the radiative effect of atmospheric CO₂ concentration (±20%) based on multi-model simulations [7,8]. These models come from the fifth stage of the coupled model comparison project (CMIP5) [9] using the uniformity of atmospheric CO₂ concentration without spatial varying. Although the spatial heterogeneity of atmospheric CO₂ concentrations varies widely on a regional scale [10,11], the contribution of spatially non-uniform CO₂ to the Earth’s climate and the terrestrial carbon cycle through CO₂-radiative forcing is largely unknown [12].

Several studies have reported the impacts of varying atmospheric CO₂ concentrations directly through CO₂-radiative forcing on climate state [13,14]. Furthermore, it can also disturb the carbon balance [15] and accumulation of carbon on land [16]. For instance, the variations in surface temperature, precipitation and soil moisture [17] influence the rate of photosynthesis and limit net primary production (NPP) and soil respiration [15]. Thus, changes in the Earth’s climate through CO₂-radiative forcing therefore constrain the processes of the carbon cycle [2].

However, predictions of the rise in temperature under historical conditions remain highly uncertain, ranging from about 0.4 to 1.3 °C, mainly due to historical radiative forcing of greenhouse gases [18]. One of the uncertainties results from the lack of a clear understanding about the feedback of spatial variations in atmospheric CO₂ and the lack of reliable measurements of this process in the field. Modeling approaches are therefore vital to constrain climate variables and the carbon fluxes resulting from the spatial varying of CO₂ concentration.

On the other hand, most conventional methodologies do not consider the influence of spatial varying of CO₂ concentrations on the Earth’s climate system and biogeochemical cycles. We would focus to quantify the response of terrestrial NPP to the spatial varying of CO₂ concentration through CO₂-radiative forcing. To further assess the significance of the non-uniform spatial distribution of CO₂ on NPP, we compared simulations with fully coupled atmosphere–land–ocean–sea ice components (simulation A) and coupled atmosphere–land components (simulation B) in the Flexible Global Ocean–Atmosphere–Land System Model Spectral Version 2 (FGOALS-s2) [19] coupled with the Atmosphere–Vegetation Interaction Model (AVIM) [20]. Specifically, we evaluated NPP in FGOALS-AVIM with the implementation of these two different simulations forced by a uniform or non-uniform spatial distribution of atmospheric CO₂ with or without active ocean and sea ice during the time period 1956–2005. We addressed three principal questions: (1) how much does the climate-driven feedback caused by the non-uniform spatial distribution of CO₂ through CO₂-raiative forcing affect climate variables and NPP; (2) how does the estimated NPP with and without ocean and sea ice dynamics differ as a result of the non-uniform spatial distribution of CO₂ through the radiative forcing associated with spatial variations in the climate variables; and (3) what are the implications of a non-uniform spatial distribution of CO₂ on the global NPP during the time period 1956–2005.

2. Methods and Simulations

The Student’s t-test used in this study is used to estimate of the statistical significance between the variables from different simulations, which are calculated as follows:

$${\mathrm{t}}_{\mathrm{i}}=\frac{ \overline{{\mathrm{x}}_{\mathrm{i}}}-\overline{{\mathrm{y}}_{\mathrm{i}}} }{{\mathrm{S}}_{\mathrm{Pi}}\sqrt{\frac{1}{{\mathrm{n}}_1}+\frac{1}{{\mathrm{n}}_2}}}$$

(1)

where $\overline{{\mathrm{x}}_{\mathrm{i}}}$ or $\overline{{\mathrm{y}}_{\mathrm{i}}}$ is the averaged NPP from simulations with spatial varying of CO₂ concentration or without spatial varying of CO₂ concentration of the i-th grid cell during 1956−2005, ${\mathrm{S}}_{\mathrm{Pi}}$ is the weighted average of the standard deviations of NPP from the two simulations during the same period in the i-th grid cell. ${\mathrm{n}}_1$ or ${\mathrm{n}}_2$ is the equivalent sample sizes of two simulations.

Based on the absolute contribution of each region (each grid cell or PFT) to global NPP, we further classify global NPP into different categories (Figure 1). To this end, we adopted an index to score each geographic location based on consistency. The magnitude of the absolute change in local NPP flux relative to the global NPP under CO₂ uniform distribution was calculated as the following, which was also estimated by [21]:

$${\mathrm{C}}_{\mathrm{j}}=\frac{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{n}}|\Delta {\mathrm{NPP}}_{\mathrm{jt}}|}{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{n}}{\mathrm{NPP}}_{\mathrm{t}}}$$

(2)

where |ΔNPP_jt| is the absolute change in NPP caused by the non-uniform spatial distribution of CO₂ relative to a uniform spatial distribution of CO₂ through CO₂-radiative forcing for the jth region at the tth time (units: g C yr⁻¹) and NPP_t is the NPP for the tth time period using forcing by a uniform spatial distribution of CO₂ for the whole terrestrial ecosystem. C_j is the contribution of a non-uniform spatial distribution CO₂ through CO₂-radiative forcing for the jth region relative to the global NPP using forcing by a uniform spatial distribution of CO₂.

FGOALS that we used here participated in the CMIP6 and was used as an assessment tool in the 6th Assessment Report of the Intergovernmental Panel on Climate Change. It is a fully coupled earth system model and is composed of four separate models simultaneously simulating the Earth’s atmosphere (SAMIL2), oceans (LICOM2), land surface (AVIM) and sea ice (CSIM). FGOALS includes one central coupler component (CPL6). It has an interactive carbon cycle model in the land component and an ecosystem–biogeochemical module in the ocean component. The simulated atmospheric CO₂ concentrations are fully coupled to the land CO₂ fluxes and are used directly to compute radiative forcing. The release of methane from the melting of permafrost potentially has a huge impact on warming, but FGOALS-AVIM currently has only very simple carbon permafrost models and no release of marine methane is included.

We used FGLOAL-AVIM at a resolution of (2.81° × 1.66°) [22]. This version of FGOALS has previously been used to assess the effects of CO₂-radiative forcing on the total cloud fraction, temperature and water vapor under different abruptly quadrupling CO₂ concentrations [23]. Evaluation of the mean summer evapotranspiration and mean annual runoff in the Lake Baikal basin [24], and major global biogeochemical fluxes and pool sizes of carbon [25,26,27,28].

Time series of the atmospheric CO₂ concentrations (i.e., fossil fuel burning, cement manufacturing and gas flaring in oilfields) from CMIP6 are available from 1850 to 2014 at a monthly resolution. The atmospheric CO₂ variable from CMIP6 only considers spatial variations in terms of latitude rather than longitude. We therefore used the Open-Data Inventory for Anthropogenic Carbon dioxide (ODIAC) [29] at a resolution of (1° × 1°) to establish a quantitative relationship between the carbon emissions in each pixel and the latitudinally averaged carbon emissions to fully reflect the spatial heterogeneity of atmospheric CO₂ and to describe the impact of anthropogenic carbon emissions on its spatial distribution in the atmosphere. We assumed that this quantitative relationship is also applicable to the relationship between the atmospheric CO₂ variable of each pixel and the latitudinally averaged atmospheric CO₂ from CMIP6. The CO₂ variable was produced by fully considering the spatial heterogeneity, which reflected the effect of anthropogenic carbon emissions. The CO₂ variables for 1850–2005 were generated by re-gridding from a spatial resolution of (1° × 1°) to (2.81° × 1.66°). The CO₂ variable was then used as input data for our model.

To quantify the effects of a non-uniform spatial distribution of CO₂ on the carbon flux through CO₂-radiative forcing, we carried out two different simulations using FGOALS-AVIM with one of the two different components used to estimate NPP. We ran each experiment for a simulated 156 years from 1850 to 2005 (

Table 1. Forcing data for simulations.

Simulation	CO2 Forcing	Component
A1	Spatial and temporal variations	Fully coupled atmosphere–land–ocean–sea ice
B1	Only temporal variations	Fully coupled atmosphere–land–ocean–sea ice
A2	Spatial and temporal variations	Only atmosphere–land coupling
B2	Only temporal variations	Only atmosphere–land coupling

). Simulation A1 includes the non-uniform space–time-varying atmospheric concentrations of CO₂ based on the active atmosphere–land–ocean–sea ice components. Simulation A2 fixes only one atmospheric CO₂ level at the global scale with only the time-varying atmospheric concentrations of CO₂ without the spatial variations described by [30] and fully considers the interactions of the atmosphere–land–ocean–sea ice components. Simulation B1 uses the non-uniform time–space-varying atmospheric concentrations of CO₂ and the coupled atmosphere–land interactions without ocean and sea ice dynamics. Simulation B2 uses the uniform time-varying atmospheric concentrations of CO₂ and the coupled atmosphere–land interactions without ocean and sea ice dynamics.

The difference in the land NPP between simulations A1 and A2 or simulations B1 and B2 represents the effect of the non-uniform spatial distribution of atmospheric concentrations of CO₂. It is well known that physiological forcing is a robust driver influencing terrestrial plants and the accumulation of carbon on land [17,31], but we did not consider physiological forcing in the land component, in which atmospheric concentrations of CO₂ were fixed at the current level of ~355 ppmv.

Simulated NPP from TRENDY used in this study can be freely downloaded from TRENDY|Dynamic Global Vegetation Model Projects (ceh.ac.uk) (see

Table 2. Datasets used in this study.

Variable	Project	Period	Website	Ref.
NPP	TRENDY	1956–2005	TRENDY|Dynamic Global Vegetation Model Projects (ceh.ac.uk)	[32]
CO2 emission	ODIAC	1980–2014	https://db.cger.nies.go.jp/dataset/ODIAC/DL_odiac_v1.7.html (accessed on 19 August 2019).	[29]

). The selected models mainly include CLM4C, CLM4CN, LPJ, LPJGUESS, OCN, SDGVM, TRIFFID and CABLE. The natural response of NPP across the terrestrial ecosystem is one of our focuses. Thus, S2 experiment from TRENDY was chosen, which did not include the impact of land use. For more detailed information about TRENDY simulation protocol could be seen in [32,33]. Previous studies have proved that results of TRENDY perform well and are also widely used to estimate the global carbon budget [6]. Compared with the observed data, TRENDY reproduced variations in carbon fluxes of the terrestrial ecosystem to a certain extent [34,35].

3. Results

Under historical conditions from 1956 to 2005, NPP from simulations A1, A2, B1, and B2 was estimated. The spatial distribution of PFT is shown in Figure 1. In the following results, at the PFT level, the changes in NPP caused by the non-uniform distribution of CO₂ would be given.

3.1. Estimates of CO₂ Concentrations and NPP

There was a large spatial heterogeneity in atmospheric concentrations of CO₂ between the northern and southern hemispheres in the time period 1956–2005. The non-uniform spatial distribution of CO₂ was lower in Australia, most of Africa and southern South America, but higher in southern Europe, the eastern USA and southeast Asia. During the same time period, simulations A1 and B1 were driven by the same CO₂-radiative forcing and simulations A2 and B2 were enforced by the same CO₂ data. The spatial variations in CO₂ ranged between 339 and 351 ppmv (Figure 2a). Greater differences between simulations A1 and A2 or B1 and B2 were mainly located in southern Europe, the eastern USA and southeast Asia (Figure 2b).

The NPP was estimated for simulations A1, A2, B1 and B2 during the time period 2000–2005 under current conditions. This period was chosen so that our results were comparable with NPP data from the TRENDY datasets. The TRENDY dataset consists of an ensemble of seven process-based terrestrial ecosystem model simulations [36]. We selected experiment S2 without land use changes and fires and including CLM4C, CLM4CN, LPJ, LPJ-GUESS, OCN, SDGVM and TRIFFID [37] from 1901 to 2010.

Our estimated land NPP for the four simulations of 51.39–56.47 Pg C yr⁻¹ was lower than the estimates of about 70 Pg C yr⁻¹ of [38], but higher than the 41 Pg C yr⁻¹ estimated by [39]. The estimated NPPs of the four simulations were closer to the observed estimate of about 56.02 Pg C yr⁻¹ of [40]. Globally, the mean annual NPP in the historical period from CMIP5 was about 62.6 Pg C yr⁻¹ [41], whereas it was 55.53 ± 1.69 Pg C yr⁻¹ from the TRENDY dataset. The mean annual NPP in the historical period simulated by FGOALS-AVIM was 53.80 ± 2.54, 53.75 ± 2.95, 58.75 ± 1.46 and 59.36 ± 2.73 Pg C yr⁻¹ from simulations A1, A2, B1 and B2, respectively (Figure 3).

3.2. Estimates of Climatic Variables under Present Conditions

Over the simulation period 1956–2005, the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing resulted in a global increase in low cloudiness of about 0.002% ± 0.005% (Figure 4a), a net solar flux of 0.52 ± 0.55 W m⁻² including the atmosphere–land–ocean–sea ice dynamics (Figure 4e) and a terrestrial latent heat flux of 0.36 ± 0.97 W m⁻² at the global scale (Figure 4c). At the global scale, the climate variables varied considerably among the different components in response to the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing. During the time period 1956–2005, the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing based on only atmosphere and land interactions caused an increase in the net solar flux of about 0.45 ± 0.31 W m⁻² and a decrease in global low cloudiness of 0.001% ± 0.003%, especially over the ocean (Figure 4b;

Table 3. Response of variables to CO₂-radiative forcing as a result of the non-uniform spatial distribution of CO₂.

Variable	Range	A1–A2	B1–B2
CO2	Global	−0.0007 ± 0.0037
Precipitation (mm yr−1)	Land	6.77 ± 21.45	−9.26 ± 14.72
	Ocean	5.89 ± 14.15	−9.58 ± 9.92
Surface temperature (°C)	Land	0.57 ± 0.45	0.02 ± 0.05
	Ocean	0.45 ± 0.38	0.07 ± 0.08
Latent heat flux (W m−2)	Land	0.59 ± 1.23	−0.64 ± 0.93
	Ocean	0.38 ± 0.80
Sensible heat flux (W m−2)	Land	−0.29 ± 0.46
	Ocean	−0.20 ± 0.42
Net radiation flux at Earth’s surface (W m−2)	Land	0.79 ± 0.69
	Ocean	0.59 ± 0.63
Soil evaporation (kg m−2 yr−1)	Land	3.99 ± 7.52
Vegetation evaporation (kg m−2 yr−1)	Land	0.43 ± 4.90
Soil transpiration (kg m−2 yr−1)	Land	1.91 ± 5.70
Surface runoff	Land	3.00 ± 6.18
Soil moisture (kg m−2)	Land	0.03 ± 0.04

). By contrast, the simulated global latent heat flux based on atmosphere and land interactions decreased by about −0.64 ± 0.93 W m⁻² (Figure 4d).

To further analyze the causes of the differences in the estimated NPP resulting from the non-uniform spatial distribution of CO₂ through CO₂-radiatvie forcing, we estimated the changes in the vital variables (e.g., surface temperature, precipitation, evapotranspiration and soil moisture in the different components) for the historical time period 1956–2005. Non-uniform CO₂ from simulations fully coupled atmosphere–land–ocean–sea ice components caused an increase in surface temperature ranging from 0.57 ± 0.45 °C in the terrestrial ecosystem, where it from simulations only coupled atmosphere and land without active ocean and sea ice components led to less increased surface temperature of 0.02 ± 0.05 °C. The changes in the surface temperature at the regional level correspond well to the changes in the non-uniform spatial distribution of CO₂ considering only the atmosphere–land interactions (Figure 5b). The responses of the surface temperature in tropical regions to the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing based on the fully coupled atmosphere–land–ocean–sea ice components were up to >0.2 °C, especially in the Amazon basin, central Africa and tropical Asia. By contrast, the simulations with only the atmosphere–land interactions simulated a small increase or even decrease in these regions.

In the time period 1956–2005, precipitation over the land surface increased in most of the northern hemisphere, but not in southern Europe, the eastern USA and southeast Asia. This was a response to the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing (Figure 5c) compared with the uniform spatial distribution of CO₂ in the same period. The simulated increase in precipitation influenced by non-uniform spatial distribution of CO₂ considering the full atmosphere–land–ocean–sea ice interactions was much greater than that which considered only atmosphere–land interactions, apart from in the Amazon basin and tropical Asia. This was a result of the larger increase in the net radiation flux from surface warming and latent heat flux over land, especially in the northern hemisphere.

Correspondingly, the responses of the simulated change in evapotranspiration by the two differently coupled components differed considerably in both the direction and magnitude of change for most terrestrial ecosystems (Figure 5d,e). The two different components predicted an increase in evapotranspiration affected by the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing relative to the uniform spatial distribution of CO₂ for most land surfaces. The simulated response of evapotranspiration to changes in the non-uniform spatial distribution of CO₂ differed at the regional scale. For example, evapotranspiration north of 50° N in the northern hemisphere was predicted to increase as a response to the uniform spatial distribution of CO₂ through CO₂-radiative forcing based on the fully coupled atmosphere–land–ocean–sea ice model, but to decrease in southern Europe, the eastern USA and southeast Asia. By contrast, evapotranspiration influenced by the non-uniform spatial distribution of CO₂ considering only the atmosphere–land interactions showed a much greater increase than the fully coupled atmosphere–land–ocean–sea ice interactions at low latitudes and in the eastern USA.

3.3. Responses of NPP to Spatial Varying of CO₂ Concentrations through Radiative Forcing

To quantify the impact of the non-uniform spatial distribution of CO₂ on NPP, we simulated the NPP for the time period 1956–2005 with or without the uniform spatial distribution of CO₂ after 1850. Specifically, we ran FGOALS-AVIM with CO₂ held at the same value at the global level for each land grid point in simulation A2 or B2 and compared the results with the non-uniform spatial distribution of CO₂ for each land grid point in simulation A1 or B1 for the time period 1956–2005.

We estimated that the decrease in NPP caused by the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing during the time period 1956–2005 was −3.39 ± 12.33 g C m⁻² yr⁻¹ or about −0.5 ± 1.88 Pg C yr⁻¹ for the fully coupled atmosphere–land–ocean–sea ice components (simulation A1 minus A2). The increase was about 0.88 ± 8.86 g C m⁻² yr⁻¹ or 0.1 ± 1.35 Pg C yr⁻¹ when considering only atmosphere–land interactions (Figure 6a,b).

The fully coupled atmosphere–land–ocean–sea ice components estimated a decrease in NPP caused by the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing for terrestrial ecosystems (Figure 6a), whereas the coupled active atmosphere–land components predicted an increase (Figure 6b). The mechanism contributing to the estimated change in NPP was similar between the two simulations. With the simulated decrease in soil moisture (Figure 6c), the model based on the fully coupled atmosphere–land–ocean–sea ice components estimated a decrease in NPP for southern Europe, the eastern USA, southeast Asia, the Amazon basin and central Africa (Figure 6a). However, the simulation considering only the atmosphere–land interaction estimated an increase in soil moisture (Figure 4d) for the eastern USA, the Amazon basin and central Africa (Figure 6d). The estimated increase in NPP for these regions using only the active atmosphere–land components was attributed to the relatively higher increase in NPP (Figure 6b).

We used the present vegetation cover for the period of the 1990s based on the plant functional types (PFTs) of the International Geosphere Biosphere Program (IGBP) data [42]. At the PFT-level, the NPP simulated by simulations A due to spatial varying in CO₂ concentrations were reduced for EBF, DBF, C4 and crop, but increased for tundra. The variations in NPP simulated from simulations A and B were feeble for ENF. Globally, the changes of annual NPP from 1956 to 2005 resulting from varying in CO₂ concentrations was ~−0.5 Pg C yr⁻¹ from simulations A or estimates of ~0.1 Pg C yr⁻¹ from simulations B (Supplementary Materials Figure S1).

The non-uniformity of CO₂ leads to substantial changes in surface temperature and water fluxes, and the changes caused by non-uniformity of CO₂ in NPP depend on changes in surface temperature and soil moisture, rather than the impact of CO₂, which may lead to changes in the spatial distribution of NPP (Figure 7). Figure 7b,e show how much change in NPP is locally related to changes in surface temperature, especially in the tropics and the Southern Hemisphere. Obviously, these two variables are significantly correlated, which indicates that changes in surface temperature can explain the changes in NPP in these regions. For example, considering that the fully coupled land-atmosphere-ocean-ice interaction, rising surface temperature leads to negative changes of NPP in tropics due to non-uniform CO₂. On the contrary, in most terrestrial ecosystems, there is a significantly positive correlation between soil moisture and local NPP (p < 0.05) (Figure 7c,f). Therefore, our results indicate that the competition between the impacts of surface temperature and soil moisture on NPP is one of the mechanisms leading to changes in NPP in the tropics.

Although the absolute change in CO₂ at the global scale was estimated to be only about 0.3%, the model based on coupled components A1 minus A2 or B1 minus B2 estimated an additional change in NPP caused by the non-uniform spatial distribution of CO₂ by 14.26% or 11.51% during the time period 1956–2005 (Figure 8). This reflects that the impact of CO₂ heterogeneity on carbon cycle was amplified by the climate effect. In addition, the difference between the two different coupled components probably reflects the contribution of the ocean and sea ice to the climate system. With global warming caused by the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing, the model based on the fully coupled components A1 minus A2 estimated absolute variations in NPP caused by the non-uniform spatial distribution of CO₂ relative to the global NPP driven by a uniform spatial distribution of CO₂ ranging from 0.97 to 1.53% for C3 and from 2.92 to 4.03% for crops. The simulation with only atmosphere–land interactions ranged from 1.76 to 2.23% for C3 and from 2.56 to 3.32% for crops. For different PFTs (e.g., crops), the absolute changes in CO₂ from 0.05 to 0.10% were estimated to be higher than those for other PFTs. The impact of the non-uniform spatial distribution of CO₂ was largest for deciduous broadleaf forest, for which simulation A1 minus A2 estimated higher absolute changes in NPP from 2.53 to 4.65% for evergreen broadleaf forest due to the non-uniform spatial distribution of CO₂, compared with simulation B1 minus B2 from 1.93 to 3.05%. For crops, both simulations A1 minus A2 and B1 minus B2 estimated high absolute changes ranging from 2.56 to 3.32%.

4. Discussion

At present, few studies have considered effect of radiative forcing of spatial varying of CO₂ concentrations in the atmosphere on NPP. Our results showed that the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing has important implications for carbon fluxes. We estimated that the non-uniform spatial distribution of CO₂ contributed to a decrease in NPP, although most other studies have not considered this. The CMIP5 provides an ensemble mean annual NPP of 62.6 Pg C yr⁻¹ over the time period 1995–2004 [41]. However, this analysis did not consider the feedbacks from the non-uniform spatial distribution of CO₂ to the Earth’s climate system. Our analysis showed an effect of the non-uniform spatial distribution of CO₂ on NPP without considering land use of up to 14.26 ± 2.22% of the NPP simulated by a fully coupled model. This highlights the tremendous leverage of a non-uniform spatial distribution of CO₂ on the global carbon cycle and we recommend a more substantial focus on understanding this process in the biosphere through the influence of the climate system based on coupled models that include dynamical sea ice and ocean components.

In the past 30 years, there has been an increase in NPP at the global scale, which is likely driven by CO₂ fertilization [34,35]. However, the increase in CO₂ radiative forcing limits this positive effect [1]. The exact impact of spatial varying in CO₂ concentration on the carbon cycle through radiative forcing is not yet fully understood. Our study illuminates that from 1956 to 2005, the spatial varying in the CO₂ concentration resulted in a reduction of NPP, which mainly occurred in the tropics and the middle-latitudes of the Northern Hemisphere. Nevertheless, as shown in Figure 7, a relationship between variations in CO₂ concentration and NPP could not pass a significant level. It might result due to not considering the impact of CO₂ fertilization. On the contrary, spatial varying in CO₂ concentration through the radiation effect influences temperature, precipitation and soil moisture, which ultimately affects NPP. The northern hemisphere accounts for nearly 68% of the global area, which has a higher CO₂ concentration than the globally averaged level. Higher temperature has a negative impact on NPP in tropics and mid-latitudes of the Northern Hemisphere [33,38,43,44]. In addition, considering ocean feedback, an increase in temperature could reduce NPP in these regions. Conversely, when the ocean feedback is shielded, increases in temperature is slight in tropics. At the same time, in this condition, the increase in the soil moisture would offset the negative impact of warming in tropics.

To estimate the contribution of the ocean and sea ice influenced by the spatial varying of CO₂ concentrations through CO₂-radiative forcing to the terrestrial NPP, we carried out simulations B1 and B2, in which the effect of CO₂-radiative forcing was assessed using fixed the sea surface temperature (SST) in 1980. The additional increase in SST was estimated to be 0.45 ± 0.38 °C using the fully coupled model, compared with 0.07 ± 0.08 °C based on only the coupled active atmosphere–land components. This can be understood by the positive feedback from a greater increase in the SST (Figure 9a) to surface temperature from simulations A1 minus A2 relative to from simulations B1minus B2 without varying SST. The response of NPP in low-latitudes of terrestrial ecosystem to changes in surface temperature might be asymmetric. This may also lead to asymmetric NPP changes related to SST changes. For example, in some tropical terrestrial regions, increased SST has a negative impact on NPP under warmer changes [45]. Therefore, more detailed quantification is needed to understand the asymmetric response of NPP to SST changes.

The spatial response of soil moisture differed widely, although the SST had a key influence soil moisture [46]. Under warmer conditions, a non-uniform spatial distribution of CO₂ decreased soil moisture across the eastern USA, southern Europe and southeast Asia in the fully coupled model. This is consistent with previous research [47], which showed that an increased SST in the tropics would increase the depletion of soil moisture in the eastern USA, southern Europe, southeast Asia and the tropics.

At high latitudes (e.g., the Arctic), the decrease in the area of ice as a result of the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing is >2% (Figure 9b). An additional increase in precipitation has been shown in most of this region. This change can be explained by emergent positive feedbacks between the increasing area of melted ice and the increased water vapor at high latitudes in the fully coupled model. Recent studies have found that a decrease in Arctic sea ice has an important role in the Earth’s climate and ecosystems [48]. Specifically, stronger warming results in an increased loss of sea ice, which, in turn, simulates more melting through a reduced albedo [48]. Ice-loss feedback could amplify changes in the Earth’s climate at high latitudes [48,49]. We estimated larger increases in surface temperature and precipitation for tundra regions using the fully coupled model including sea ice dynamics than the model without sea ice dynamics. There is a positive feedback from an increase in melted sea ice in a warmer and wetter climate, which enhances the NPP by the non-uniform spatial distribution of CO₂ through CO₂-radiative forcing in the tundra regions.

This finding is consistent with a previous study of CO₂-radiative forcing and has important implications for our understanding of the changes in carbon flux associated with climate change [50]. Changes in the non-uniform spatial distribution of CO₂ concentrations influenced changes in the heat budget (Figure 4), then continued to regulate the spatial distribution of surface climate variables. Changes in the atmospheric circulation (Figure 10) introduced by the non-uniform spatial distribution of CO₂ also could influence the water balance and, in turn, affect soil moisture. Therefore, the patterns and responses of the water budget at the regional scales triggered by the non-uniform spatial distribution of CO₂ determine the patterns and responses in NPP. In particular, near-global increases in specific humidity do not mean that this increase necessarily results in an increase in precipitation. The simulation using the fully coupled model estimated that weaker westerlies reduced the input of humid air from the Atlantic Ocean to southern Europe, which reduced precipitation and then increased the depletion of soil moisture. The eastward radial circulation in the eastern USA was strengthened. Such a circulation pattern does not favor precipitation. The eastward radial circulation was also enhanced in East Asia and suppressed the monsoon from the Eastern Pacific Ocean, limiting precipitation and inducing droughts.

CO₂-physiological forcing could also have significant impact on the spatial variation of NPP [51]. However, it was not included in our results, because our main focus is on the impacts of CO₂ radiative forcing on NPP. For terrestrial processes, atmospheric CO₂ was uniformly applied to all PFTs at a fixed level of 355 ppmv. The impact of a spatial variation in CO₂ concentrations was masked in our analysis. Results may be different if we include a variation in the spatiotemporal distribution of atmospheric CO₂ through CO₂-physiological forcing in the land component. This needs further investigation.

5. Conclusions

This study investigated the regulations of spatial varying in atmospheric CO₂ concentration through radiative forcing to annual NPP. The results implied that soil moisture resulting from spatial varying in CO₂ concentration dominates the changes of annual NPP, and the increased SST substantially reduced NPP in tropics. As a consequence, spatial varying of atmospheric CO₂ concentration did change the spatial pattern of NPP under the historical conditions from 1956 to 2005. In addition, our results showed the current offline terrestrial models only considering land and atmosphere interactions did not reproduce the feedback of ocean and sea ice to changes in NPP. In general, this study represented the key role of spatial varying in CO₂ concentration in regulating variations in NPP.