## Materials and Methods
### Ecosystem Quality Estimation  
The FVC was estimated using the pixel dichotomy model [31], which assumes that image pixels are composed of two parts: vegetation-covered and non-vegetation-covered areas [46], and thus the FVC can be defined as follows:  
\[\mathrm{FVC}=\frac{\mathrm{NDVI}-\mathrm{NDVI}_{\mathrm{soil}}}{\mathrm{NDVI} _{\mathrm{veg}}-\mathrm{NDVI}_{\mathrm{soil}}} \tag{2}\]  
where \(\mathrm{NDVI}_{\mathrm{veg}}\) and \(\mathrm{NDVI}_{\mathrm{soil}}\) are the Normalized Difference Vegetation Index (NDVI) value for pure vegetation pixels and pure bare soil pixels, respectively. We used MOD13A3 data, monthly MODIS NDVI at a 1 km spatial resolution as a gridded level-3 product, as the calculation input after aggregation through the Maximum Value Composite (MVP) approach [47] to obtain the annual maximum value and after time domain processing using a Savitzky-Golay filter [48] to minimize the influence of existing noise.  
The LAI was derived from the MODIS MOD15A2 product, which is composited every eight days at a 1 km resolution. We improved our LAI data quality using the temporal spatial filter (TSF) [49] to fill in the spatial and temporal discontinuity of original datasets. We then used MVP to obtain the annual maximum value of LAI.  
NPP was evaluated by applying the GLOPEM-CEVSA model [50], which couples a satellite-based model and a process-based model. The model has been validated by carbon flux observations within forest, grassland, and cropland. In the model, NPP is the difference of gross primary productivity (GPP) and autotrophic respiration (R\({}_{\mathrm{a}}\)), where GPP is estimated by the GLoPEM model, and R\({}_{\mathrm{a}}\) is divided into maintenance respiration (R\({}_{\mathrm{m}}\)) and growth respiration (R\({}_{\mathrm{g}}\)) from four plant parts--thick root, thin root, stem, and leaf--with reference to the process-based BEPS model and Biome-BGC model:  
\[\mathrm{GPP}=\mathrm{PAR}\times\mathrm{FPAR}\times\epsilon_{\mathrm{g}} \tag{3}\]  
\[\mathrm{R_{a}}=\sum_{\mathrm{i=1}}^{4}\mathrm{R_{m,i}}+\mathrm{R_{g}} \tag{4}\]  
\[\mathrm{R_{m,i}}=\mathrm{M_{i}r_{m,i}Q_{10}}^{(\mathrm{T-Tb})/10} \tag{5}\]  
\[\mathrm{R_{g}}=\mathrm{r_{g}GPP} \tag{6}\]where PAR (photosynthetically active radiation) was estimated using climatology radiation calculation methods [51] and FPAR (fraction of photosynthetically active radiation) was adapted from satellite remote sensing inversion data based on the canopy radiative transfer model [52]. \(\epsilon_{\mathrm{g}}\), the vegetation light utilization efficiency in the GPP concept, was the product of the maximum possible light utilization efficiency and environment stresses including vapor pressure deficit (VPD), PAR, air temperature, atmospheric CO\({}_{2}\), and soil moisture [53]. M\({}_{\mathrm{i}}\) and r\({}_{\mathrm{m,i}}\) are the biomass and maintenance respiration coefficients of the four plant organs, respectively. Q\({}_{10}\) represents the change in the maintenance respiration rate when the temperature changes by 10 degrees; T and T\({}_{\mathrm{b}}\) are the actual temperature and base temperature, respectively; and r\({}_{\mathrm{g}}\) is the proportional coefficient of total growth respiration to total plant growth. The input meteorological data, including temperature, precipitation, evaporation, and total radiation, were provided by China's State Meteorological Administration. They were interpolated into grids using thin plate smoothing splines (ANUSPLIN) [54; 55].