How the Updated Earth System Models Project Terrestrial Gross Primary Productivity in China under 1.5 and 2 ◦ C Global Warming

: Three Earth system models (ESMs) from the Coupled Model Intercomparison Project phase 6 (CMIP6) were chosen to project ecosystem changes under 1.5 and 2 ◦ C global warming targets in the Shared Socioeconomic Pathway 4.5 W m − 2 (SSP245) scenario. Annual terrestrial gross primary productivity (GPP) was taken as the representative ecological indicator of the ecosystem. Under 1.5 ◦ C global warming, GPP in four climate zones—i.e., temperate continental; temperate monsoonal; subtropical–tropical monsoonal; high-cold Tibetan Plateau—showed a marked increase, the smallest magnitude of which was around 12.3%. The increase was greater under 2 ◦ C of global warming, which suggests that from the perspective of ecosystem productivity, global warming poses no ecological risk in China. Speciﬁcally, in comparison with historical GPP (1986–2005), under 1.5 ◦ C global warming GPP was projected to increase by 16.1–23.8% in the temperate continental zone, 12.3–16.1% in the temperate monsoonal zone, 12.5–14.7% in the subtropical–tropical monsoonal zone, and 20.0–37.0% on the Tibetan Plateau. Under 2 ◦ C global warming, the projected GPP increase was 23.0–34.3% in the temperate continental zone, 21.2–24.4% in the temperate monsoonal zone, 16.1–28.4% in the subtropical–tropical monsoonal zone, and 28.4–63.0% on the Tibetan Plateau. The GPP increase contributed by climate change was further quantiﬁed and attributed. The ESM prediction from the Max Planck Institute suggested that the climate contribution could range from − 12.8% in the temperate continental zone up to 61.1% on the Tibetan Plateau; however, the ESMs differed markedly regarding their climate contribution to GPP change. Although precipitation has a higher sensitivity coefﬁcient, temperature generally plays a more important role in GPP change, primarily because of the larger relative change in temperature in comparison with that of precipitation.


Introduction
Gross and net primary productivity (GPP and NPP, respectively) are representative indicators that reflect ecosystem production capacity [1][2][3]. Many previous studies have considered future GPP/NPP change. For example, Huang et al. [4] evaluated NPP variations in the 21st century under various climate scenarios using the Lund-Potsdam-Jena dynamic global vegetation model. They found that total NPP in China is projected to increase continuously under different scenarios, with CO 2 concentration playing the dominant role. Using a machine learning model to constrain the spatial uncertainty in GPP projections, Schlund et al. [5] predicted a higher increase in GPP in northern high latitudes over the 21st century under the Representative Concentration Pathway [6] 8.5 W m −2 (RCP8.5) in comparison with regions closer to the equator. Under 1.5 • C of global warming, the GPP in China is expected to increase by 15.5% ± 5.4% on a stabilized pathway and by 11.9% ± 4.4% on a transient pathway [7]. Zhang et al. [3] explored the trend features of GPP/NPP in the 21st century under the Shared Socioeconomic Pathway [8] 24.5 W m −2 (SSP245) with the Beijing Climate Model. Their results predicted the overall trends of increase in both the near-term and long-term terrestrial GPP/NPP. However, in certain districts, the trend of GPP/NPP showed an initial increase followed by a decrease. Wang et al. [9] investigated the variation in NPP over the 21st century using the Earth system models (ESMs) of the Coupled Model Intercomparison Project phase 5. The results obtained under all four RCP scenarios suggested an increasing trend of NPP over China, especially in western areas.
In summary, GPP/NPP in China under different scenarios is expected to show a trend of increase in the 21st century. However, large uncertainties exist in the various ESMs [5,9]. Under the global warming targets of 1.5 and 2 • C above preindustrial levels set by the Paris Agreement, many regional impacts wait to be assessed. In particular, as the Coupled Model Intercomparison Project enters into the 6th phase (CMIP6), more and more ESMs have distributed their latest climate simulation under the Shared Socioeconomic Pathways (SSPs). How the latest ESMs will project the future ecosystem change in China and the corresponding climate attribution remains to be determined and revealed. In a recent study on performance in presenting historical terrestrial GPP in China, three out of seven ESMs evaluated were found to perform well in terms of climatological GPP, spatial pattern, and the ecosystem-climate relationship [10]. Consequently, these three ESMs were chosen in this study to predict ecosystem change under the warming targets. The ecological indicator of annual GPP was applied to measure the general state of the ecosystem. Changes in annual GPP predicted using the different ESMs were quantified with respect to the different climate zones in China. Furthermore, the relationship between the ecosystem and climate variables was tested and built through linear correlation and multiple regression. Relying on the model-specific parameters of the ecosystem response to the climate, the climate-related GPP changes were revealed and quantified.

Data
In this study, three ESMs that performed well in historical GPP reproduction [10] were chosen to project future GPP in China: (1) the Beijing Climate Center Climate System Model (BCC-CSM2-MR) [11], (2) the Euro-Mediterranean Centre on Climate Change coupled climate model (CMCC-CM2-SR5) [12], and (3) the Max Planck Institute for Meteorology Earth System Model version 1.2 (MPI-ESM1.2-HR) [13]. Specifically, BCC-CSM2-MR and MPI-ESM1-2-HR, out of seven ESMs, gave the best estimation of climatological GPP in China from 1980 to 2013. MPI-ESM1-2-HR performed best in characterizing the spatial structure. BCC-CSM2-MR and CMCC-CM2-SR5 best captured the response of the ecosystem to the climate [10]. The land surface models used for the three ESMs were BCC-AVIM2.0, CLM4.5, and JSBACH3.2. Major improvements or parameterizations have been made to these models in comparison with their predecessors; they make use of new scientific understanding to better simulate vegetation phenology [12][13][14]. These ESMs could provide not only the monthly GPP, but also the monthly surface air temperature and precipitation. In CMIP6, new SSPs were employed for climate modelling. The SSPs included five narratives describing alternative socio-economic developments, such as sustainable development, fossil-fueled development, etc. [8]. The middle of the road development-i.e., the SSP2 scenario-featured a continuation of the current fossil fuel-dominated energy mixed with intermediate challenges for both mitigation and adaptation, which resembled the historical pattern most [8]. SSP245, as the sole scenario of SSP2 implemented in CMIP6, was thus chosen to represent the most possible future world. Historical data  were applied to determine the ecosystem-climate relationship-i.e., correlation and multiple linear regression. Data from the BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1.2-HR ESMs were output as 1.125 • × 1.125 • , 0.9375 • × 1.25 • , and 0.9375 • × 0.9375 • grids, respectively. Because the grids were not uniform, they were first transformed to a 1 • × 1 • grid through bilinear interpolation for comparative purposes.
A climate division map of China was applied for regional analysis. It divided the country into four climate zones-i.e., temperate continental, temperate monsoonal, subtropicaltropical monsoonal, and high-cold Tibetan Plateau, as in He et al. [15] and Zhang et al. [16].

Bilinear Interpolation
Bilinear interpolation can produce a smoother interpolation than that achieved using the nearest neighbor method [17]. Thus, it was applied to transform fields from various grids of the ESMs into the formal 1 • × 1 • grid. In this approach, g (n 1 , n 2 ) is defined as a linear combination of the values of its four nearest neighbors. Given the four neighbors with coordinates f (n 10 , n 20 ), f (n 11 , n 21 ), f (n 12 , n 22 ), and f (n 13 , n 23 ) (i.e., the four nearest neighbors of f (n 1 , n 2 )), the geometrically transformed field g (n 1 , n 2 ) is computed as: (1) The bilinear weights A 0 , A 1 , A 2 , and A 3 are found by solving: 1 n 10 n 20 n 10 n 20 1 n 11 n 21 n 11 n 21 1 n 12 n 22 n 12 n 22 1 n 13 n 23 n 13 n 23

Area Weighting
Regional and global mean variables-e.g., temperature, precipitation, and GPP-on the 1 • × 1 • grid are calculated through area weighting: where θ represents the latitude of the grid, R is the Earth's radius, and V is the variable.

Linear Correlation and Multiple Regression
The correlation coefficient r is used to test the relationship between ecosystem productivity and climate factors. The formula can be expressed as follows: where E and C represent ecosystem productivity and climate factors, respectively. The interannual variation in GPP reflects year-to-year differences attributable mainly to climate variations [18,19]; therefore, the relationship between GPP and climate-i.e., precipitation and surface air temperature-was explored using a multiple regression approach [20]: where y is the detrended anomaly of the carbon flux GPP, variable x T is the detrended annual temperature anomaly, and x P is the detrended annual precipitation anomaly. The fitted regression coefficients a and b define the apparent carbon flux sensitivity to interannual variations in temperature and precipitation, and ε is the residual error term. The use of the detrended time series instead of the original nonstationary time series in the above linear correlation and regression analysis provides a robust estimate of their relationship [21][22][23]. Some definitions set 1986-2005 as a reference period when the global surface air temperature was 0.61 • C warmer than preindustrial levels [24,25]. We adopted this definition and defined the 1.5 and 2 • C warming periods as the first time when the 20-year-moving-average global temperature was 0.89 and 1.39 • C warmer, respectively, than that from 1986-2005 in the models. The corresponding changes in ecosystem and regional climate were based upon the reference period of 1986-2005. It is also necessary to point out that 20 years is a duration that is commonly applied to represent a climate state in the scientific world [26][27][28].

GPP Distribution and Projected Changes
The climatological GPP distribution produced by each of the three ESMs from 1986 to 2005 is shown in Figure 1. The ESMs all produce a similar spatial pattern of GPP, showing high (low) values in the southeast (northwest) of China. Regionally, the GPP in the subtropicaltropical monsoonal zone is largest, followed in descending order by the temperate monsoonal zone, Tibetan Plateau, and temperate continental zone. The three ESMs produced comparable estimates in the climate zones except in the monsoonal regions, where CMCC-CM2-SR5 produced larger estimates, especially in the subtropical-tropical monsoonal zone. The GPP change under 1.5 °C of global warming is shown in Figure 2. Throughout China, the GPP of all three ESMs showed a positive anomaly except over certain individual grid points. In the subtropical-tropical monsoonal with MPI-ESM1-2-HR, the negatively changed grids tend to concentrate (Figure 2e). The areal GPP reduction may be related to reduced local precipitation. However, many factors could contribute to the GPP change in addition to precipitation and temperature, such as land-use change, soil mois- The GPP change under 1.5 • C of global warming is shown in Figure 2. Throughout China, the GPP of all three ESMs showed a positive anomaly except over certain individual grid points. In the subtropical-tropical monsoonal with MPI-ESM1-2-HR, the negatively changed grids tend to concentrate (Figure 2e). The areal GPP reduction may be related to reduced local precipitation. However, many factors could contribute to the GPP change in addition to precipitation and temperature, such as land-use change, soil moisture, wind speed, humidity, solar radiation, nitrogen deposition, etc. Thus, it is really hard to be conclusive. Moreover, the aggregated negative grids tend to dissipate under 2 • C global warming (Figure 3e). The GPP change patterns differ among the models (Figure 2). For example, the largest anomaly in the output of BCC-CSM2-MR appears over the southeastern Tibetan Plateau, whereas the largest anomalies in the output of CMCC-CM2-SR5 and MPI-ESM1-2-HR appear in the central and southern parts of the subtropical-tropical monsoonal zone, respectively. The absolute GPP change is largest in the subtropical-tropical monsoonal zone; however, the relative change is rather small and stable among the models (Figure 2 right). The absolute change is smallest in the temperate continental zone owing to its low base value in GPP. The relative change is large over the Tibetan Plateau, and there are strong differences in the magnitude of the GPP increments among the different models.
The GPP anomaly under 2 • C of global warming shows a spatial pattern similar to that found under 1.5 • C global warming but with a stronger intensity ( Figure 3). Regional statistics indicate that the regional GPP changes will be larger under 2 • C of global warming. The projected increment of GPP in China under the different warming targets is consistent with previous findings [5,7,29]. This suggests that from the perspective of GPP, there is no ecological crisis in the projected future climate within the studied domain [3]. As with 1.5 • C warming, the subtropical-tropical monsoonal zone with the highest GPP value contributed the most to the increment in China's GPP under 2 • C of warming. However, the increase rate does not show much difference in magnitude when compared with that of other regions. It is worth noting that the rate of increase in GPP is substantial on the Tibetan Plateau-i.e., the increase is nearly 63% with regard to BCC-CSM2-MR. Thus, the Tibetan Plateau would appear to be the region most susceptible to the effects of climate warming, although the influence could be considered positive and beneficial.
The seasonal GPP anomalies under the 1.5 and 2 • C warming scenarios are shown in Figures 4 and 5, respectively. The spatial modes between the two warming scenarios are similar, noting that the magnitude in the 2 • C warming is much larger than in the 1.5 • C warming. In spring and summer, the GPP anomalies are the most prosperous, as they correspond to the growing season in China, while they drop to become the weakest in winter. In spring and winter with all ESMs, the GPP all over China generally shows a positive anomaly, with only sporadic negative points. Some negative changes occur in summer and autumn, especially with BCC-CSM2-MR and MPI-ESM1-2-HR. For BCC-CSM2-MR in summer, the negative GPP anomalies concentrate in the Huaihe River, which divides the subtropical-tropical monsoonal and temperate monsoonal regions. However, there were no negative changes in GPP at the zone scale. For MPI-ESM1-2-HR in summer and autumn, we observed some negative changes over the grassland in the temperate continental, which is similar to the results of Ma et al. [30]. They found that large areas in Northern China showed a decreasing trend in NPP under global warming, although the overall NPP increased significantly. The fact that only one ESM obtained similar results also indicates the large inter-model spread in representing the future GPP change. The negative changes in the temperate continental were weakened under the 2 • C warming scenario ( Figure 5).    The GPP anomaly under 2 °C of global warming shows a spatial pattern similar to that found under 1.5 °C global warming but with a stronger intensity (Figure 3). Regional statistics indicate that the regional GPP changes will be larger under 2 °C of global warming. The projected increment of GPP in China under the different warming targets is consistent with previous findings [5,7,29]. This suggests that from the perspective of GPP, there is no ecological crisis in the projected future climate within the studied domain [3]. As with 1.5 °C warming, the subtropical-tropical monsoonal zone with the highest GPP value contributed the most to the increment in China's GPP under 2 °C of warming. How-

Climate Attribution
The variation in GPP is closely related to climate [20,31], and the correlation parameters within the studied ESMs are provided in Table 1. It can be seen that GPP is correlated significantly with at least one climate variable. There are cases in which GPP correlates negatively with temperature, such as in the temperate continental zone with MPI-ESM1-2-HR and in the temperate monsoonal zone with BCC-CSM2-MR, which imply inherent differences in ecological modeling between the different ESMs [9,10]. On the Tibetan Plateau, it is unanimous within the ESMs that GPP is most closely related to temperature rather than to precipitation. Because the ESMs substantially overestimate precipitation over the Tibetan Plateau [32][33][34], it is possible that precipitation is not the primary climate factor constraining the regional ecosystem.  The corresponding climate changes-i.e., precipitation and surface air temperature-under the warming targets are shown in Tables 2 and 3, respectively. The mode of temperature change in the four regions is consistent among the models. The hottest region-i.e., the subtropical-tropical monsoonal zone-increases least under the effects of warming. Both CMCC-CM2-SR5 and MPI-ESM1-2-HR produced similar estimates of temperature change, whereas the estimates from BCC-CSM2-MR were larger, especially over the monsoonal regions.

Climate Attribution
The variation in GPP is closely related to climate [20,31], and the correlation parameters within the studied ESMs are provided in Table 1. It can be seen that GPP is correlated significantly with at least one climate variable. There are cases in which GPP correlates negatively with temperature, such as in the temperate continental zone with MPI-ESM1-2-HR and in the temperate monsoonal zone with BCC-CSM2-MR, which imply inherent differences in ecological modeling between the different ESMs [9,10]. On the Tibetan Plateau, it is unanimous within the ESMs that GPP is most closely related to temperature rather than to precipitation. Because the ESMs substantially overestimate precipitation over There is a greater model variety regarding the change in precipitation. Under 1.5 • C of warming, there are negative changes-e.g., in the temperate monsoonal zone with CMCC-CM2-SR5 and in the subtropical-tropical monsoonal zone with both BCC-CSM2-MR and CMCC-CM2-SR5. Conversely, under 2 • C of warming, there are no negative changes, but the incremental differences for one certain region are huge. These results indicate the large uncertainty in the precipitation projections made by the ESMs. It is also worth noting that in comparison with their variabilities, the change in temperature under the warming scenarios is reasonably large, whereas the precipitation change is rather limited [10].
The apparent sensitivity of climate to the ecosystem of each of the ESMs is shown in Table 4. The response of the ecosystem to climate varies strongly among the models. For some ESM regions, climate plays a crucial role, such that the variation in climate explains more than half of the variation in GPP. However, for certain other ESM regions, the degree of explanation attributable to climate is rather small-e.g., MPI-ESM1-2-HR in the subtropical-tropical monsoonal zone and BCC-CSM2-MR on the Tibetan Plateau. For one particular region, the same climate factor might affect the ecosystem differently in the various models. Taking the subtropical-tropical monsoonal zone as an example, precipitation is the major influencing factor and affects the ecosystem positively with BCC-CSM2-MR and CMCC-CM2-SR5. However, with MPI-ESM1-2-HR, the correlation between precipitation and GPP is insignificant and negative (Tables 1 and 4). Moreover, the overall climate contribution to ecosystem variation with MPI-ESM1-2-HR is very small ( Table 4). These features further reflect the inherent differences of ecological modeling within the ESMs. Table 1. Correlation parameter r between GPP and climate variables during the historical period of 1980-2013. * denotes correlation that is significant at the 0.1 level; ** denotes correlation that is significant at the 0.05 level; *** denotes correlation that is significant at the 0.01 level.   Based on the ecosystem-climate relationship (Table 4) and the known climate changes (Tables 2 and 3), the GPP change over the climate zones with the different ESMs is attributed quantitatively in Figure 6. As mentioned before, some contributions from the climate factors are negative, for which there are two major reasons. First, the climate change is negative-e.g., the negative precipitation change leads to a negative contribution to GPP. Second, the correlation between the climate factor and GPP is negative-e.g., temperature and GPP in the temperate continental zone in MPI-ESM1-2-HR. A positive anomaly in temperature could also lead to GPP reduction. The climate contribution to the variation in GPP changes among the studied ESMs. Even with the same model-e.g., MPI-ESM1-2-HR-it can be −12.8% in the temperate continental zone and 61.1% on the Tibetan Plateau. Generally, under the effects of global warming, the influence of temperature on the ecosystem is larger than that of precipitation. This is mainly because the relative change in temperature is much larger than that in precipitation. On the Tibetan Plateau, where temperature is the most constraining factor (Tables 1 and 4), temperature plays a more dominant role than precipitation in the increase in GPP ( Figure 6). In addition, it is observed that the relative GPP increase in the Tibetan Plateau is much larger than that in other regions. On the one hand, it is related to the low baseline value of GPP in the Tibetan Plateau (Figure 1d). A light increase in the GPP of the Tibetan Plateau is salient in relative values compared to the respective change in the subtropical-tropical monsoonal region. On the other hand, this may be related to the vegetation structure on the Plateau. Demonym plants can be divided into three types based on their photosynthesis patterns (i.e., C3, C4, and crassulacean acid metabolism). C3 photosynthesis produces a threecarbon compound during the Calvin cycle, while C4 photosynthesis makes an intermediate four-carbon compound that splits into a three-carbon compound for the Calvin cycle. They favor different conditions of nature. The conditions on the frigid Tibetan Plateau are unsuitable for the growth of C4 plants [35]. Consequently, the plateau is dominant by C3 plants [36]. C3 plants are more efficient in vegetative growth than C4 plants in response to the increasing air CO 2 [37]. As a result, GPP increases more rapidly with increased air CO 2 in the Tibetan Plateau than in other regions containing both C3 and C4 plants.

Conclusions
To evaluate the GPP change under 1.5 and 2 °C of global warming, this study selected three CMIP6 ESMs (i.e., BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1-2-HR) that performed well in historical GPP modeling; the principal conclusions derived are as follows:

Conclusions
To evaluate the GPP change under 1.5 and 2 • C of global warming, this study selected three CMIP6 ESMs (i.e., BCC-CSM2-MR, CMCC-CM2-SR5, and MPI-ESM1-2-HR) that performed well in historical GPP modeling; the principal conclusions derived are as follows: 1.
Under 1.5 and 2 • C of global warming, the projections of the ESMs indicate that global warming introduces no ecological risk in China. Although certain individual grid points showed negative GPP changes, regional GPP showed a marked increase, the smallest magnitude of which was more than 10% greater than that from 1986 to 2005.

2.
Specifically under 1.5 • C warming, the GPP in the temperate continental zone is projected to increase by 16. Climate change is projected to contribute positively to GPP change, except in the temperate continental zone with MPI-ESM1-2-HR. Although precipitation has larger sensitivity parameters, temperature generally plays a more important role in GPP change because of the larger change relative to its own variability in comparison with that of precipitation.
The output of the three studied ESMs showed a marked spread, not only in GPP change but also in the accountability of climate to the ecosystem. In addition, the change in climate, especially precipitation, differed strongly within the models, which indicates the large uncertainty in the climate projections of the ESMs. All of these add to difficulties in attributing future GPP change to climate. Moreover, this study analyzed the influence of annual precipitation and temperature upon the ecosystem productivity. However, GPP variation depends not only on these, but also on wind speed, humidity, solar radiation, nitrogen deposition, etc. Future studies should be more comprehensive in building the regressed equations between GPP and the impact factors. This study failed to analyze the contribution of CO 2 to the GPP increase. This was due to the lack of gridded/regional CO 2 concentration data. Future research should take into account the CO 2 effect when analyzing the GPP change and be more specific about vegetation of C3 and C4 types. This study is more general in that it focuses on the general productivity of the climate zones. Future studies should be refined to specific vegetation covers, such as forest, grass, etc. The fact that different ESMs lack consensus in the response mechanism of the ecosystem to climate, even over one specific climate zone, highlights that there is still a long way for ecological modeling in China to go.