Regional Contribution and Attribution of the Interannual Variation of Net Primary Production in the Yellow River Basin, China

Abstract

Net primary production (NPP) serves as a crucial indicator of the ecosystem’s capacity to capture atmospheric CO₂. Gaining insights into the dynamics of NPP and its driving mechanisms is pivotal for optimizing ecosystem carbon sink resource management. Since the implementation of the Grain-for-Green Program (GFGP) in 1999, the Yellow River Basin (YRB) has been one of the most significant areas for ecological restoration in China. However, our knowledge regarding the interannual variability (IAV) of NPP and the underlying driving forces in this region remains incomplete. In this study, we utilized a light use efficiency model to assess the spatiotemporal dynamics, IAV, and driving factors of NPP in the YRB during the period from 1999 to 2018. Our findings revealed that the average annual NPP in the YRB approximated 189.81 Tg C. Over the study duration, NPP significantly increased in 79.63% of the basin with an overall increasing rate of 6.76 g C m⁻² yr⁻¹. The most prominent increase was observed in the key GFGP implementation area, predominantly in the semi-humid region. Notably, the middle altitude region (1–1.5 km), semi-humid region, and grassland emerged as the primary contributors to the basin’s total vegetation carbon sequestration. However, it is worth emphasizing that there was substantial IAV in the temporal trends of NPP, with the semi-humid region being the most influential contributor (62.66%) to the overall NPP IAV in the YRB. Further analysis of the driving mechanisms unveiled precipitation as the primary driver of NPP IAV in the YRB with a contribution of 62.9%, followed by temperature (23.07%) and radiation (14.03%). Overall, this study deepened our understanding of the IAV and driving mechanisms of NPP in the YRB under ecological restoration, and provided scientific support for optimizing the management of regional carbon sequestration resources.

1. Introduction

The quantity of tree biomass per unit area of land constitutes the primary inventory data needed to understand the flow of materials and water through forest ecosystems [1,2]. Net primary production (NPP) is the net amount of carbon absorbed from the atmosphere via photosynthesis, and it is an important component of the carbon balance of all terrestrial ecosystems [3,4,5,6]. NPP is crucial for ecosystem services, agriculture, and forestry [7,8,9], as well as global carbon cycling, which in turn affects water cycling and climate. NPP can intuitively reflect the production capacity of surface vegetation under natural environmental conditions, and is often used as an important factor for determining carbon sources/sinks and regional ecological carrying capacity [10]. Therefore, chronosequence assessment of NPP and its driving forces would help to predict the vegetation development and dynamics of the species [11,12,13,14].

Although many studies have evaluated the magnitudes and dynamics of regional and even global NPP based on various approaches [15,16,17,18,19], large uncertainties and conflicts still exist among different estimates. Taking the Multi-scale Synthesis and Terrestrial Model Intercomparison Project (MsTMIP) [20] as an example, estimates of the NPP in China using 13 terrestrial biosphere models (TBMs) forced by the same drivers in this project exhibit substantial inconsistency in their magnitudes, spatiotemporal patterns, and dynamics [21]. The largest NPP in China estimated using the SiB model (4.4 Pg C) is almost twice of that (2.3 Pg C) simulated by the ISAM model [21]. Light use efficiency models, such as the Carnegie–Ames–Stanford approach (CASA) [22], utilize remote-sensing-based vegetation indices and climate–hydrological data to simulate vegetation productivity. Compared to process-oriented models, their principles are relatively straightforward and involve fewer parameters; therefore, light use efficiency models do not require complex model parameter calibration [23]. The input data for these models are easily obtainable, making them well suited to large-scale applications [24]. Previous applications have indicated that the simulation results from light use efficiency models align reasonably well with real-world conditions [23,24,25,26,27]. In addition, most studies have mainly focused on the trends of NPP, and our attention on the climate-driven NPP interannual variability (IAV) is still limited. The IAV of NPP reflects the annual variation in plant net photosynthesis, which is a key indicator that affects the variation of the atmospheric CO₂ concentration [28,29,30] and ecosystem supply function [31,32]. Although the IAV of terrestrial carbon sequestration flux has attracted increasing interest, debates still exist when involving the regional contributions and driving mechanisms [21,33,34,35]. Understanding the IAV of land NPP is crucial for estimating the future evolution of atmospheric CO₂ as carbon–climate feedbacks emerge [28,29,30,36]. At a global scale, tropical [33,34,37,38] and semiarid ecosystems [33,39,40] play an important role in influencing the IAV of the global land carbon sink. The major climate factors driving the annual land carbon sink changes also exhibit different patterns between regional and global scales [21,41,42,43]. Therefore, the primary climate factors driving NPP IAV may differ at different scales. Deciphering the regional contributions and the driving forces of the NPP IAV can help us identify the IAV hot spots of the ecosystem carbon sink function and provide effective regulatory measures under the urgent need of carbon neutrality.

The Yellow River is the second longest river in China and the Yellow River Basin (YRB) is one of the most sensitive areas to climate change due to severe water scarcity and intense human activities in China [6,44,45,46,47]. In recent decades, under the combined influence of climate change and human activities, frequent dry riverbeds and sharp decreases in surface runoff have occurred in the downstream area of the YRB [48]. This has led to a growing conflict between water supply and demand, significant changes in the ecological environment, and to a large extent, affected biodiversity and ecological balance [49,50,51,52]. However, in recent decades, a series of ecological projects (e.g., the Grain-for-Green Program is one of the most successful programs) has substantially improved the ecosystem situation in the YRB and Loess Plateau [6,53]. In this context, the academic community has shown increasing interest in the vegetation status and functionality, especially the carbon sequestration function, of the YRB [54]. Many studies have been conducted to investigate the spatiotemporal variations and driving forces of NPP in the YRB [55,56,57]. However, compared to changes in NPP, its IAV may be even more crucial, as it reflects the stability and resilience of the regional vegetation carbon sequestration function [21,35]. This index can better assess the long-term value and carbon sequestration benefits of ecological restoration projects. However, our understanding of the spatiotemporal characteristics of NPP IAV in the YRB and its climate-driven mechanisms is still limited. Under the urgent need of ecological protection and carbon neutrality in the YRB, detecting the spatiotemporal dynamics of the NPP IAV in the YRB and its responses to climate change can provide scientific references for evaluating the ecological benefits of ecological restoration programs and deepening our understanding of the resistance of artificially restored ecosystems to climate disturbances.

This study utilized a light use efficiency model (i.e., CASA), vegetation index, and climate–hydrological data to estimate the spatiotemporal dynamics of NPP and its IAV in the YRB from 1999 to 2018 and evaluated the accuracy of the model by using multiple data sources. We quantified the regional contributions to the holistic NPP IAV of the YRB and identified the IAV hotspots. We extended our study by deciphering the driving mechanisms of NPP IAV in different regions. The goal of this research is to provide scientific support for optimizing future ecological restoration measures and the management of regional terrestrial carbon sink resources.

2. Materials and Methods

2.1. Study Area

The YRB is located in North China (96°E~119°E, 32°N~42°N) and covers an area of 79.5 million hectares (Figure 1). Most parts of the YRB are located in the Tibetan Plateau, arid, semi-arid, and semi-humid regions. The annual precipitation varies from 110 mm to 912 mm and the annual mean temperature is between −19 °C and 15 °C [58]. The YRB is characterized by a temperate continental monsoon climate. The distribution of precipitation shows large spatial variations with a decreasing gradient from south to north. High-temperature regions are mainly located in the east of the basin and the low-temperature regions are mainly located in the western part of the basin (mostly in the Tibetan Plateau). The elevation in the YRB ranges from −1 m to 6142 m with an increasing gradient from the east to west. The YRB is one of the most sensitive areas to climate change in China due to severe water scarcity and intense human activities. During the last century, many large-scale ecological restoration measures have be implemented in the YRB to control soil erosion and ecological degradation [53]. The Grain-for-Green Program is one of the most successful ecological restoration projects in China and has significantly increased the vegetation coverage in this region [59,60].

2.2. Data Sources

In this study, we utilized several spatiotemporal datasets, including meteorological and hydrological data (monthly precipitation, monthly average temperature, monthly actual evapotranspiration, and monthly potential evapotranspiration), monthly total solar radiation, a Normalized Difference Vegetation Index (NDVI), a Digital Elevation Model (DEM), and a vegetation type map. The monthly precipitation and monthly average temperature with a spatial resolution of 1 km × 1 km were collected from an open meteorological database [58]. The monthly total solar radiation, monthly actual evapotranspiration (ET), and potential evapotranspiration (PET) with a spatial resolution of 4 km × 4 km were extracted from the TerraClimate dataset [61]. The high-spatial-resolution (1 km × 1 km) NDVI data with a 10-day time frequency were collected from the Satellite Pour l’Observation de la Terre (SPOT) VEGETATION Collection 3 (SPOT/VGT-C3 NDVI) dataset [62,63,64]. We generated the monthly NDVI using the monthly maximum compositing method. The DEM data with a 1 km × 1 km spatial resolution were extracted from the Global Land One-km Base Elevation Project [65]. The vegetation type map and land use and land cover maps with a spatial resolution of 1 km × 1 km from 2000 to 2018 were collected from the Resource and Environment Data Cloud Platform (REDCP, www.resdc.cn, accessed on 21 May 2020). The climate zone boundary was provided by Li et al. [21]. Data used to calculate NPP in this study are listed in

Table 1. Data used to calculate NPP in this study.

Data	Spatial Resolution	Time Span	Resources
Precipitation	1 km × 1 km	1999~2018	Peng et al. [58]
Temperature
Actual Evapotranspiration	4 km × 4 km	1999~2018	TerraClimate [61]
Potential Evapotranspiration
Solar Radiation
NDVI	1 km × 1 km	1999~2018	Baret et al. [62]
Vegetation Type	1 km × 1 km	/
Land Use and Land Cover	1 km × 1 km	2000~2018	www.resdc.cn, accessed on 21 May 2020

.

2.3. NPP Estimation

This study used the light use efficiency model, i.e., the Carnegie–Ames–Stanford approach (CASA) [22], to estimate the NPP of the YRB. The specific technical process is as follows:

The vegetation NPP can be represented using two factors: absorbed photosynthetically active radiation (APAR) and actual light use efficiency (ε):

$$\mathrm{N}\mathrm{P}\mathrm{P}\left(x,t\right)=\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}\left(x,t\right)\times \epsilon \left(x,t\right)$$

(1)

In Equation (1), APAR(x, t) represents the photosynthetically active radiation absorbed by pixel x in month t (unit: MJ m⁻²), and ε(x, t) represents the actual light use efficiency of pixel x in month t (unit: g C MJ⁻¹).

Potter et al. [22] indicated that vegetation has maximum light use efficiency under ideal conditions, but in reality, the maximum light use efficiency is mainly affected by temperature and water availability, which can be calculated using Equation (2):

$$\epsilon (x,t)=f_1(x, t)\times f_2(x, t)\times W(x, t)\times {\epsilon }_{max}$$

(2)

where f₁(x, t) and f₂(x, t) represent the low-temperature and high-temperature stress on light use efficiency (unitless), W(x, t) is the water stress coefficient (unitless) reflecting the impact of water conditions, and ε_max is the maximum light use efficiency under ideal conditions (unit: g C MJ⁻¹). The formulas of f₁, f₂, and W can be obtained from Zhu et al. [66].

Combining Equations (1) and (2), Equation (3) can be obtained:

$$\mathrm{N}\mathrm{P}\mathrm{P}\left(x,t\right)=\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}\left(x,t\right)\times f_1\left(x, t\right)\times f_2\left(x, t\right)\times W\left(x, t\right)\times {\epsilon }_{max}$$

(3)

For a given pixel, when the NPP, APAR, f₁, f₂, and W are known, ε_max on that pixel can be calculated using Equation (3). Then, the calculated ε_max values of numerous pixels are classified by vegetation type, and the maximum light use efficiency of each vegetation type is simulated based on the principle of minimum error.

For a certain vegetation type, the error between its observed NPP and simulated NPP can be expressed by a function:

$$EE\left(x\right)={\sum }_{i=1}^j{(m_i-n_{i}x)}^2, x\in \left[l,u\right]$$

(4)

In the equation, i represents the number of samples for a certain vegetation type. j is the maximum number of samples for a certain vegetation type. m is the observed NPP data. n is the product of the APAR, temperature, and water stress factors, which means APAR × f₁ × f₂ × W. x is the maximum light use efficiency for a certain vegetation type. l and u are the minimum and maximum values of the maximum light use efficiency for each vegetation type. Then, we obtained Equation (5):

$$EE\left(x\right)={\sum }_{i=1}^j{n_i}^2{x}^2-2{\sum }_{i=1}^j{m}_i{n}_{i}x+{\sum }_{i=1}^j{m}_i^2, x\in \left[l,u\right]$$

(5)

which is a quadratic equation with an upward-facing parabola. There must be a minimum value in the interval [l, u], and the corresponding x is the simulated maximum light use efficiency ε_max for a certain vegetation type.

Estimating the fraction of absorbed photosynthetically active radiation (APAR) by vegetation in the photosynthetically active radiation (PAR, 0.4~0.7 μm) is achieved based on the reflection characteristics of the vegetation in the red and near-infrared bands in the remote sensing data. The absorbed PAR by vegetation depends on the total solar radiation and the characteristics of the plant itself, and can be calculated using Equation (6):

$$\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\mathrm{S}\mathrm{O}\mathrm{L}(x,t)\times \mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)\times 0.5$$

(6)

In the equation, SOL(x, t) represents the total solar radiation (MJ·m⁻²) at pixel x in month t; FPAR(x, t) is the fraction of absorbed PAR by the vegetation layer (unitless); and the constant 0.5 represents the proportion of photosynthetically active radiation (wavelength of 0.38~0.71 μm) that the vegetation can utilize out of the total solar radiation.

Within a certain range, there is a linear relationship between the FPAR and NDVI, which can be determined based on the maximum and minimum NDVI values and the corresponding FPAR maximum and minimum values for a certain vegetation type:

$${\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}\left(x,t\right)=\frac{(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(x,t\right)-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{i,min})}{({\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{i,max}-{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{i,min})}\times \left({\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{max}-{\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{min}\right)+{\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{min}$$

(7)

In the equation, NDVI _i,max and NDVI _i,min correspond to the maximum and minimum NDVI values for the ith vegetation type.

Further research has shown that there is also a good linear relationship between FPAR and the Ratio Vegetation Index (SR), which can be expressed as follows:

$${\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{\mathrm{S}\mathrm{R}}\left(x,t\right)=\frac{(SR\left(x,t\right)-{SR}_{i,min})}{({SR}_{i,max}-{SR}_{i,min})}\times \left({\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{max}-{\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{min}\right)+{\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{min}$$

(8)

In the equation, FPAR_min and FPAR_max are independent of the vegetation type, and are 0.001 and 0.95; the SR is determined by Equation (9), where SR_i,max and SR_i,min correspond to the maximum and minimum NDVI_i,max and NDVI_i,min for the ith vegetation type.

$$SR\left(x,t\right)=\frac{1+\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(x,t\right)}{1-\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\left(x,t\right)}$$

(9)

Comparing the estimated results of FPAR-NDVI and FPAR-SR, it was found that the FPAR estimated using the NDVI is higher than the observed value, while the FPAR estimated using the SR is lower than the observed value, but the error is smaller than that of the direct estimation using the NDVI. Taking this into account, based on the flux observation and remote sensing monitoring data, this study optimized the model parameters and finally combined Equations (7) and (8), taking their average value as the estimated value of FPAR:

$$\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(x,t)=\alpha \times {\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}+(1-\alpha )\times {\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}}_{\mathrm{S}\mathrm{R}}$$

(10)

In the equation, FPAR_NDVI is the estimated value by Equation (7), FPAR_SR is the estimated value by Equation (8), and α is the adjustment coefficient between the estimated value and the observed value.

2.4. Model Accuracy Validation

We assessed the accuracy of our results using MODIS NPP data, in situ observations, and simulations of 14 terrestrial biosphere models (TBMs). The MOD17A3 product has undergone validation for NPP studies at both regional and global scales [16,67]. Numerous scholars regard these data as the standard for assessing the reliability of their results [21,68,69]. We conducted pixel-scale correlation analysis between our estimates and the MODIS NPP product. We also evaluated the performance of our estimates by comparing them with in situ observations [70]. In addition, we further compared our results with simulations of these TBMs. We only evaluated the overall magnitude and temporal trend consistencies of the NPP in the YRB between our result and the TBMs simulations due to the coarser spatial resolutions of the process-oriented models. The CLM5.0 simulations (1999~2018) originate from Friedlingstein et al. [71], while other TBM simulations (1999~2010) are derived from the MsTMIP [20]. We used four criteria to evaluate the quality of our estimates, including the coefficient of determination (R²), Pearson correlation coefficient (r), root mean square error (RMSE), and percent bias (PB).

$$r=\frac{{\sum }_{i=1}^n\left(Y_{i,CASA}-{\overline{Y}}_{CASA}\right)\left(Y_{i,com}-{\overline{Y}}_{com}\right)}{\sqrt{{\sum }_{i=1}^n{\left(Y_{i,CASA}-{\overline{Y}}_{CASA}\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(Y_{i,com}-{\overline{Y}}_{com}\right)}^2}}$$

(11)

$$R^2=1-\frac{{\sum }_{i=1}^n{(Y_{i,com}-Y_{i,CASA})}^2}{{\sum }_{i=1}^n{(Y_{i,com}-{\overline{Y}}_{CASA})}^2}$$

(12)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(Y_{i,CASA}-Y_{i,com}\right)}^2}$$

(13)

$$PB=\frac{1}{n}{\sum }_{i=1}^n\left(\frac{Y_{i,CASA}-Y_{i,com}}{Y_{i,com}}\times 100\%\right)$$

(14)

where Y_i,CASA and Y_i,com are the CASA-based NPP and comparative data, respectively. ${\overline{Y}}_{CASA}$ and ${\overline{Y}}_{com}$ are the mean values of the CASA-based NPP and comparative data, respectively.

2.5. Spatiotemporal Variation Analysis

Univariate linear regression analysis can simulate the change trend of each grid, and comprehensively characterize the regional pattern evolution of a certain time series using the spatial change characteristics of individual pixels at different periods. The advantage of this method is that it can eliminate the influence of abnormal factors on the trend analysis of NPP and truly reflect the evolutionary trend of NPP over a long time series. In this study, the annual total NPP was analyzed for its spatial evolution trend using the univariate regression trend analysis method on a per-pixel basis. The NPP of the YRB was regressed against time for each pixel. The slope of the regression equation for each pixel reflects the trend of NPP evolution, and a positive slope indicates an overall upward trend, while a negative slope indicates a downward trend. The magnitude of the slope reflects the degree of increase or decrease in the pixel, and the larger the absolute value of the slope, the more drastic the change. The slope was estimated as:

$$Slope=\frac{n\times {\sum }_{i=1}^n\left(i\times {NPP}_i\right)-\left({\sum }_{i=1}^{n}i\right)\left({\sum }_{i=1}^n{NPP}_i\right)}{n\times {\sum }_{i=1}^n{i}^2-{\left({\sum }_{i=1}^{n}i\right)}^2}$$

(15)

where Slope represents the slope of the linear fitting equation, i represents the current year, n represents the study period, and NPP_i represents the NPP for the ith year. The significance of the trend was tested using the Mann–Kendall (MK) detection method [72,73].

2.6. Quantification of Regional Contributions to NPP IAV

This study used a pixel-based partitioning contribution calculation method developed by Ahlstrom et al. [33] to evaluate the contribution of different regions to the overall IAV of NPP. The equation is shown as:

$$f_j=\frac{\sum _{t=1}^N\frac{x_{jt}\left|X_t\right|}{X_t}}{{\sum }_{t=1}^N\left|X_t\right|}$$

(16)

where x_jt represents the IAV of the NPP for pixel j in year t, X_t represents that of the entire study area, and f_j represents the contribution rate of pixel j to the study area. Pixels with larger absolute contribution rates contribute more to the whole study area, while pixels with smaller absolute contribution rates have less influence on the IAV of the study area. Positive values indicate that the pixel has a positive promotion effect, while negative values indicate that the pixel has a negative inhibitory contribution for the entire study area. When we calculate according to different partitions, the sum of the contribution rates of all pixels in each partition indicates the total contribution rate of that partition to the study area.

2.7. Attribution Analysis

In order to assess the correlation between climate change factors and NPP, we used the Pearson correlation coefficient to evaluate the correlation between the time series of NPP and its influencing factors. Additionally, we used the Lindeman–Merenda–Gold (LMG) model [74] to quantitatively evaluate the relative importance of different drivers.

$$\mathrm{L}\mathrm{M}\mathrm{G}\left(x_k\right)=\frac{1}{p!}{\sum }_{S\subseteq \{x_{1,\dots ,x_p}\}\backslash \{x_k\}}n\left(S\right)!\left(p-n\left(S\right)-1\right)!seq{R}^2\left(\right\{x_k\left\}\right|S)$$

(17)

where x_k is the explanatory variable and S means a set of variables.

$$seq{R}^2\left(\left\{x_k\right\}|S_k\left(r\right)\right)=R^2\left(\left\{x_k\right\}\cup S_k\left(r\right)\right)-R^2\left(S_k\left(r\right)\right)$$

(18)

where r represents the order of the regressors (r₁, …, r_p), S_k(r) denotes the set of regressors entered into the model before regressor x_k in the order r. This evaluation method can effectively avoid the order effects of regression variables and accurately provide relative importance measures of each factor. Vegetation growth is influenced by multiple factors, among which precipitation, temperature, and radiation are crucial. To further clarify the dominant influencing factors in different regions of the YRB, this study used the IAV time series of annual mean temperature, precipitation, and radiation as independent variables, and the IAV time series of annual mean NPP as the dependent variable. The LMG method was used to quantitatively evaluate the relative contributions of each factor to the NPP IAV. The relative contributions of the three factors were normalized, and then the RGB combination was used to obtain the spatial distribution pattern of the relative contributions of the three factors to visually depict the spatial patterns of the dominant factors in different regions. R software (V4.3.1) and Python software (V3.9.13) were used in this study.

3. Results

3.1. Model Performance

Our simulation showed a significant positive correlation with the MODIS NPP, with an R² of 0.76 (p < 0.001) and RMSE (root mean square error) of only 65.57 g C m⁻². The PB (percentage bias) between the NPP simulated using the CASA model and the MODIS NPP was only 8.31% (Figure 2a), indicating good consistency between our simulation results and mainstream satellite retrievals in terms of the spatial pattern.

We further compared our results with 220 in situ annual NPP measurements over the YRB published in an open access database and found that our results were in good consistence with the measurements (R² = 0.75, p < 0.001) (Figure 2b). The RMSE and PB between the estimated and observed NPP were 72.77 g C m⁻² and 15.72%, respectively. Although the MODIS NPP also kept consistent with the observations (R² = 0.51, p < 0.001), our estimates showed a better performance than the MODIS NPP when compared with the in situ observations in terms of the PB and RMSE (Figure 2b).

We further compared our results with the simulations of 14 TBMs (Figure 3). The comparison revealed significant differences in magnitudes among different TBM simulations. Between 1999 and 2010, the annual mean NPP of the YRB simulated by these TBMs ranged from approximately 40 gC m⁻² (VEGAS2.1) to 450 gC m⁻² (GTEC), with the maximum value being about 11 times the minimum simulation. During this period, our estimate fell within the middle of the distribution of these model simulations. The mean NPP we simulated (233.63 gC m⁻² yr⁻¹) closely aligned with the ensemble mean of these models (222.26 gC m⁻² yr⁻¹), with an RMSE of only 18.82 gC m⁻² yr⁻¹ and a PB of only 5.09%. Additionally, our estimate shows high consistence in the temporal trend (r = 0.61, p < 0.05) with the ensemble mean of these TBMs. Furthermore, when compared to the simulation of CLM5.0 over a longer time span (1999—2018), it was found that the CLM5.0 estimates were lower than those of most (11 out of 14) of the TBMs, as well as lower than our results. However, we observed a good consistency in the temporal trends between our estimates and CLM5.0 simulations (r = 0.69, p < 0.01). Overall, the above results indicated that our results based on the CASA can reasonably reflect the spatiotemporal patterns of the NPP in the YRB.

3.2. Spatial Pattern of NPP

The spatial distribution of the average annual NPP in the YRB from 1999 to 2018 is shown in Figure 4a. The high NPP values are concentrated in the southern part of the basin, the Qinling Mountains, and the northeastern part of the Tibetan Plateau, while the low values are distributed in the upstream Tibetan Plateau and the northern part of the Loess Plateau in the middle reaches. Cell grids with an annual NPP between 200 and 400 g C m⁻² yr⁻¹ account for the largest proportion of the vegetation area in the basin (47.43%), followed by areas with an annual NPP below 200 g C m⁻² yr⁻¹ (38.62%) and areas with an NPP between 400 and 600 g C m⁻² yr⁻¹ (12.51%). The multi-year average NPP in the YRB was about 257.02 g C m⁻² yr⁻¹, with an annual total primary production of 189.81 Tg C.

According to the distribution of the NPP at different altitudes (Figure 4b), areas with altitudes below 1 km and between 2 and 4 km have a higher annual mean NPP (>300 g C m⁻² yr⁻¹) among all altitude gradients. Areas with an altitude above 4 km show the lowest annual mean NPP (132.55 g C m⁻² yr⁻¹). In terms of the annual total NPP, regions between 1 and 1.5 km (the major implementation areas of the GFGP) account for the largest proportion (33.8%) among all altitude gradients, followed by the regions between 1.5 and 2 km (14.52%).

From the perspective of climate zones (Figure 4c), the semi-humid region accounts for the largest area proportion (45.07%) of the YRB. The annual mean NPP in the humid region is the largest (482.83 g C m⁻² yr⁻¹). However, this region only covers a very small area (0.23%) of the YRB, and its contribution to the total NPP of the basin is quite small (0.82 Tg C yr⁻¹, or 0.43%). Although the annual mean NPP of the semi-humid region (325.17 g C m⁻² yr⁻¹) is only 67.4% of that of the humid region, its contribution to the total NPP of the basin is the largest (57.03%), and the annual mean NPP in the arid region is the lowest (133.62 g C m⁻² yr⁻¹). In terms of the total NPP, the Tibet Plateau and the semi-arid region also show non-negligible contributions to the YRB with rates of 24.30% and 14.80%, respectively. The above results highlight the important role of the semi-humid region in controlling the vegetation carbon sink capacity of the YRB.

In terms of vegetation types (Figure 4d), woodland shows the largest annual mean NPP (348.57 g C m⁻² yr⁻¹), followed by cropland (292.40 g C m⁻² yr⁻¹) and grassland (224.95 g C m⁻² yr⁻¹). Grassland accounts for the largest contribution (44.23%) to the total NPP in the YRB, followed by cropland (31.12%) and woodland (19.31%).

3.3. Spatiotemporal Dynamics of NPP

Figure 5a illustrates the spatial pattern of annual NPP change rates in the YRB. The NPP change rates in the study area range from −24 to 35 g C m⁻² yr⁻¹ and the whole study area shows an increasing trend in NPP with a rate of 6.76 g C m⁻² yr⁻¹. Regions with the largest increasing rates (>12 g C m⁻² yr⁻¹) are concentrated around the Qinling Mountains and the central Loess Plateau. Most parts (79.63%) of the YRB shows significantly (p < 0.05) increasing trends in NPP (Figure 5b), especially the Qilian Mountains and the key implementation areas of the GFGP. Areas where the NPP significantly decrease are mostly urban areas, including the Guanzhong Plain urban agglomeration and cities (e.g., Luoyang) downstream of the YRB.

Figure 5c shows the temporal variations of the NPP in the whole YRB and different segmentations. The NPP in the upstream, midstream, and downstream basins of the YRB all show fluctuating upward trends. The largest increasing rate (8.4 g C m⁻² yr⁻¹, p < 0.01) occurs in the midstream basin, followed by the downstream basin (3.81 g C m⁻² yr⁻¹, p < 0.01) and upstream basin (2.34 g C m⁻² yr⁻¹, p < 0.05).

Figure 5d reveals the NPP IAV in the whole YRB and different segmentations. The NPP IAV in the downstream basin is relatively stable compared with other segmentations. It is clear that the NPP IAV in the midstream basin showed the largest contribution (62.36%) to that of the whole YRB, followed by the upstream basin (36.25%). It is worth noting that in 2011, the whole YRB, midstream basin, and upstream basin all showed the largest negative IAV during the study period with magnitudes of −25.34 Tg C, −14.75 Tg C, and −9.94 Tg C, respectively. This large negative variation may be attributed to the extreme drought event during this year [75].

3.4. Regional Contributions to the Holistic NPP IAV

Figure 6a shows the spatial pattern of the cell-scale contributions of the NPP IAV to the holistic NPP IAV of the YRB. Most regions of the basin (96%) showed a positive contribution to the basin holistic NPP IAV. Only a few regions, for instance, the urban agglomeration in the south and southeast of the midstream basin, showed negative impacts on the holistic NPP IAV of the YRB. As for different climate regions (Figure 6b), the contribution rate of the semi-humid region was the highest (62.66%), followed by the semi-arid region (18.19%) and the Tibetan Plateau (15.74%). We further analyzed the contribution rate per unit area and found that the contributions of the semi-humid (1.31 × 10⁻⁴% km⁻²) and humid regions (1.30 × 10⁻⁴% km⁻²) were very close, both being substantially higher than those of other climate regions. Overall, the semi-humid area dominated the NPP IAV of the YRB.

In terms of the contributions of different altitudes (Figure 6c), areas with altitudes between 1 and 1.5 km contributed the largest (43.07%) to the basin holistic NPP IAV, followed by areas with altitudes between 1.5 and 2 km (17.77%) and 0.5 and 1 km (12.24%). In terms of the contribution rate per unit area, the areas between 0.5 and 3 km shared similar contributions rates (around 1.1 × 10⁻⁴% km⁻²).

We further evaluated the contribution rates of different vegetation types to the NPP IAV in the YRB (Figure 6d). The results showed that the contribution rate of grassland was the largest, reaching 45.69%, followed by cropland and woodland, which were 31.21% and 45.69%, respectively. However, woodland showed the largest contribution rate per unit area (1.26 × 10⁻⁴% km⁻²), followed by cropland (1.08 × 10⁻⁴% km⁻²) and grassland (0.85 × 10⁻⁴% km⁻²). Figure 6e shows the contribution rates of the area under the GFGP and the other area. The GFGP area contributes to 19.57% of the holistic NPP IAV in the YRB and its contribution rate per unit area (1.25 × 10⁻⁴% km⁻²) is larger than that (0.89 × 10⁻⁴% km⁻²) of the other area.

3.5. Attribution of the NPP IAV

We quantified the relative contributions of the key drivers (radiation, precipitation, and temperature) of the NPP IAV based on the LMG method (Figure 7a). It can be found that the NPP IAV in most parts of the YRB was influenced by multiple factors. Specifically, the southern upstream and southeastern midstream regions were mainly affected by both temperature and radiation, while the central and northern upstream and northeastern midstream regions were mainly affected by precipitation and radiation. The central upstream region was mainly affected by precipitation and temperature. The areas affected only by precipitation were mainly concentrated in the midstream region, while areas dominated by temperature were concentrated in the upstream region, and the southern midstream region was primarily affected by radiation.

By extracting the dominant factors of each grid cell (Figure 7b), we found that water availability was the most important factor for NPP variability in the entire basin, with a relative contribution rate of 51.31%, followed by temperature (27.73%) and radiation (20.96%) (Figure 7d). Precipitation dominated the NPP IAV in 62.90% of the YRB, followed by temperature (23.07%) and radiation (14.03%) (Figure 7b). As for the upstream, midstream, and downstream basins, precipitation was the dominant factor for the NPP IAV in these regions, showing relative contributions that ranged from 43.31% to 56.69% (Figure 7d) and dominant areas of 66.31%, 59.81%, and 47.22%, respectively. Temperature was the second important factor in the upstream and downstream region, showing dominant areas of 28.96% and 38.55%, respectively. Radiation was the second important factor affecting the NPP IAV in the midstream region, dominating an area of 25.15%, especially in the Guanzhong Plain (Figure 7b).

As for different climate zones, precipitation is the dominant factor in the semi-humid, semi-arid, and arid zones, showing relative contributions of 51%, 80%, and 75%, respectively. In other words, water availability controls the NPP IAV in water-limited regions. In the humid zone, radiation shows the largest relative contribution (73%), which means energy controls the NPP IAV in regions without water and heat limitation. In the Tibetan Plateau region of the YRB, temperature and precipitation are both important drivers of the NPP IAV, sharing a similar relative contribution (about 40%).

4. Discussion

4.1. Spatiotemporal Changes of NPP in the YRB

This study simulated the spatiotemporal patterns of the NPP in the YRB using the CASA model. We found that the spatial distribution of the NPP in the basin exhibited a gradually decreasing pattern from south to north, which was consistent with the spatial pattern of rainfall. These results confirm the previous findings on the spatial distribution of vegetation growth and carbon sinks in the YRB [21,59,60]. Comparison of our simulation results with the MODIS NPP and in situ measurements at different spatiotemporal scales also confirmed the spatial distribution pattern of the NPP in the basin (Figure 2).

We further analyzed the temporal evolution of the NPP and found a significant increasing trend in NPP for more than 98% of the basin area after 1999. The vegetation in the middle stream of the YRB showed the fast recovery rate, with an NPP growth rate of 8.4 g C m⁻² yr⁻¹ (p < 0.001). Our results confirmed the vegetation carbon sink effect of the ecological restoration project in the YRB.

In addition, we systematically analyzed the IAV of the NPP and found that the largest positive IAV occurred in 2004 for the overall basin and the middle stream. The extreme precipitation in 2003, the largest since 1965 [76], may have contributed to the high positive IAV in 2004 via the lagged effects of vegetation response to water availability [77]. The abundant precipitation in 2003 may have stimulated the growth and productivity of vegetation in the following year. Meanwhile, the overall basin, upstream, and middle stream all exhibited the largest negative IAV in 2011, possibly due to the extreme drought that occurred in this year [75].

Furthermore, through analysis of the regional contributions of different vegetation types in the YRB, we found that grasslands contribute the most to the NPP IAV due to their extensive coverage, while forests exhibit the highest per-unit-area contribution rate. Figure 6e also demonstrates the significant impact of the GFGP on NPP IAV. It is worth noting that previous studies indicated that the forest net carbon gain (NPP) increases rapidly in youth, peaks in middle age, and then decreases in old age [78,79]. Forest aging may lead to a future decline in the NPP [80]. Therefore, selecting the right forest species for maximal carbon stock and progressive forestation may be an appropriate way to reduce CO₂ and compensate for the carbon sink decline in existing forests [81].

4.2. Regional Contribution of NPP IAV

Our results showed that most of the areas in the YRB positively contributed to the basin’s holistic NPP IAV. The semi-humid region showed the largest contribution (62.66%). Our results differ from previous global-scale studies [33,37,38], which indicated that semi-arid areas or tropical regions controlled the vegetation productivity worldwide. However, our findings are consistent with a previous study, which also identified the semi-humid area as a major contributor to the NPP IAV in China [21]. For the YRB, the GFGP was mainly implemented in the semi-humid region, which is a key reason for the rapid recovery of vegetation in this region over the past two decades [6,49,53]. Our results further highlight the important role of this region on the NPP IAV in the YRB.

4.3. Response of NPP IAV to the Key Drivers

Our study found that in most areas of the YRB (including the western upstream and eastern midstream regions), especially in the Loess Plateau region, the NPP IAV was mainly affected by precipitation. This is because that these regions are located in the arid, semi-arid, and semi-humid regions, which are characterized by an uneven spatiotemporal distribution of precipitation and high evapotranspiration, leading to increased sensitivity of the NPP IAV to water availability [44,47]. The southern upstream region of the basin is located in the Tibetan Plateau, where radiation is strong and there is a clear division between the wet and dry seasons. However, the temperature in this region is low throughout the year, and water is limited; hence, the NPP IAV is mainly influenced by temperature and precipitation [43,82]. The southern midstream region of the basin has sufficient water and heat conditions; thus, it is mainly controlled by energy [21]. Other influencing factors include the soil and aboveground biomass. Soil provides essential nutrients and moisture for plant growth [83]. Soil properties such as texture, pH, and organic matter content directly affect plant growth and root expansion [84]. Additionally, aboveground biomass, including trees, herbs, and other vegetation, influences vegetation competitiveness, growth, and diversity by intercepting sunlight, reducing water evaporation, and providing a habitat [85]. In summary, changes in soil quality and aboveground biomass can also significantly impact vegetation productivity, thereby exerting a widespread influence on ecosystems and agricultural yields.

4.4. Uncertainties and Limitations

Despite the effective evaluation of the simulation results, there are still some uncertainties in this study. First, uncertainties may come from the satellite observations and model estimations of the NPP IAV, which are crucial when trying to accurately estimate changes in vegetation production and the impacts of different factors. Reducing the uncertainties of satellite observations and model estimations should be the focus of further research. Second, our study did not fully consider the impact of extreme weather events on the NPP IAV. For example, heavy precipitation, high temperature, non-extreme precipitation, and low-temperature events may show considerable impacts on the IAV of vegetation production [21,86,87,88].

5. Conclusions

This study, employing the CASA model, conducted a comprehensive assessment of the NPP in the YRB for the period spanning from 1999 to 2018. Our findings demonstrated a remarkable consistency with the spatial patterns observed in the MODIS NPP, and notably, our estimates exhibited a better performance when compared with in situ observations. The average annual total NPP in the YRB was 189.81 Tg C yr⁻¹ during the study period. The NPP in the YRB exhibited distinct spatial patterns with elevated NPP concentrated in the southern region of the basin, the Qinling Mountains, and the northeastern part of the Tibetan Plateau. In contrast, lower values were predominantly situated in the northern part of the Loess Plateau. Key contributors to the basin’s overall carbon sink were areas with altitudes ranging from 1 to 1.5 km, the semi-humid region, and grasslands. Remarkably, approximately 79.63% of the YRB displayed statistically significant (p < 0.05) increasing trends in NPP, with particular prominence observed in the Qilian Mountains and the primary implementation areas of the GFGP. The NPP exhibited fluctuating upward trends in the upstream, midstream, and downstream basins, with the midstream basin registering the most substantial increasing rate of 8.4 g C m⁻² yr⁻¹ (p < 0.01). This region also made the most significant contribution (62.36%) to the NPP IAV of the entire YRB. Spatially, 96% of the YRB exhibited a positive contribution to the holistic NPP IAV, with the semi-humid region exhibiting the largest contribution (62.66%). These results underscore the crucial role of the semi-humid region in governing both the magnitude and IAV of NPP in the YRB. In addition, we identified varying drivers influencing the NPP IAV across the YRB. Precipitation emerged as the dominant factor, exerting a positive influence on 62.90% of the entire basin, especially in the water-limited regions. In summary, this study provides a framework for detecting interannual variation in NPP and understanding the multifaceted impacts of future climate fluctuations on the persistence of vegetation production in the YRB.