*Article* **Effects of Environmental Factors on the Changes in MODIS NPP along DEM in Global Terrestrial Ecosystems over the Last Two Decades**

**Zhaoqi Wang 1,\*, Hong Wang 1, Tongfang Wang 1, Lina Wang 1, Xiaotao Huang 2,3, Kai Zheng <sup>1</sup> and Xiang Liu <sup>1</sup>**


**Abstract:** Global warming has exerted widespread impacts on the terrestrial ecosystem in the past three decades. Vegetation is an important part of the terrestrial ecosystem, and its net primary productivity (NPP) is an important variable in the exchange of materials and energy in the terrestrial ecosystem. However, the effect of climate variation on the spatial pattern of zonal distribution of NPP has remained unclear over the past two decades. Therefore, we analyzed the spatiotemporal patterns and trends of MODIS NPP and environmental factors (temperature, radiation, and soil moisture) derived from three sets of reanalysis data. The moving window method and digital elevation model (DEM) were used to explore their changes along elevation gradients. Finally, we explored the effect of environmental factors on the changes in NPP and its elevation distribution patterns. Results showed that nearly 60% of the global area exhibited an increase in NPP with increasing elevation. Soil moisture has the largest uncertainty either in the spatial pattern or inter-annual variation, while temperature has the smallest uncertainty among the three environmental factors. The uncertainty of environmental factors is also reflected in its impact on the elevation distribution of NPP, and temperature is still the main dominating environmental factor. Our research results imply that the carbon sequestration capability of vegetation is becoming increasingly prominent in high-elevation regions. However, the quantitative evaluation of its carbon sink (source) functions needs further research under global warming.

**Keywords:** net primary productivity (NPP); global warming; digital elevation model (DEM); uncertainty

#### **1. Introduction**

Since the nineteenth century, the global near-surface temperature has continued to increase according to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC). The temperature in the last 10 years after 2000 resulted in the hottest 10 years in history, and 1983–2012 was the hottest 30 years in the past 800 years. The widespread impact of global warming has caused a series of negative ecological consequences, such as drought [1,2], melting [3], rising sea levels [4], and frequent extreme climates [5,6]. The accelerating global warming has become a major challenge that is restricting the sustainable development of human society [7–9].

Vegetation is an important part of the terrestrial ecosystem and plays a crucial role in sequestering carbon and mitigating climate change [10]. The vegetation ecosystem is found to be more vulnerable and sensitive to climate change than the other ecosystems. As a variable that reflects the efficiency of vegetation fixation and conversion of light

**Citation:** Wang, Z.; Wang, H.; Wang, T.; Wang, L.; Huang, X.; Zheng, K.; Liu, X. Effects of Environmental Factors on the Changes in MODIS NPP along DEM in Global Terrestrial Ecosystems over the Last Two Decades. *Remote Sens.* **2022**, *14*, 713. https://doi.org/10.3390/rs14030713

Academic Editors: Xuejia Wang, Tinghai Ou, Wenxin Zhang and Youhua Ran

Received: 16 December 2021 Accepted: 31 January 2022 Published: 2 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

energy, net primary productivity (NPP) is widely used in the monitoring of vegetation dynamics [11,12]. It is related to life activities such as vegetation growth, development, and reproduction, and it also provides an indispensable material basis for the life activities of other biological members in the entire ecosystem. Most of the studies focus on the spatiotemporal changes in NPP [13,14] and the impact of climate and human activities on NPP [15,16], modeling organic carbon storage with NPP as input data, and changes in carbon sources and sinks of the terrestrial ecosystem [12,17–19]. Studies on the changes in NPP with increasing elevation gradients (EG) are mainly conducted in local regions, and all come from instantaneous surveys [20–24]. Therefore, we still lack the knowledge and a full picture of the changes in NPP along EG at the global scale.

Typically, NPP declines with increasing elevation, which can largely be explained by the decreasing temperature as elevation increases [20,22]. However, the impact of temperature on the NPP of vegetation may be more significant in high-elevation areas with the intensification of global warming [25,26], because the warming rate in high-elevation areas is often greater than that in low-elevation areas [25,27]. Thus far, the warming rate has been intensifying at the global scale for decades [26,28–30], which is likely to cause a completely different spatial pattern of NPP on EG [31]. A few studies have focused on determining the spatial pattern of changes in vegetation greenness along EG in recent years [29,32,33], and their results indicated that the signal of the vegetation greenness increases with increasing elevation is found at the global scale and regional scale. However, vegetation greenness is not completely equal to vegetation productivity, and it does not directly participate in the process of the carbon cycle. This situation raises a scientific question as to whether the elevation pattern of NPP has changed under the influence of global warming.

Considering the importance of the effect of environmental variation on vegetation NPP, and the defects and gaps in current scientific research, we analyzed the inter-annual variation of MODIS NPP product (MOD17A2HGF, version 6.1) and environmental factors (air temperature, solar radiation, and soil moisture that derived from the reanalysis data of ERA5, MERRA2, and NCEP2). We calculated the changes in NPP, temperature, radiation, and soil moisture along EG by using the digital elevation model (DEM), which was also used to verify whether the elevation pattern of NPP has been altered under environmental variation. Furthermore, we analyzed the effect of environmental factors on the spatiotemporal changes and elevation distribution of the NPP from 2001 to 2020. We hope that the results of this study can provide references for the evaluation of terrestrial ecosystem carbon source and sink functions and the improvement and development of carbon cycle models. This study will contribute to our understanding of the impact of environmental factors on the elevational distribution of NPP in recent decades and will help us formulate strategies for mitigating climate change. We also expect that the outcomes of this study can provide references for the evaluation of terrestrial ecosystem carbon sink (source) functions and the improvement of current carbon cycle models.

#### **2. Materials and Methods**

*2.1. Datasets*

#### 2.1.1. NPP Data

The NPP data used in this study come from the MODIS (Moderate Resolution Imaging Spectroradiometer) remote sensing product of MOD17A3HGF Version 6.1 [34], which will be generated at the end of each year when the entire yearly 8-day MOD15A2H is available. MOD17A3HGF has two data fields named Npp\_500m and Npp\_QC\_500m, which represent NPP data and their quality control (QC). The poor-quality inputs were cleaned from the 8-day leaf area index and the fraction of photosynthetically active radiation based on the QC label for every pixel. The data type of Npp\_QC\_500m is an unsigned 1-byte integer (uint8) with a valid range from 0 to 100 (Units = 100%). Finally, we use the QC file to select pixels with good quality to participate in the analysis. The annual NPP, with a spatial resolution of 500 m, is derived from the sum of all 8 days of NPP in a certain year from 2001

to 2020. We further converted the coordinate system to World Geodetic System 1984 and resampled the spatial resolution to 0.008◦ (≈1 km) by using the bicubic method to match the resolution of the DEM data.

#### 2.1.2. DEM Data

The Shuttle Radar Topography Mission (SRTM) is a joint project between the National Geospatial-Intelligence Agency and the National Aeronautics and Space Administration. The SRTM30\_PLUS DEM dataset with a spatial resolution of 30-arc second (≈1 km) was developed by the Scripps Institute of Oceanography, University of California San Diego [35]. The coordinate system is the World Geodetic System 1984. The land data are based on the 1 km averages of topography derived from the United States Geological Survey (USGS) SRTM30 (30 arc-sec) grided DEM data product created with data from the NASA SRTM [36]. Global 30 Arc-Second Elevation (GTOPO30) data are used for high latitudes where SRTM data are not available. GTOPO30 is a global DEM with a horizontal grid spacing of 30 arc seconds (approximately 1 km) [37]. SRTM has been extensively used and provides a good representation of the topography [38], which performs better than Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) (the spatial resolution is 500 m) in terms of micro-topography, hydrologicnetwork and structural information characterization [39,40]. Considering the high quality of the SRTM data and the computational efficiency of the data, we chose SRTM30\_PLUS as the DEM data in this study.

#### 2.1.3. Reanalysis Data

We selected three sets of reanalysis data (Table 1) to study the spatiotemporal changes of environmental factors (temperature, radiation, and soil moisture), their distribution along EG, and the uncertainties of the correlation between environmental factors and NPP. The three sets of reanalysis data are the fifth generation of European ReAnalysis (ERA5), the recent Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA2), and National Centers for Environmental Prediction Reanalysis 2 (NCEP2). These reanalysis data can completely cover the data period of MODIS NPP and contain three kinds of environmental factor data at the same time. In data processing, we convert the unit of temperature from K to ◦C by subtracting 273.15, and the unit of radiation was converted from W/m<sup>2</sup> to MJ/m2 because the unit of absorbed photosynthetically active radiation is MJ/m<sup>2</sup> in the NPP simulation. Therefore, we can intuitively observe the magnitude of the spatial distribution of radiation. The soil moisture has different names in each reanalysis dataset, and we selected the variables with similar meanings and consistent units. Furthermore, we supplemented the monthly data with mean annual temperature (MAT), radiation (MAR), and soil moisture (MASM) raster data. Then, we extracted these environmental data from 2001 to 2020 for analysis. To keep the coordinate system of the reanalysis datasets consistent with the DEM, we set the resolution of the MERRA2 and NCEP2 data to 0.5◦ and 1.875◦, respectively. After that, we upscaled the DEM data to match them, and then performed the moving window operation. We explain the reason for this in Section 2.2. However, trend and correlation analysis were performed at the resolution of 0.008◦ to keep the data resolution consistent throughout the main text, and also to be consistent with the NPP data.

#### *2.2. Calculation of the Changes in NPP and Environmental Factors along EG*

The changes in NPP along EG (NPPEG) were calculated in three steps that are the same as those for the change of environmental factors along EG [29,32,33]. Step 1: We selected a 9 × 9 moving window to traverse the DEM and NPP (or environmental factors) raster data (Figure 1). The size of the moving window determines the amount of data used for analysis for each moving step. A previous study indicated that the difference caused by the window size is extremely small and recommended choosing a window size of 9 × 9 for calculation [29,33]. We also tested the window sizes of 5, 7, 9, and 11, and then expanded the size of the moving window to 19, 29, 39, and 49. We found that the changes in the areas of positive NPPEG were 0.0039% from the window size of 5 to 11, while the changes in the areas were reached 1% from the window size of 19 to 49. Therefore, we followed the recommendation of previous studies and chose a window size of 9 × 9 for the study. Step 2: The data extracted by the moving window have two dimensions. Thus, we need to reduce the dimensions and arrange them in pairs for linear regression analysis. A 9 × 9 moving window contains 81 data after dimensionality reduction. Step 3: We assign the regression coefficient obtained by linear regression analysis to the center pixel of the moving window and traverse the entire raster data accordingly to obtain the distribution map of the change of the variable along EG. The statistical significance of the slope is tested by a *t*-test. It should be noted that we set the NPP value of the ocean to no data. Therefore, when the moving window contains pixels of the ocean, the center pixel of the moving window will be filled with no data. In addition, we upscaled the DEM data to match the original resolution of environmental factors derived from three reanalysis datasets, such that we can avoid the regression analysis was established between the same environmental factors (coarser grid) and different DEM data (fine grid) in a moving window.

**Table 1.** Summary of reanalysis data products used in this study.


#### *2.3. Trend Analysis*

We use the ordinary least squares regression method [15,44–46] to determine the variation of NPP and environmental factors in time series:

$$\text{Slope} = \frac{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{i} \times (\text{VAR})\_{i} - (\sum\_{i=1}^{n} \mathbf{i}) (\sum\_{i=1}^{n} (\text{VAR})\_{i})}{\mathbf{n} \times \sum\_{i=1}^{n} \mathbf{i}^{2} - (\sum\_{i=1}^{n} \mathbf{i})^{2}} \text{ } \tag{1}$$

where VAR can be NPP, MAT, MAR, MASM, and their changes along EG; i is the sequence number of the year (from 2001 to 2020); n represents the total number of years. The significance of the trend was determined by the *T*-test.

#### *2.4. Correlation Analysis*

We use Pearson correlation analysis to calculate the correlation between NPP and environmental factors (MAT, MAR, and MASM). The correlation coefficient (r) of the two variables can be calculated by Equation (2).

$$\mathbf{r\_{xy}} = \frac{\sum\_{i=1}^{n} (\mathbf{x\_i} - \overline{\mathbf{x}})(\mathbf{y\_i} - \overline{\mathbf{y}})}{\sqrt{\sum\_{i=1}^{n} (\mathbf{x\_i} - \overline{\mathbf{x}})^2 \sum\_{i=1}^{n} (\mathbf{y\_i} - \overline{\mathbf{y}})^2}},\tag{2}$$

The correlation coefficient between NPP and MAT was taken as an example. Variables xi and yi denote the NPP and MAT in year i, respectively; x and y are the mean values of NPP and MAT from 2001 to 2020, respectively.

#### **3. Results**

#### *3.1. Trends of NPP and Environmental Factors*

NPP exhibited an increasing trend (1.75 gC/m2/yr) from 2001 to 2020 (Figure 2a), which was almost below the mean value for the first 10 years, and above the mean value for the next 10 years. The trends of the three kinds of reanalysis data showed that MAT was significantly increased. However, ERA5 was significantly decreased, while both MERRA2 and NCEP2 showed significant increasing trends. By contrast, the trends of MAR exhibited a decreasing trend, but only MERRA2 was statistically significant (Figure 2b–d). The spatial pattern of mean annual NPP is consistent with our perception that tropical forests have the highest NPP, while NPP is relatively low in alpine and arid regions (Figure 3a). The spatial pattern of the NPP trend showed that 70.63% of the global areas (16.21% were statistically significant) presented an increasing trend from 2001 to 2020, which were mainly found in central and western Canada, parts of central and northern China, central and southern South America, and central Africa (Figure 3b). By contrast, the regions with a decreasing trend of NPP occupied 29.37% of the global areas, and 3.02% of them were statistically significant, which were mainly distributed in northern South America.

The uncertainties of the trend of environmental factors were reflected in the spatial distribution (Figure 3.2). MAT showed a significant increasing trend in high−latitude regions of Asia and Europe, but opposite trends were found in central Africa (Figure 3.2a,d,g). The consistency of the positive trend of MAT accounted for 56.80% of the global area, whereas the areas with a negative trend accounted for 7.28%, mainly in northeastern North America, the Iranian plateau, and parts of the region across central India (Figure 3.2 j). Strong spatial heterogeneities are observed in the spatial pattern of MAR trends of the three sets of reanalysis data (Figure 3.2b,e,h). The inconsistent trends of MAR accounted for 61% of the global area. By contrast, the consistency of negative trends (27.55%) was greater than that of positive trends (11.65%). It is difficult to find areas where the MASM trend is consistent (Figure 3.2c,f,i), and MASM trends are inconsistent across most regions of the globe (73.84%) (Figure 3.2l). We found similar spatial patterns of MAT, MAR, and MASM (Figure 5). MAT is lower in high−elevation areas (Figure 5a,d,g). However, the cold climate of the Qinghai–Tibet Plateau is not reflected in NCEP2 MAT data (Figure 5g). Therefore, we speculated that the NCEP2 MAT probably fail to capture the temperature distribution. The higher uncertainty of MAT includes western North America, the Andes in South America, Greenland, the Mongolian Plateau, and the Qinghai–Tibetan Plateau (Figure 5j). MAR in high latitudes is lower than that in other regions (Figure 5b,e,h), which is highly uncertain in Greenland and the Sahara Desert (Figure 5k). MASM has relatively high values in the high latitudes of the northern hemisphere, the eastern United States, southeastern China, and tropical forests (Figure 5c,f,i). By contrast, it is lower in the western United States, the Mongolian plateau, the Sahara Desert, southern Africa, and Australia. The higher uncertainty of MASM is mainly found in the high latitudes of the northern hemisphere, eastern China, and tropical forests (Figure 5l).

**Figure 2.** Trend of the NPP (**a**), MAT, MAR, and MASM (**b**–**d**) anomaly from 2001 to 2020. The dashed straight lines denote the trendlines in subfigure (**a**). Asterisks (\*\*) denote that the slope is statistically significant at the 0.01 level.

**Figure 3.** Spatial pattern of global NPP (**a**) and its trend (**b**) from 2001 to 2020. The frequency of the uncertainty value is in the left of each subfigure. The regions with black dots in (**b**) indicate that the trend is statistically significant at the 0.05 level.

**Figure 4.** Spatial patterns of the trends of MAT (**a**,**d**,**g**), MAR (**b**,**e**,**h**), MASM (**c**, **f**,**i**) (derived from ERA5, MERRA2, and NCEP2), and their spatial consistencies (**j**,**k**,**l**) from 2001 to 2020. The frequency of the uncertainty value is in the left of subfigures (**j**,**k**,**l**). The regions with black dots indicate that the trend is statistically significant at the 0.05 level. "Pos" and "Neg" denote the regions with the positive and negative agreement, respectively. "Non" denotes that the regions have not reached an agreement (the same meaning thereafter).

#### *3.2. Elevational Distribution of NPP and Environmental Factors*

Changes in NPP and its slope along EG (NPPEG and NPPslope EG ) are shown in Figure 6. The positive NPPEG accounted for 59.98% of the global area, and 31.32% of them reached a significance level of 0.05. The positive NPPEG was scattered around the Great Lakes, the eastern foothills of the Andes, the eastern Brazilian plateau, the sub−Saharan African continent, the Indochina Peninsula, and eastern Australia. By contrast, the region with negative NPPEG accounted for 40.02%, and the significantly negative region occupied 11.18%. These regions are mainly found in northeastern North America, southern South America, and Central Asia (Figure 6a). The NPPslope EG showed strong spatial heterogeneity during the study period (Figure 6b), and there is no obvious spatial distribution feature exists. The region with positive NPPslope EG accounted for 53.88%. By contrast, 46.12% of the global area had negative NPPslope EG .

**Figure 5.** Spatial patterns of the MAT (**a**,**d**,**g**), MAR (**b**,**e**,**h**), MASM (**c**,**f**,**i**) (derived from ERA5, MERRA2, and NCEP2), and their spatial consistencies (**j**,**k**,**l**) from 2001 to 2020.

**Figure 6.** Spatial pattern of the changes in NPP (**a**) and its slope (**b**) along EG from 2001 to 2020. The value frequency is in the left of each subfigure. The regions with black dots indicate that the trend is statistically significant at the 0.05 level.

The spatial pattern of MAT, MAR, and MASM along EG (MATEG, MAREG, MASMEG) is shown in Figure 7. MATERA5 EG and MATERA2 EG follow the natural law that temperature decreases with the increase in elevation, and the area of the two reached 84.36% and 83.67% (Figure 7a,d), with the statistically significant area even reaching up to 78.57% and

74.35%, respectively. By contrast, only 32.90% of the area of MATNCEP2 EG conforms to the natural law, which indicates that the temperature data of NECP2 are likely to be wrong in elevational distribution (Figure 7g). Therefore, we excluded NCEP2 when calculating the spatial uncertainty of MAT. The results showed that the spatial consistency of MATERA5 EG and MATMERRA2 reached 78.84%. MAREG does not show obvious spatial distribution characteristics (Figure 7b,e,h), and the areas with positive MAREG (59.73%) were larger than those with negative MAREG (40.27%). The spatial consistency of the MAREG is 39.06%, which is mainly distributed in central and eastern North America, most of Europe, and central and western Russia. However, the inconsistent regions still accounted for most of the global area (Figure 7k). The areas with positive MASMEG have strong spatial heterogeneity, although it occupied (64.09 ± 7.02)% of the global area (Figure 7c,f,i), and the spatial inconsistency of MASMEG is up to 62.73% (Figure 7l).

**Figure 7.** Spatial patterns of MAT (**a**,**d**,**g**), MAR (**b**,**e**,**h**), MASM (**c**,**f**,**i**) (derived from ERA5, MERRA2, and NCEP2) along EG and their spatial consistencies (**j**,**k**,**l**) from 2001 to 2020. The frequency of the uncertainty value is in the left of subfigures (**j**,**k**,**l**).

We further investigated the spatial pattern of changes in the slope of MAT, MAR, and MASM along EG (MATslope EG , MARslope EG and MASMslope EG ) (Figure 8). The positive and negative MATslope EG accounted for 52.03% and 47.97% of the global area, respectively (Figure 8a,d,g). The area with consistent changes of MATslope EG accounts for 56.28% of the global area, and the remaining regions are highly uncertain (Figure 8j). The areas with positive MARslope EG occupied (45.50 ± 3.00)% of the global area, which is less than that of

negative MARslope EG (54.50 <sup>±</sup> 3.00)% (Figure 8b,e,h). Anyway, MARslope EG is highly uncertain in most of the global area (68.11%) (Figure 8k). The areas of positive MASMslope EG (51.49%) are slightly more than that of negative MASMslope EG (48.51%) (Figure 8c,f,i). The inconsistent region covered up to 72.09% of the global area (Figure 8l).

**Figure 8.** Spatial patterns of the slope of MAT (**a**,**d**,**g**), MAR (**b**,**e**,**h**), and MASM (**c**,**f**,**i**) (derived from ERA5, MERRA2, and NCEP2) along EG and their spatial consistencies (**j**,**k**,**l**) from 2001 to 2020. The frequency of the uncertainty value is in the left of subfigures (**j**,**k**,**l**).

#### *3.3. Effect of Environmental Factors on the Elevational Distribution of NPP*

Figure 9 illustrates spatial patterns of the dominating environmental factor on the changes in NPP from 2001 to 2020. ERA5 data showed that MAT and MAR seem to be the dominating factors on the changes in NPP in the north of 30◦N, central and northern parts of South America, central Africa, and Southeast Asia, whereas MASM mainly affects the changes in NPP in the central and eastern United States, the eastern part of the Brazilian plateau, southern Africa, and Australia (Figure 9a). By contrast, the effects of environmental factors of MERRA2 on the changes in NPP have a clear spatial distribution pattern. MAT mainly affects the changes in NPP in the Qinghai–Tibet Plateau and parts of central South America. MAR is the dominant factor in the high latitudes of the northern hemisphere, southern South America, and Southeast Asia. MASM has a wider range of influence, including southern North America, Eurasia from 30◦N to 60◦N, the Indian peninsula, most of Africa, and Australia (Figure 9b). The MAT and MAR of NCEP2 are the dominant factors in the changes in NPP in the north of 30◦N. MAR is also the main environmental factor that affects NPP changes in northwestern South America, central Africa, and Southeast Asia. The areas where MASM presented a dominant environmental factor include most of the United States, northern South America and most of Brazil, the southern edge of the Sahara

Desert and southern Africa, most of Eurasia from 30◦N to 60◦N, the Indian Peninsula, and Australia (Figure 9c). MASM has the highest spatial consistency (25.15%). These regions are mainly distributed in the United States, eastern and southern South America, southern Africa, Eastern Europe, Central Asia, parts of East Asia, India, and Australia, followed by MAR (9.74%) and MAT (7.70%), which can be found in high latitudes in the northern hemisphere (Figure 9d).

**Figure 9.** Spatial patterns of the effect of MAT, MAR, and MASM derived from ERA5 (**a**), MERRA2 (**b**), and NCEP2 (**c**) on the changes in NPP and their spatial consistencies (**d**) from 2001 to 2020.

We found that NPP has an obvious turning point at the elevation of 3050 m and divided the elevation into high and low gradients based on the turning point (Figure A1 in Appendix A). Then, we explored the impact of environmental factors on NPP in high− and low−elevation areas (Figure 10a). The three kinds of reanalysis datasets showed that MAT still maintained the highest R<sup>2</sup> in the three environmental factors. We further fitted the regression equation curve and uncertainty range of NPP and environmental factors (Figure 10). The effect of MAT on NPP is linear and non−linear at low and high elevations, respectively (Figure 10a). The R<sup>2</sup> between MAT and NPP exceeds 0.9 and is statistically significant at the 0.01 level. In general, the uncertainty of the effect of MAT on NPP gradually decreases as the temperature increases, and it decreases more at low elevations. The effect of MAR on NPP is nonlinear at both low and high elevations (Figure 10b). NPP increases as MAR increases at low elevations, whereas opposite trends were observed at high elevations. The uncertainty of the impact of MAR on NPP gradually decreases with the increase in MAR at low elevations. By contrast, the uncertainty does not decrease significantly as the MAR decreases at high elevations (Figure 10c). We found a linear and nonlinear effect of MASM on NPP at low and high elevations, respectively. NPP increases as MASM increases at low elevations. However, the upper and lower limits of uncertainty trend toward the opposite direction. The same situation occurs at high elevations that larger uncertainty of the effect of MASM on NPP, and MASM remains stable with the increase in NPP (Figure 10d). The R2 of all fitted equations is statistically significant at both high and low elevations.

#### **4. Discussion**

#### *4.1. Elevational Patterns of Environmental Factors and NPP*

In this article, we use the moving window method to explore the changes in NPP and environmental factors along EG at the global scale. Unlike field surveys, this method is based on remote sensing and DEM data, which can quickly and economically monitor the changes of variables along EG. In addition, the moving window method used in this article is intended to investigate the change of NPP and environmental factors at the relative EG within the moving window, which is different from the studies on the absolute EG conducted through a field survey. The accuracy of this method can be indirectly proved by the natural law that the temperature decreases with the increase in elevation. Nearly 80% of the regions showed a decrease in temperature with increasing elevation, and more than 70% of the regions are statistically significant. This result means that the moving window method can effectively monitor the change of the variable in the EG and has high accuracy. Undeniably some regions (although very few) have a local micro-climate, the temperature of which will not decrease with the increase in elevation. We can also evaluate the data quality of temperature based on the natural law of temperature–elevation and the moving window method. For example, NCEP2 fails to capture the elevation distribution of MAT because it is contrary to the natural law that the temperature decreases with the increase in elevation. Unfortunately, we are still unable to assess the data quality of radiation and soil moisture because the spatial distributions of the reanalysis data and their trends are highly uncertain, and there is no corresponding natural law to follow. Typically, NPP decreased with increasing elevation (negative NPPEG) because of the limitation of low temperature. However, we found that the areas with positive NPPEG accounted for 59.98% of the global area, which means the elevation pattern of NPP in our perception has changed. The spatial pattern of NPP slope demonstrated that, compared with NPP in low-elevation areas, NPP in high-elevation areas has a higher increase rate, and this phenomenon has become more common worldwide.

#### *4.2. Uncertainty of Environmental Factors and Their Effect on NPP*

We found that MAR and MASM had larger uncertainties than MAT in the interannual variations and the spatial distribution along EG. Such a large uncertainty makes it difficult to assess their effect on the changes in NPP and further forms a lower spatial

correlation consistency. Factors affecting MAR mainly include cloud height, thickness, shape, and aerosols. [47]. The factors that affect soil moisture include soil texture [48], the data assimilation process [49], the sample size of observation data [50], and land use type [51].

We found that the effect of radiation on NPP is nonlinear regardless of whether it is in high- or low-elevation regions. The positive correlation between NPP and radiation in low-elevation is likely to be that, as the elevation increases, temperature and water inhibit the photosynthesis of plants, thereby reducing the use of photosynthetically active radiation. Moreover, the radiation available for plants decreases, thereby forming a nonlinear positive correlation between NPP and radiation. However, the non-linear negative correlation between radiation and NPP in high-elevation regions suggests that radiation is not the main driving force for changes in NPP. We suppose that radiation cannot be fully utilized by plants in high-elevation areas, and temperature is still the main environmental factor that determines changes in NPP along EG. The effect of soil moisture on NPP shows a non-linear positive correlation in low elevation areas, but it is highly uncertain because of the large differences in the trend of soil moisture among ERA5, NCEP2, and MERRA2. Anyway, the positive correlation between NPP and soil moisture has also been confirmed in the arid region [52], whereas soil moisture in high-elevation areas is higher than that in low-elevation areas, it has not caused an increase in NPP. Therefore, we hold the view that soil moisture determines the lower limit of NPP, and the upper limit of NPP is determined by both soil moisture and temperature.

Overall, the temperature is still the dominating climate factor that determines the spatial patterns of NPP along EG, and it is generally positively correlated with NPP whether in high- or low-elevation areas. However, NPP is more sensitive to temperature in low-elevation areas, which indicates that temperatures are more conducive to vegetation photosynthesis in this region. By contrast, the temperature gradually decreases with the increase in high-elevation, and the photosynthesis of vegetation is also limited by low temperature. Thus, NPP appears to rise first and gradually stabilizes. The effect of temperature on NPP along EG has been confirmed in regional-scale studies [53,54].

#### *4.3. Implications for Carbon Cycle*

NPP is the source of energy in the carbon cycle of the terrestrial ecosystem. Our research shows that nearly 60% of the global area exhibits an increase in NPP with increasing elevation, which means that vegetation in high-elevation areas plays an increasingly prominent role in absorbing atmospheric CO2 and mitigating climate change. However, this increase in NPP is caused by global climate change, especially the increase in temperature in high mountain areas [33], which improves the photosynthesis capacity of plants. However, when the temperature exceeds the optimal temperature of the plants, it will inhibit the photosynthesis of plants, and even cause the death of local species because they cannot adapt to the rapid warming. This kind of research on vegetation degradation caused by warming has been widely reported [12]. What we want to emphasize is that this phenomenon of increasing NPP with elevation may be beneficial to CO2 fixation in the short term, and an uncontrolled increase in temperature will inevitably lead to vegetation degradation and even ecosystem collapse. The increase in temperature will also cause the increase in plant autotrophic respiration and soil heterotrophic respiration, and the CO2 produced by the respiration process is directly discharged out of the vegetation–soil system. With the high degree of uncertainty in soil respiration, much uncertainty exists in the quantitative evaluation of the carbon source and sink functions of the ecosystem.

Environmental factors have a strong influence on the terrestrial carbon cycle. Temperature is the basic input data to establish its impact on the soil carbon cycle. In this study, the difference in the spatial distribution of ERA5 and MERRA2 temperature is very small. However, the algorithm for the effect of temperature on the decomposition of soil organic carbon (shown as a nonlinear positive correlation effect) has slight differences. The differences in the algorithms are reflected in space, mainly located in high latitude

areas and the Qinghai–Tibet Plateau with great uncertainty. This condition means that the uncertainty of the temperature algorithm in the alpine region will greatly affect the decomposition of soil organic carbon in the carbon cycle, which in turn affects the spatial distribution of soil organic carbon storage in the alpine region, leading to distortions in the evaluation of the carbon source and sink functions of these regions.

Numerous studies have shown that soil moisture is crucial to the carbon cycle process of terrestrial ecosystems [55–57]. However, the three sets of reanalysis data show that the spatial distribution of soil moisture has great uncertainty. Similarly, there are large discrepancies in the algorithm for the effect of soil moisture on the decomposition of soil organic carbon. The spatial uncertainty of soil moisture involved tropical rain forests and the Sahara Desert [17]. The fundamental reason is that the understanding of the effect of soil moisture on the decomposition of soil organic carbon varies, especially as to whether highhumidity environments will inhibit the decomposition of soil organic matter. An anaerobic environment created by high humidity will inhibit the decomposition of soil organic matter in the carbon cycle model, which will lead to the accumulation of soil organic carbon (high soil organic carbon), and vice versa, leading to the decomposition of soil organic carbon (low soil organic carbon). The uncertainty of the original input data of soil moisture and the algorithm makes it difficult for us to evaluate its effect on the carbon cycle. Considering the importance of soil moisture in the carbon cycle of terrestrial ecosystems, we propose that future studies need to further strengthen the observation of soil moisture and improve its simulation accuracy to more accurately assess the carbon source and sink functions of terrestrial ecosystems.

#### **5. Conclusions**

This study analyzed the spatiotemporal changes of the MODIS NPP product and environmental factors (temperature, radiation, and soil moisture data derived from the reanalysis data ERA5, MERRA2, and NCEP2) and their distribution along EG. We also identified the spatial uncertainty of environmental factors and their effects on the elevation distribution of NPP. We found that nearly 60% of the global area presented an increase in NPP with increasing elevation, which implied that the elevation pattern of NPP has changed, and the carbon sequestration capacity of vegetation is increasing elevation. However, soil respiration is likely to increase as well. Quantitatively evaluating the carbon sink (source) function of vegetation remains to be further studied in high-elevation regions. The temperature of NCEP2 failed to capture the alpine environment of the Qinghai–Tibet Plateau, and it does not clearly show the natural law that the temperature decreases with the increase in elevation. Soil moisture has the largest areas of spatial consistency in affecting the spatiotemporal changes in NPP among the three environmental factors. However, its spatial pattern and variation are the most uncertain among the three environmental factors, even though it is essential to the carbon cycle of terrestrial ecosystems. NPP has obvious elevation differentiation with an elevation of 3060 m as the demarcation point, which divides the elevation into low and high. MAT is the main driving force that affects the elevation distribution of NPP, with its effect on NPP exhibiting a significant linear and nonlinear positive correlation at low and high elevations, respectively. The results of this study are expected to contribute to our understanding of the changes in NPP along EG and provide references for the development of terrestrial ecosystem carbon cycle models.

**Author Contributions:** Conceptualization, Z.W. and X.H.; data curation, Z.W., K.Z. and X.L.; Formal analysis, Z.W.; investigation, Z.W., H.W. and T.W.; methodology, Z.W. and L.W.; writing—original draft, Z.W.; writing—review and editing, Z.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant No. 32160278); the Open Project of State Key Laboratory of Plateau Ecology and Agriculture, Qinghai University (Grant No. 2020-ZZ-14, 2021-KF-08); the National Natural Science Foundation of China

(Grant No. NSFC32001188); the Natural Science Foundation of Qinghai Province of China (Grant No. 2021-ZJ-973Q); the National Natural Science Foundation of China (Grant No. 42067070).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are openly available in (NASA EOSDIS Land Processes DAAC) at (doi.org/10.5067/MODIS/MOD17A3HGF.061), reference number [33].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Changes in NPP and environmental factors along EG from 2001 to 2020. The elevation of 3060 m is used to distinguish high and low elevation (the vertical black solid line in Figure (**a**–**j**)). The dashed line represents the average value of the variables at high and low elevations. The shaded areas with different colors in each figure represent ±1 SD.

#### **References**

