Global Vegetation-Temperature Sensitivity and Its Driving Forces in the 21st Century

Abstract

It has been projected that climatic warming will contribute to vegetation productivity variability at the global scale. With a continued warming, to what extent and where the vegetation productivity is most affected by warming has still not been adequately quantified. Herein, based on 11 earth system model outputs, we predict the characteristics of vegetation-temperature sensitivity (S_vpt, defined as higher/lower temperature produce more/less vegetation) changes under different CO₂ emission scenarios and various vegetation types, further assessing the relationship of the S_vpt to socio-ecosystems. At the end of the 21st century, the area proportion with the global temperature increases >2°C under the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios are 3%, 40%, 97%, and 99%, respectively. The largest contribution to the global terrestrial gross primary productivity growth is at low latitudes. 33%–63% of global terrestrial ecosystems show a significantly negative trend in the S_vpt, indicating the weakened promotion effect of warming on vegetation growth. In particular, in 2060, there will be a clear reversal of the trend from carbon sink to carbon source under the SSP3-7.0 scenario, spatially distributed mainly in the Amazon rainforest, tropical Africa, and southern North America. Precipitation is also an important factor affecting the S_vpt change, and there is an inverted U-shaped relationship between them. Precipitation thresholds are higher under the high emission scenario when S_vpt reaches its highest value. Moreover, socio-demographic pressures in places like Central Africa and East Africa will offset the promotion of vegetation growth by warming; in the future, these countries should develop appropriate population and land management strategies to achieve socio-ecosystems sustainable development.

Introduction

As an important part of the surface ecosystem, vegetation plays a key role in regulating climate change, storing carbon sinks, and maintaining the surface energy balance by absorbing CO₂ from the atmosphere (Ballantyne et al., 2017; Forzieri et al., 2017). Previous studies have shown that increased CO₂ will aggravate global warming to a certain extent, and global warming can enhance vegetation productivity by lengthening the active growing season and improving the maximum photosynthetic rate (Bastos et al., 2019; Nemani et al., 2003; Thomas et al., 2016), that is, warming plays a positive role in improving vegetation productivity (Song et al., 2022). However, recent studies have found that in comparison to the pre-industrialization level, the temperature in the Amazon region in the dry season increased by nearly 3°C, which is three times the global annual mean level, and warming has accelerated the transformation of the Amazon tropical forest from a “carbon sink”—which absorbs CO₂—to a “carbon source”—which emits CO₂ (Gatti et al., 2021); the latter has a negative effect on vegetation growth and increases tree mortality (Green et al., 2020). If urgent measures are not taken to protect the tropical rainforest, its carbon sink will be completely transformed into a carbon source by mid-2030 (Hubau et al., 2020). Furthermore, some studies on tree ring records, satellite observations, and simulated vegetation agents have also shown that the vegetation productivity increase caused by warming in the northern hemisphere may be reversed, and will eventually inhibit the growth trend of vegetation primary productivity (Koch, 2022; Zhang et al., 2022). In other words, global warming caused by increased CO₂ emissions changes the function and balance of vegetation ecosystems, and the vegetation-temperature sensitivity (S_vpt, defined as higher/lower temperature produce more/less vegetation) is a major factor affecting vegetation functioning (Zhang et al., 2022). Thus, a quantification of the vegetation-temperature relationship has been established as a successful proxy for the S_vpt, which can be used as a simplified indicator for vegetation functioning and is applied here accordingly.

The S_vpt has been widely researched over the past few decades. For example, Piao et al. (2014) found a weakening relationship between interannual temperature variability and northern vegetation activity between 1982 and 2011. Shestakova et al. (2019) indicated that the sensitivity of vegetation productivity to summer temperature decreased over the entire 20th century, which exhibited obvious spatial heterogeneity. Furthermore, changes in the S_vpt can also be affected by other climatic forces, such as droughts. Increased drought stress can modify the response of plant growth to temperature changes, and could be a potential cause for the declined temperature sensitivity of tree growth (Angert et al., 2005; D'Arrigo et al., 2004); increases in the vegetation water use efficiency can mitigate the impact of drought on the S_vpt (Field et al., 1995; Peñuelas et al., 2011). These studies have determined how and where the impact of the S_vpt has changed over the past three decades. However, under different CO₂ emission scenarios at the end of the 21st century, the change trend of global S_vpt and its socio-ecological driving forces remains uncertain.

Here, we use empirical remotely sensed data, including the gross primary productivity (GPP), the carbon flux entering plants via photosynthesis, and the mean annual temperature to represent vegetation productivity and atmospheric temperature. To ensure the accuracy of GPP and temperature data, we select 11 earth system models (ESMs) and obtain the mean value of GPP and temperature corresponding to these models, further revising them (see Section 2), and then investigate their spatio-temporal change trends under different scenarios in the future (2015–2100) using the Coupled Model Intercomparison Project Phase 6 (CMIP6). Based on this, we couple future GPP to temperature projections to assess how vegetation productivity changes with the temperature under different CO₂ emission scenarios and vegetation types through sequential linear regression slopes (SeRGS, see Methods) by the end of this century, quantified here using the S_vpt. We then relate the observed changes in the S_vpt to both social and ecological stressors, estimating the relative importance of these potential stressors and further assessing their interconnectivity (Figure 1).

Figure 1. Overview of the framework in this study.

Material and Methods Data Acquisition and Pre-Processing Earth System Model Simulations

GPP is the gross carbon fixed by plants through photosynthesis, which is used as an important metric to represent plant activity and function. Under climate change conditions, to conduct a projection analysis on how ecological–climatic interactions affect terrestrial carbon cycles, the historical temperature, future temperature, GPP, and precipitation were extracted from the CMIP6 outputs. CMIP6 contains a series of ESMs, and each ESM provides combination scenarios that link the socioeconomic and technological development, named the Shared Socioeconomic Pathways (SSPs), with the future greenhouse gas representative concentration pathway (RCPs) (Kriegler et al., 2014; O'Neill et al., 2014). Quantitative elements that constitute SSPs is population and gross domestic product, while the qualitative elements include seven aspects: environment and natural resources, policies and institutions, and technological development, etc. For SSPs, it can be divided into five types that span potential futures of green or fossil-fueled growth (SSP1 van Vuuren et al., 2017, and SSP5 Kriegler et al., 2017), high inequality between or within countries (SSP3 Fujimori et al., 2017, and SSP4 Calvin et al., 2017), and a “middle-of-the-road” scenario (SSP2 Fricko et al., 2017). For RCPs, it contains seven radiative forcing levels of 1.9 W/m², 2.6 W/m², 3.4 W/m², 4.5 W/m², 6.0 W/m², 7.0 W/m² and 8.5 W/m². Based on the socio-economic conditions and radiative forcing in different scenarios, the global climate, vegetation, and other related data has been predicted in the ESMs (Gidden et al., 2019). In this study, the final 11 ESMs were selected in this study (Table S1 in Supporting Information S1) considered four combination scenarios, including SSP1-2.6 (low social vulnerability and radiative forcing, low-emission scenario), SSP2-4.5 (medium social vulnerability and radiative forcing, medium-emission scenario), SSP3-7.0 (medium-to high social vulnerability and radiative forcing, medium-to high emission scenario), and SSP5-8.5 (high social vulnerability and radiative forcing, high-emission scenario). To decrease the differences in the simulation data among different models, only monthly data from the r1i1p1f1 series were used for all models, which are the only collections available from 1982 to 2100. Initially, the selected ESMs have different spatial resolutions; after bias correction and spatial disaggregation (BCSD, see Methods), all of these models provide a resolution of 0.5° × 0.5° in longitude and latitude. The raw monthly data of the three variables were converted into the mean annual temperature, annual GPP, and annual total precipitation. Since the simulation results of GPP (Figure S1 in Supporting Information S1) and temperature (Figure S2 in Supporting Information S1) from 11 ESMs varied, this study adopted the multi-mode ensemble averaging method to calculate the arithmetic mean of all model simulation values to reduce the uncertainty of multi-mode simulations.

Observation Data

The primary monthly temperature and precipitation climate records set was obtained from the Climatic Research Unit at the University of East Anglia (CRU TS v4.05). Based on an archive of monthly mean temperature provided by more than 4,000 weather stations distributed around the world (Harris et al., 2021), these data were produced using angular-distance weighting interpolation and formed of land-only gridding at a spatial resolution of 0.5° from January 1901 to December 2020.

We used FLUXCOM GPP as the observation GPP data, produced through upscaling of local eddy covariance carbon flux tower measurements from 224 globally distributed sites (Tramontana et al., 2016). The predictor variables were based exclusively on satellite-based remote sensing data (Sun et al., 2018). Since we considered vegetation changes under different greenhouse gas emission scenarios in the future, the interannual variability and trend patterns of the selected FLUXCOM GPP product were derived from time-varying meteorological input variables exclusively, only the seasonal cycles of plant growth were constrained by satellite vegetation data (Gampe et al., 2021), and, thus, excluded the associated effects of CO₂ fertilization.

Additional Data

In this study, vegetation types were used to quantify the S_vpt values under different ecosystem types, which were obtained through MCD12C1 global land cover products originating from annual MODIS Terra and Aqua data. After processing the raw data, 17 global land cover types are included (referring to global vegetation classification scheme from the international geosphere-biosphere program), as shown in Figure S3 in Supporting Information S1.

To further understand the relationship between S_vpt and the socio-ecosystem nexus, we divided the world into different population density and income level regions. Population density classification at the spatial scale is based on the fourth version of the Gridded Population of the World (Figure S4 in Supporting Information S1), as lower, low, middle, and high population density regions, and the classification criteria are shown in Table S2 in Supporting Information S1. In addition, gross national income (GNI) is widely used to measure income levels and the quality of life across countries; based on the GNI in US$ from 2015, as provided by the World Bank (Figure S5 in Supporting Information S1), we divided the income groups into four classes: low, lower-middle, upper-middle, and high income groups, and the criteria are shown in Table S3 in Supporting Information S1.

Methods BSCD for CMIP6 Outputs

Before we provide details, we give an overview of our approach as follows. We use ESMs to project future changes, and extract diagnostics of temperature, precipitation and GPP from them, and use contemporary measurements of these three variables to remove ESM biases. We also create partial correlations between GPP and temperature and precipitation, and then use these correlation values to explain their relationships on the pixels. We applied a sequential linear regression slope as explanatory statistics to determine the potential impacts of an altered climate change on vegetation productivity.

CMIP6 outputs cannot be used directly to analyze climate changes due to their coarse spatial resolution and substantial biases (Zhang, Shen, et al., 2019). BCSD has been shown to be an effective way to reduce inherent errors and increase the resolution of the outputs of climate models (Thrasher et al., 2012); it includes two relatively independent parts, namely bias correction and spatial downscaling. The process of BCSD refers to studies conducted by Wood et al. (2004), Nahar et al. (2017), and Xu et al. (2021).

Bias correction adopts the area-weighted average method to interpolate the high-resolution observations to the grids of the CMIP6 model being processed. For this study, observation-based data from 1982 to 2014 was used as the basis to correct the cumulative distribution function of ESM data on each grid. We selected the cubic smoothing spline to build the model for the transformation function to ensure that the observations have the same probability distribution as ESMs. Considering seasonal variations in the biases and the spatial patterns of precipitation and temperature, all calculations were performed separately for each month of the year (396 values, 12 values for each year). After all ESM values for a particular variable were bias-corrected, a spatial disaggregation was performed. The observation-based climate data was aggregated to the scale of the ESMs using an area-weighted average method. Then, the correction factor was interpolated to the original high-resolution observation-based climate data through bilinear interpolation (the correction factor is the ratio of the ESM data after bias correction to the monthly average observations in ascending order). The final product of spatial downscaling was obtained by multiplying the interpolated correction factor field by the high-resolution observed observations. We compared the cumulative distribution function of the ESMs before and after BSCD, and the analysis reveals that the distribution of the corrected data is closer to the true values (Figure S6 in Supporting Information S1).

The Sensitivity of Vegetation Productivity to Temperature (S_vpt) Calculation

We calculated the correlations between GPP and temperature (R_GPP-Temp) for the period of 1982–2014 through Pearson's partial correlation analysis under different scenarios (Sha et al., 2022). In the partial correlation analysis, we controlled for the precipitation variable to ensure significant spatial correlation between vegetation and temperature for each grid used in the calculation.

The sensitivity of vegetation productivity to temperature (S_vpt) defined as higher/lower temperature produce more/less vegetation. Here, we introduced a method called SeRGS as a proxy of S_vpt. SeRGS is based on an ordinary least squares regression with the assumption of a direct correlation between vegetation growth and climate dynamics, regardless of the time-lagged and cumulative effects (Abel et al., 2019). The advantage of this method is that by moving sequentially along the time axis of the study area time series, the changes in vegetation productivity caused by per unit temperature changes are understood better, thereby better exposing the signal of the underlying causal processes (Figure 2). The combination over a spatial moving window of 7 × 7 pixels and a 4-year temporal moving window was chosen in this study to generate the spatial and temporal dynamics of the S_vpt. This combination of temporal windows ensured a high number of effective pixels and the maximum significance of the regression-fitted results (Abel et al., 2021). A linear regression between vegetation (using GPP as a proxy) and the temperature was fitted within a given pixel window, and the regression slope was assigned to the central pixel of the spatially moving window and the central year of the temporal window, respectively. By sequentially moving the temporal window in 1-year steps along the entire time series, a time series of slopes was eventually generated, referred to as the S_vpt time series in this study. To exclude abnormal values, the regression was not fitted and no data were assigned to the central pixel when more than three-quarters of the pixels were of low quality.

Figure 2. Example of sequential linear regression sequential linear regression slopes (SeRGS) calculation for one spatial moving window. 3-dimensional input data (a) for the combined spatial window (7 × 7 pixels) at central (x, y) and temporal window (4 years) at central position t. Linear regression between gross primary productivity and temperature (b) and assignment of the regression slope S (xy, t)to the central pixel of the spatial window (c). Sequentially repeating b and c results in a time series of SeRGS.

Change Trends and Tests on the S_vpt

The Theil–Sen slope estimator was used to determine the trend of the SeRGS time series and indicate the directions of change in the S_vpt. In general, positive trends suggest that vegetation becomes more responsive to the available temperature over time. Negative trends indicate lower vegetation productivity per unit temperature, indicating that vegetation might be inhibited in its function and, thus, converted. The Mann–Kendall (MK) test was applied to assess the sensitivity trends at the 95% significance level (P < 0.05). This approach is known as the distribution-free test, which eliminates the need for the sample to be followed with a certain pattern without the interference of a few outliers, and has been widely used in the trend analysis of various environmental indicators (Mann, 1945; Montibeller et al., 2022; Zhang, Shen, et al., 2019). However, with over 50,000 raster points involved in the region we studied, the significance level of just one hypothesis test may not be sufficient to ensure that the overall error rate is kept within acceptable limits. Therefore, to ensure that we can draw more reliable and accurate conclusions, and to reduce the risk of incorrectly rejecting hypotheses, we adopted the Benjamini-Hochberg method. This method controls for false discovery rate (FDR) by sorting and adjusting p-values and identifying hypotheses that are significant under FDR control. Figures S7 and S8 in Supporting Information S1 show the operation of the FDR procedure for the ACCESS-ESM1-5 model for the historical period and the four scenarios, respectively (see Duque et al., 2019; Wilks, 2016 for more discussion and details). In subsequent analyses, we used the adjusted p-values by Benjamini-Hochberg method to represent the significance of the S_vpt trend.

The observation-based GPP and temperature data from 1982 to 2014 were extracted to assess the observation-based S_vpt, and the results were compared with the ESMs-based S_vpt for this period to verify the reliability of the selected ESMs. The results of the ESMs were found to be in substantial accordance with the spatial pattern of the observation-based S_vpt (Figure S9 in Supporting Information S1), along with a similar distributional pattern throughout the entire value range (Figure S10 in Supporting Information S1), indicating that the all 11 ESMs selected for this study are reliable (correlation coefficients are all greater than 0.4, and the proportion of grids with the same trend is more than 50%). Furthermore, 14% of the grids have the same trend among the 11 ESMs (Figure S11 in Supporting Information S1). We retain the 11 ESMs for our subsequent analysis.

Results Vegetation and Temperature Conditions Under Historical and Future Scenarios Spatiotemporal Changes in Temperature

In the future (2015–2100), the mean temperature will rise at 0.2°C, 0.3°C, 0.5°C, and 0.6°C per decade under the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, respectively (Figure 3). In terms of the spatial distribution, SSP5-8.5 experiences a 163% increase in the grids with annual mean temperature >30°C in comparison to SPP1-2.6, which is mainly distributed in the tropical desert region of Africa. A 46% decrease was observed in the area with annual mean temperature <−20°C, which is mainly distributed in Greenland (Figure S12 in Supporting Information S1).

Figure 3. Temporal changes in the mean annual temperature during the period of 1982–2100 from the 11 earth system models. Gray = historical, pink = SSP1-2.6, blue = SSP2-4.5, green = SSP3-7.0, red = SSP5-8.5, and the shaded area indicates the standard deviation.

At the end of this century, the global temperature increases >2°C under the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios are 3%, 40%, 97%, and 99%, respectively (the Paris Agreement proposed that the global mean temperature increase should not exceed 2°C by the end of this century); warming is most obvious in the middle- and high-latitude regions (Figure 4). At present (1982–2014), over 66% of the global regions exhibit a temperature increase >0.2°C per decade; it is identified that the global temperature dynamics are in the mid-range of SPP2-4.5 (50%) and SSP3-7.0 (93%) (Figure S13 in Supporting Information S1), which is worthy of research attention.

Figure 4. Global spatial distribution of changes in the mean annual temperature between 1982–2015 and 2091–2100 in the four scenarios from the 11 earth system models. The right panel shows the changes in the mean annual temperature variation along the latitudinal gradient, the gray dotted line is global mean annual temperature variation (common to all panels).

Spatiotemporal Changes in Global Terrestrial Vegetation Productivity

The global total terrestrial GPP exhibits significant increases throughout the 21st century. In comparison to the historical period, the total GPP will increase by 6%, 11%, 14%, and 32% under the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios, respectively, which indicates that more CO₂ emissions lead to stronger increases (Figure 5). In terms of the spatial distribution, where with the higher total annual GPP is mainly distributed in parts of tropical Asia, Africa, and South America, and the vegetation in the mid- and high-latitude regions of the Northern Hemisphere (30°–75°N) will show better growth under SSP2-4.5, while the vegetation in the low-latitude regions (30°S–30°N) will yield higher GPP under SSP3-7.0 and SSP5-8.5 (Figure S14 in Supporting Information S1).

Figure 5. Temporal changes in the gross primary productivity under the four scenarios during the period of 1982–2100 from the 11 earth system models. Pink = SSP1-2.6, orange = SSP2-4.5, purple = SSP3-7.0, blue = SSP5-8.5, and the shaded area indicates the standard deviation; the green solid line indicates the mean value for the observation period from 1982 to 2014; the colored dotted line is the trend line, and the black dotted line is the mean value for the historical period (common to all panels).

By the end of the 21st century, the largest contribution to the total global terrestrial GPP growth is at low-latitudes, followed by mid-latitudes in the Southern Hemisphere (30°–60°S) (Figure 6). The global GPP values in the sustained growth regions for SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 are 83%, 82%, 98%, and 99%, respectively. During the observation period, only 55% area of the GPP is growing continuously (Figure S15 in Supporting Information S1).

Figure 6. Global spatial distribution of changes in the gross primary productivity (GPP) between 1982–2015 and 2091–2100 under the four scenarios from the 11 earth system models. The right panel shows the changes in the GPP variation along the latitudinal gradient, the gray dotted line is global mean annual GPP variation (common to all panels).

Future Vegetation–Temperature Relationship Change in Partial Correlation Caused by Future CO₂ Emissions

Spatial patterns of multi-model correlations between the GPP and temperature across the historical period and future exhibit significant spatial heterogeneity, which indicates that the vegetation in each region is severely vulnerable to climate change. During the historical period, warming promoted the growth of vegetation for regions >60°N, but induced photosynthetic suppression of vegetation in some regions at low-latitude regions (45°S–45°N), such as southern Africa, Australia, central South America, and southern North America (Figure 7a). The latitudinal distribution under SSP1-2.6 is similar to the historical period (Figure 7b); however, for the central African and peninsular Indian regions, the effect of warming on vegetation productivity is reversed, with negative effects. The warming forces on vegetation productivity weaken under SSP2-4.5 (Figure 7c), which could be attributed to the increased contribution of other drivers to vegetation growth (e.g., precipitation, soil nitrogen fixation) when CO₂ increases. The warming-induced deterioration of vegetation productivity is more severe only in southern North America and eastern South America under SSP3-7.0 and SSP5-8.5 (Figures 7d and 7e), but the spatial heterogeneity of the R_GPP-Temp distribution under SSP5-8.5 becomes smaller over global land. In comparison to the observation period (76% of global land grids), the area proportions with a significantly positive R_GPP-Temp (P < 0.05) under the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 scenarios are 42%, 83%, 75%, and 77%, respectively.

Figure 7. Spatial pattern of RGPP-Temp based on (a) observations during the period of 1982–2014 and (b–e) simulations from the 11 earth system models during the period of 2015–2100. Black dots indicate grid cells with significant correlations (P [less than] 0.05); the right panel shows the changes in the RGPP-Temp variation along the latitudinal gradient, and the bottom panel shows the changes in the RGPP-Temp variation along the longitudinal gradient.

Differences in the S_vpt Trend for Multiple Scenarios

Figure 8 depicts the long-term trends in the S_vpt for the future periods of the four scenarios. Globally, the S_vpt significantly declines by the end of the 21st century (P < 0.1), which indicates that continued warming will slow vegetation growth. However, the attenuation degrees of the S_vpt vary for the future scenarios of the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 from 2015 to 2100, and are −0.0008/decade, −0.0001/decade, −0.0019/decade, and −0.0008/decade, respectively. In addition, it is indicated that the suppression of vegetation growth by warming is not constant; the timing of the inflexion of the trend shift differs for all four scenarios. Under the SSP1-2.6 and SSP5-8.5 scenario (Figures 8a and 8d), the S_vpt inflexion points occur in 2080 and 2085, respectively, after exhibiting an increasing trend. However, negative trends will strengthen after the inflexion under the SSP2-4.5 and SSP3-7.0 scenario (Figures 8b and 8c)—in particular, under the SSP3-7.0 scenario, the S_vpt is lower than the average S_vpt of the four scenarios after 2060, suggesting a possible carbon sink-to-source tipping point in 2060.

Figure 8. Svpt trends based on the 11 Coupled Model Intercomparison Project Phase 6 models (black dashed line). The mean Svpt values under the four scenarios (blue line) and the Svpt under each scenario (red line) from 2015 to 2100: (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5.

The variability trends of the S_vpt under the four emission scenarios were also analyzed. 33%–63% of global terrestrial ecosystems change significantly and primarily exhibit negative trends (adjusted P < 0.05; Figure 9). Further, the spatial distribution of the change trends differs slightly under different scenarios. From the SSP1-2.6 to the SSP5-8.5 scenario, a negative trend of the S_vpt is shown at mid-latitudes (35°S–35°N), while change trends weaken toward the poles. Under the SSP1-2.6 scenario, a positive trend was found in the S_vpt in central Africa; however, this trend reverses as the CO₂ concentration increases, while the regions with a positive trend continue to expand in Australia. The negative trends in the Amazon, tropical Africa, and the Indian Peninsula tend to strengthen under SSP3-7.0 in comparison to other scenarios, indicating that the vegetation greening of these areas will be affected negatively by global warming. As the concentration of CO₂ emissions increases, more attention should be paid to the vegetation growth dynamics in these areas in the future.

Figure 9. Spatial patterns of the positive and negative trends in the Svpt for 2015–2100 based on the 11 Coupled Model Intercomparison Project Phase 6 models. Black dots indicate significant grid cells (adjusted P [less than] 0.05) using a Mann–Kendall test; the left panel shows the changes in the Svpt variation along the latitudinal gradient, and the bottom panel shows the changes in the Svpt variation along the longitudinal gradient: (a) SSP1-2.6, (b) SSP2-4.5, (c) SSP3-7.0, and (d) SSP5-8.5.

Comparisons of the S_vpt in Specific Vegetation Types

To further explore how vegetation itself affects the S_vpt, we relate the S_vpt to different vegetation types to analyze the sensitivity trends exhibited by each vegetation type under various scenarios (Figures 10 and 11). In this study, we observe that deciduous needle-leaf forests exhibit negative trends in the S_vpt in comparison to other vegetation types under SSP1-2.6, grasslands. Closed shrublands exhibit a decline trend in the S_vpt, and are more sensitive to warming under SSP2-4.5. Croplands exhibit negative trends, while permanent wetlands have positive trends in the S_vpt under SSP3-7.0. Deciduous needle-leaf forests exhibit greater sensitivity to warming under SSP5-8.5. In general, with increases in CO₂ emissions; significant decline trends in the S_vpt are identified for grasslands, croplands, and evergreen broadleaf forests. In comparison, more permanent wetlands, closed shrublands, and deciduous needle-leaf forests exhibit mild positive trends in the S_vpt.

Figure 10. Area percentages of positive and negative trends in the Svpt (adjusted P [less than] 0.05) across different vegetation types under the four scenarios from the 11 earth system models.

Figure 11. Area percentages of Svpt values (adjusted P [less than] 0.05) across different vegetation types under the four scenarios from the 11 earth system models.

Future Variations in S_vpt Owing to Other Forces

Although the S_vpt in this study considered temperature to be the dominant contributor to global vegetation productivity and interpreted the trends of S_vpt under different radiative forcing scenarios, other meteorological variables may also influence the bio-physical processes of vegetation (e.g., moderately favorable precipitation conditions or extreme weather). Therefore, we analyzed the potential relationship between the S_vpt and precipitation under the four scenarios to obtain a better understanding of the future predicted trends (Figure 12). Across different scenarios, the S_vpt trend shows an inverted U-shaped relationship with precipitation, that is, as precipitation increases, the S_vpt increases; however, when the precipitation exceeds a certain threshold, the positive effect of temperature on vegetation growth is suppressed. Further, the precipitation thresholds differ for all four scenarios. Under the low-emission scenario, the threshold of precipitation is around 1600 mm, while under the high-emission scenario, it is around 2000 mm. That is, the precipitation threshold is higher under the high-emission scenario, which is closely related to the CO₂ emission concentration. In the high-emission scenario, the temperature increases rapidly, and the higher temperature enhances transpiration; therefore, plants require more precipitation to maintain natural growth.

Figure 12. Relationship between the Svpt calculated from 11 earth system models and precipitation at the grid level under the four scenarios. The green points represent mean values of the Svpt and precipitation for the seven precipitation intervals, which are 0–500 mm, 500–1,000 mm, 1,000–1,500 mm, 1,500–2,000 mm, 2,000–2,500 mm, 2,500–3,000 mm, and greater than 3,000 mm. The curve fitting in panels (a–d) is based on a quadratic polynomial model.

Precipitation influences the S_vpt trends of different vegetation types, which differ in each of the four scenarios (Figure 13). The influence of precipitation on S_vpt trends increased nonlinearly under the low-emission scenario, while under middle- and high-emission scenarios, it exhibited an inverted U-shaped relationship. Furthermore, because the water demand of different types of vegetation in the process of growth and development is different, the influence of precipitation on S_vpt differs significantly among them. For grasslands, open shrublands, and deciduous needle-leaf forests, increasing in precipitation has only a slight effect on the sensitivity, while for vegetation such as permanent wetlands, mixed forests, and deciduous broadleaf forests, the contribution of precipitation to S_vpt changes was greater.

Figure 13. Relationship between the Svpt calculated from 11 earth system models and precipitation at the grid level derived from different vegetation types under four scenarios. The colored points represent the mean values of the Svpt and precipitation, the shaded areas indicate 95% confidence interval.

Although climatic dynamics are the primary factors that influence vegetation growth, it is thought that socio-ecosystems also affect it. Population pressure and economic inequality can, to a certain extent, reflect regional differences in the S_vpt variation. Through statistical calculations, it is found that countries in different income groups exhibit different trends in the S_vpt under the four radiative forcing paths (Figure 14a, Figure S16 in Supporting Information S1). It can be seen that as CO₂ emission concentrations increase, the mean S_vpt significantly weakens in low-income regions, such as countries in Central and South Africa (Figure S5 in Supporting Information S1). However, the mean S_vpt exhibits positive trends in high-income countries, such as developed countries in North America and Europe.

Figure 14. Mean Svpt trends of the (a) per income group and (b) population density region under the four scenarios.

Differences in the population density also contribute to regional changes in the S_vpt (Figure 14b, Figure S17 in Supporting Information S1). From a low CO₂ emissions path to a high energy-intensive path, lower-population regions exhibit a significant positive trend, while the mean S_vpt apparently decreases in the other three population density level regions. As population density increases, the negative trend in the S_vpt is stronger. The strongest declining trends are found in high-population areas; these countries are mainly in peninsular India and East Asia (Figure S4 in Supporting Information S1).

In general, there is a potential correlation between economic conditions and population pressure, and their forces on the S_vpt tend to overlap (Figure 15). It should be noted that a consistent positive trend in the S_vpt with rising CO₂ emissions is present in low-income countries with lower population densities, while the S_vpt trend changes from negative to positive in low-income countries with middle and high population densities. There are consistent negative trends in the S_vpt in lower middle-income countries with lower population densities from the SSP1-2.6 to SSP5-8.5 scenarios, and a similar trend is also evident in upper middle-income regions. In high-income countries, as the intensity of radiative forcing increases, the trend of S_vpt in low-population-density areas turns from negative to positive, while the negative trend in middle and high-population-density areas will strengthen in the future.

Figure 15. Mean Svpt trends per income group for different population density levels under the four scenarios.

Discussion Spatio-Temporal Change of Global Temperature and GPP

In this study, we provide a comprehensive and timely analysis regarding past and future changes in global terrestrial's climate dynamic and vegetation productivity condition in the 21st century. It is observed different warming rates under various scenarios (Fan et al., 2020). The warming rates are 0.2°C, 0.3°C, 0.5°C and 0.6°C per decade under the SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5 scenarios, respectively, not only reflect the likely trend of future climate change, but also highlight the potential impacts of different mitigation strategies. In the SSP1-2.6 scenario, a decline in global temperatures is observed in the late 21st century, validating the effectiveness of proactive climate mitigation and adaptation strategies (Gidden et al., 2019; Kriegler et al., 2014). Furthermore, the global goal outlined in the Paris Agreement is to limit the increase in global average temperature to within 2°C of pre-industrial levels (Meinshausen et al., 2009). However, the IPCC Sixth Assessment Report indicated that to achieve this goal, the global cumulative emissions allowed between 2020 and 2100 are 2.3 trillion tons, with only a 17% chance of success (DeAngelo et al., 2021). Currently, global carbon emissions have exceeded 50% of the carbon budget. At the current emission rate, it is projected that the world will exhaust the remaining emission budget within the next 30 years (Lamboll et al., 2023). To address this urgent challenge, 97 Parties, accounting for about 82% of global greenhouse gas emissions, have already made net-zero emission commitments (Hoehne et al., 2021). If national targets for net-zero emissions are achieved in full, projected global average temperature rise could be limited to an estimated 2.0–2.4°C by the year 2100 (Hoehne et al., 2021). However, to move from current emission levels to higher emission reduction targets in the future and eventually transition to net-zero emissions, countries need to take more decisive and effective measures. This includes accelerating energy transition, developing low-carbon technologies, improving energy efficiency, and implementing comprehensive landscape management.

For terrestrial ecosystems, we predicted that most terrestrial vegetation productivity will continue to increase globally by the end of the 21st century, especially in medium- and high-emission scenarios, and particularly in low-latitude regions such as the Amazonian forests and the tropical forests of Africa. This possibly link to the CO₂ fertilization effect, where elevated atmospheric CO₂ concentrations stimulate plant photosynthesis rates, thus enhancing productivity (Keenan et al., 2023). In high-latitude regions, the slow increase in vegetation productivity may relate to earlier tree activity or prolonged growing season by warmer temperature (S. Gao et al., 2022; Grossiord et al., 2022; Wang et al., 2021). On the other hand, some studies have pointed out that the irrational logging and wildfires could cause the reduction of forest cover in the future, which poses a pronounced threat to the improvement of global vegetation productivity (Gatti et al., 2021; Hubau et al., 2020). To mitigate natural disasters, urban expansion, and the adverse effects of climate change, enhancing the resilience and recovery capacity of ecosystems is crucial. Some regions of the world have implemented policy measures, such as China's Grain for Green program, Brazil's Amazon Rainforest Protection Plan, and Africa's Great Green Wall initiative, which have effectively enhanced the carbon sequestration function of ecosystems. The future global vegetation carbon cycle will depend on the interactions between these negative and positive factors and how we manage and protect our natural resources.

Response of GPP to Global Warming

Global environmental changes, particularly climate warming, are rapidly altering the growth condition of terrestrial vegetation (Yuan et al., 2019), disrupting the carbon balance of the pre-existing terrestrial ecosystem and threatening global socio-ecological systems sustainable development (Piao et al., 2020; Shi et al., 2021). Thus, the response of vegetated ecosystems to global warming has become an important research component with fruitful results (Alkama et al., 2022; X. F. Gao et al., 2021; Wu et al., 2017). However, the time scale of these studies is mainly focused on the last a few decades, it is still unclear about the global warming trend, changes in vegetation productivity and its sensitivity to warming under different CO₂ emission scenarios in the future. Accordingly, we profile the S_vpt of different vegetation types under different CO₂ emission scenarios in this century. It is found that the global S_vpt significantly declines by the end of the 21st century, suggesting that future global warming will cause vegetation growth to slow down, and vegetation function will deteriorate persistently (Zhang et al., 2022). At the same time, the sensitivity of vegetation to warming varies significantly under different CO₂ emission intensities. The S_vpt under the SSP1-2.6 and SSP5-8.5 scenarios show a decreasing and then increasing trend. After 2080 and 2085, respectively, the S_vpt values are consistently greater than the global mean, indicating that the warming effect on vegetation growth is enhanced; The S_vpt under the SSP2-4.5 and SSP3-7.0 scenarios will keep weakening, especially under the SSP3-7.0 scenario, the S_vpt value will be less than the global mean after 2060, and there will be a reversal trend from carbon sink to carbon source, particularly in the Amazon rainforest, tropical Africa and southern North America, where continued warming will restrain vegetation growth, and more attention should be paid to the response of vegetation growth to climate change in these regions in the future.

In this study, by comparing the differences in the S_vpt between different vegetation types, it is found that future warming has an enhanced positive effect on the growth of permanent wetlands, closed shrublands, and deciduous needle-leaf forests. While the positive effect of warming on the growth of vegetation—such as grasslands, croplands, and evergreen broadleaf forests—weakened (O'Sullivan et al., 2020). On the one hand, the difference in species' structure and biomes lead to differences in the thresholds of the temperature response of different vegetation types (Piao et al., 2014, 2019; Yuan et al., 2021). As indicated by some research, the conclusion that—in southwestern Africa—elevated temperatures have hampered the survival of larger woody plants and facilitated the spread of smaller ones validates the findings of this study (Venter et al., 2018; Zhang, Brandt, et al., 2019). On the other hand, vegetation growth is actually subject to the interactions of multiple factors (Zhu et al., 2016). Warming alone may not lead to significant changes in vegetation productivity, while extreme droughts caused by warming may trigger severe vegetation decline and even death (Allen et al., 2015; Xu et al., 2019); this is confirmed through the analysis of the effect of precipitation on the S_vpt in this study. It is also found that the impact of precipitation on the S_vpt is nonlinear, with a significant threshold. Before reaching the threshold, the “synergistic effect” of precipitation and temperature will promote the vegetation photosynthesis (Gonsamo et al., 2021); while after the threshold, excessive precipitation can cause the problems of hypoxia, nutrient loss, and root rot, and inhibit vegetation growth (Chen et al., 2022). Extensive studies have indicated that the arrival of climate thresholds can result in biodiversity loss, decreased carbon sink capacity, and accelerated climate change (Li et al., 2023a, 2023b). The study of Doughty et al. (2023) on tropical forests noted that if temperatures rise by 4°C, it could stop photosynthesis in leaves, leading to death. Berdugo et al. (2020) predicted that under the SSP5-8.5 scenario, the productivity and richness of plants in over 20% of global drylands will surpass drought thresholds. The results of this study showed that the precipitation threshold is higher under the high-emission scenario, that's because vegetation in the high-emission scenario faces greater drought risks and intensity, and vegetation growth requires more water (Maurer, et al., 2020). Therefore, based on the differences in the S_vpt across vegetation types and the effect of precipitation on the S_vpt changes, suitable land management, water resource management, and ecological restoration policies can be used to offset the adverse effects of global warming on vegetation growth. On the one hand, drought-resistant tree species, mixed planting and intercropping, and landscape configuration can be adopted to improve the water use efficiency of vegetation ecosystems (Liu, Guan, et al., 2021; Vitali et al., 2018); on the other hand, water storage facilities can be built manually to increase the spatiotemporal water availability to alleviate the rigid constraints of water resources on vegetation (Labonte-Raymond et al., 2020).

Potential Relationship Between Socio-Economy and S_vpt

In addition to climate change, socio-economic factors such as population pressure and income levels are also closely related to the S_vpt changes (Venkatesh et al., 2022). The results indicate that under the low-emission scenario, the S_vpt in sparsely populated areas with lower income generally presents a positive trend, that is because the region is subjected to less human disturbance (Ge et al., 2021), and vegetation growth is mainly driven by the positive effect of climatic factors (Philippon et al., 2014). In contrast, under the low- and medium-emissions scenario, the S_vpt in low-income areas with higher population densities always exhibits a negative trend, mainly because increasing population pressures can aggravate the overexploitation of resources, such as deforestation and the shortening of the fallow period, ultimately leading to land degradation. Meantime, with the increase of CO₂ concentration, there is a clear “U-shaped” relationship between different population densities (income levels) and S_vpt, indicating a high S_vpt under SSP1-2.6 and SSP5-8.5 scenarios, while a low S_vpt under SSP2-4.5 and SSP3-7.0 scenarios. For a heavily transient climate, the oceans temporarily protect by surpressing warming. What this implies is that CO₂ can be ahead of any detrimental warming, and so vegetation responds well to the fertilization (Tian et al., 2021). This could be what is causing the high productivities under SSP5-8.5. At the SSP1-26, less land use/cover changes and appropriate temperature might be good for vegetation as well (Hou et al., 2022). Then in the medium, for SSP2-4.5 and SSP3-7.0, there is substantial warming, that “catches up” with the higher CO₂ concentration, and the respiration and transpiration caused by this offsets the positive effect of fertilization (Liu, He, et al., 2021). Furthermore, as CO₂ concentrations rise, the S_vpt tends to improve in the areas with low population densities and high income levels, mainly because developed countries have easier access to the economic resources and technological conditions required for ecological restoration and are, thus, better able to cope with the threats to vegetation from warming (Abel et al., 2021). On the whole, compared to high-income countries, vegetation growth in low-income countries will be more vulnerable to the threat of global warming (Diffenbaugh & Burke, 2019), and it is necessary to develop effective measures to enhance the resilience of socio-ecological system to resist climate changes.

Specifically, for the low-income countries, emphasis should be placed on the transformation and optimization of traditional high-carbon industries to improve energy utilization efficiency and reduce high-carbon energy consumption (Sun et al., 2017); at the same time, high-temperature resistant crop can be developed, and the ecological restoration and governance should be actively carried out to improve vegetation coverage area and environmental quality (Tong et al., 2020; Tribouillois et al., 2016). For the high-income countries, low-carbon new energy development and technological innovation should be strengthened to form a complete low-carbon technology R&D system (Shi et al., 2020); a stricter low-carbon living management systems should be implemented and the high-temperature detection and warning systems should be established to reduce unnecessary socio-economic losses (Lorenz et al., 2019). At the same time, international cooperation especially the technology transfer and knowledge sharing should be strengthened, and measures to cope with global warming should be jointly formulated; moreover, the carbon emission trading mechanism and compensation system among different countries should be improved (Zheng, 2023).

Limitation and Prospect

CMIP6 includes a series of ESMs, and there is still uncertainty in the application of these ESMs. For example, some terrestrial carbon cycle models only include the one-way process of the climate impacting GPP (Huang et al., 2021); due to the differences in parameterization and structure among different modes, there is significant internal variation when simulating the same index (Anav et al., 2013; Zhu et al., 2016); most ESMs fail to capture the temporal variability in the annual vegetation index, such as leaf area index because of large biases in both the simulated trend magnitude and temporal pattern of interannual variability (Song et al., 2021). And such biases may have further impacts on the simulations of the related vegetation productivity, carbon and land-atmosphere interaction variables for global change research. Hence, bias-correction should be made to improve reliability of vegetation structure and distribution when future projections and historical reconstructions. This study adopted the quantile mapping method for the bias-correction and achieved the ideal result. However, the selection of the correction method may affect the accuracy and applicability of the corrected data. Therefore, different bias-correction methods will be compared to obtain the optimal simulation results in the future. This study multi-mode ensemble averaging method to calculate the arithmetic mean of all model simulation values to reduce the uncertainty of multi-model simulations. Some studies have confirmed that the weight determination method based on the performance of simulating historical scenarios and the convergence degree of future predictions for each model may have a better result (Giorgi & Mearns, 2002). However, having too many sets of models can lead to computational and administrative challenges, and the difference in model predictions for different variables can also lead to weighting differences. It would be great if in the future some kind of constraints could be developed to relate current observations to the model's estimates, identifying which model is most trustworthy for a given variable. Moreover, because the spatial projection data for other environmental factors in CMIP6 models is missing—including the soil nutrient content (e.g., phosphorus, magnesium)—only the potential relationship between the S_vpt and a portion of the driving forces was analyzed in this study, and the mechanisms of multi-driving force interactions spatially is not considered. In the future, it is necessary to strengthen the comprehensive analysis of the driving factors for global S_vpt changes.

Conclusion

Vegetation productivity will increase in most regions across the globe for a long time to come, especially at low latitudes, such as the Amazon, central Africa, and Indonesia. Under the four CO₂ emission scenarios in the future, 33%–63% of global terrestrial ecosystems will exhibit significant changes in the S_vpt (adjusted P < 0.05). Overall, the global sensitivity of vegetation to temperature will decrease, but the suppression of vegetation growth by warming is not constant; this trend can be reversed or strengthened after inflexion occurs. It is found that carbon sink-to-source reversal may occur after 2060 under the SSP3-7.0 scenario. For different vegetation types, such as grasslands, croplands, and broadleaf forests, there will be a reduced availability of heat; in contrast, wetlands, shrublands, and coniferous forests will exhibit greater productivity per unit heat as warming occurs. Further, precipitation can impact the S_vpt trend; across different scenarios, the S_vpt trend shows an inverted U-shaped relationship with precipitation. In addition to climate, the S_vpt is also closely linked to social systems; growing population pressure will reduce carbon sinks during warming. However, the economic conditions and social systems of each country differ; to better address climate change and promote socio-ecological development in the future, it is necessary to formulate corresponding measures in terms of land use and population management according to national conditions.

Acknowledgments

This article is funded by the National Key Research and Development Program of China (No. 2022YFF1300701) and National Natural Science Foundation of China (No. 42001090).

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

The GPP and climate simulation data for historical and future period can be accessed from the CMIP6 archive (https://esgf-node.llnl.gov/search/cmip6/); Global monthly temperature and precipitation data from the Climatic Research Unit gridded Time Series Version 4 (CRU TS4.05) data set can be gained from the Climatic Research Unit at the University of East Anglia (UEACRU et al., 2021); The FLUXCOM GPP data provided by Jung et al. (2019) was publicly available for download from the Data Portal of the Max Planck Institute for Biogeochemistry (https://doi.org/10.17871/FLUXCOM_RS_METEO_CRUNCEPv6_1980_2013_v1); The Gridded Population of the World, Version 4 (GPWv4) was provided by Center for International Earth Science Information Network at Columbia University (CIESIN, 2018); Global income groups data can be gained from the World Bank (2015); The MODIS land cover product (MCD12C1) was available from Friedl and Sulla-Menashe (2019).