*4.1. Credibility of HANPP*

In existing studies, HANPP (%NPPpot) has attracted more attention than HANPP (gC/m<sup>2</sup> ) or HANPP (gC) due to its ability to represent regional ecological sustainability and the ecological pressure caused by human occupation. However, the lack of a uniform standard in the selection of estimation models of NPPpot and NPPact hinders the comparison of HANPP results among different regions. For instance, in related studies, a variety of models have been used to estimate NPPpot, e.g., LPJ-DGVM [16,41], the Guangsheng Zhou model [5], the new calibrated Miami model [45], and multiple-iterations modification based on the disturbed pixels in NPPact [2]. Rather than comparing the similarities of HANPP (%NPPpot) between regions, it is better to focus on whether its spatial distribution and long-term variation characteristics are similar.

In China, HANPP (%NPPpot) was found to be 60.33% in 2015, which is higher than it was in 2010 (57.80%), estimated by Chen et al. (2015) [5], as well as the global HANPP (%NPPpot) values in 2000 (23.80%), 2005 (25.00%), and 2010 (20.70%) [16–18]. Huang et al. (2020) found that HANPP (%NPPpot) was 72.4% in the Yangtze River Delta in 2015 [51], and we found a slightly lower value (68.92%) in the same region in this study. From the perspective of the spatial distribution characteristics of HANPP, the grid-scaled HANPP (%NPPpot) in this study is highly similar to the results of Haberl et al. (2007) and Kastner et al. (2021) [16,17]. In general, the simulation of HANPP among provinces and cities in China is still in its infancy and the deviations are mainly due to different methods and datasets used during NPPpot and HANPPharv estimation within limited studies.

### *4.2. Relationships between HANPP and China's Urban and Rural Development*

Significant spatial differentiation of HANPP was found under the impact of unbalanced human activities and the ecological environment in China, especially the duality of rural–urban development. In this study, 20 indictors reflecting regional economy, population, urban expansion, and rural agricultural technology were chosen, as shown in Table 1, and the correlation coefficients and nine HANPP indexes (Table 3) of 31 provinces were summarized (Figure 4). From the perspective of urban economy, significant positive correlations were identified between the added value of the secondary and tertiary industries (UE1) and HANPP (gC/m<sup>2</sup> ), HANPPharv (gC/m<sup>2</sup> ), and HANPPluc (gC/m<sup>2</sup> ). The added value of transportation, warehousing, and postal services (UE3) had a significant influence on most HANPP indexes (Figure 4a), which guarantees the circulation of regional ecological resources. However, UE2 and UE4 showed a negative correlation with the total amount of HANPP, HANPPharv, and HANPPluc, which is because a larger proportion of UE1 and a higher consumption expenditure of urban residents always appears in highly urbanized areas, which usually have a relatively weak biomass production (low HANPPharv); the HANPPluc would also not increase significantly when urban construction is saturated, so the HANPP indexes would be suppressed instead. High urban population (UP1) resulted in a great amount of biomass occupation per unit area, whereas the negative correlations between urban population proportion (UP2) and HANPP (gC), HANPPharv (gC), and HANPPluc (gC) reflected the polarization characteristics of the urban population agglomeration and ecological resource production. Among the urban spatial expansion indicators, the higher regional average nighttime light intensity (US1), the proportion of the built-up area (US2), and the proportion of the urban road area (US4) always corresponded with large HANPPluc (% NPPpot) and HANPPluc (gC/m<sup>2</sup> ) values, which also directly led to the significant positive correlations (*p* < 0.01) and relatively high correlation coefficients with HANPP (% NPPpot) and HANPP (gC/m<sup>2</sup> ). Conversely, urban ecological construction gained valuable space to counter the degradation of natural vegetation caused by urban expansion, and most of the HANPP indexes, especially the HANPPluc indexes, decreased as the proportion of urban green spaces increased (US3).


**Figure 4.** Relationship between HANPP indexes and urban (**a**) and rural (**b**) development indicators. (UE1–4, UP1–2, and US1–4 and RE1–4, RP1–2, and RA1–4; refer to Table 1; \*: *p* < 0.1, \*\*: *p* < 0.05, and \*\*\*: *p* < 0.01). **Figure 4.** Relationship between HANPP indexes and urban (**a**) and rural (**b**) development indicators. (UE1–4, UP1–2, and US1–4 and RE1–4, RP1–2, and RA1–4; refer to Table 1; \*: *p* < 0.1, \*\*: *p* < 0.05, and \*\*\*: *p* < 0.01).

HANPP indexes responded to rural development in different ways (Figure 4b). Stronger positive correlations existed between the added value of the primary industry (RE1) and all HANPPharv indexes than UE1, and the total amount of HANPP, HANPPharv, and HANPPluc in each province was more sensitive to the proportion of the added value of the primary industry (RE2) than other HANPP indexes. Higher consumption expenditure of rural residents (RE3) appeared not only in highly urbanized areas but also in the agriculture, forestry, and animal husbandry developed areas, which showed relatively lower correlations with the total amount of HANPP, HANPPharv, and HANPPluc than UE4. The higher the Engel coefficient of rural residents (RE4), the lower the HANPP (% NPPpot) and all HANPPharv indexes. This shows that the main article consumed by residents was food, the regional economy was relatively lagging in terms of development, and the influence of human activities on ecology was relatively weak. Rural population (RP1) had significant and higher positive correlations with all the HANPPharv indexes than UP1 but no statistical correlations with HANPPluc indexes. Moreover, the biomass occupied through land use per unit area and its percentage in NPPpot would be reduced if the rural population was the main component (RP2) in a region. All the rural agricultural technology indicators showed similar impacts on the nine HANPP indexes, especially in terms of the significant positive effects (*p* < 0.01) on and high correlation coefficients among HANPPharv (gC), HANPPharv (% NPPpot), and HANPPharv (gC/m<sup>2</sup> ). Adequate effective irrigated area (RA3) and total sown area of crops (RA4) are necessary for adequate food supply, the total power of the agricultural machinery (RA1) reflects the regional agricultural machinery ownership, and the harvest biomass of grain and forestry would benefit from these elements as well as the amount of agricultural fertilizer (RA2) applied.

An interaction effect always exists between regional economy improvement and biomass occupation, and industrial structure would alter the HANPP proration. Overdevelopment of the secondary and tertiary industries would restrain harvest biomass, but adequate ecological supply depends on the primary industry, so balanced industrial development should be the focus among regions. Harvested biomass (HANPPharv) is mostly determined by the rural population rather than the urban population, while the urban or rural population proportion has hardly any influence on it, indicating that urban people mainly play the role of the consumer of ecological resources. Urban sprawl, a result of land surface modification by humans, usually leads to the irrecoverable loss of NPP. Cultivated land offers additional products along with the improvement of agricultural technology, although it replaces part of the natural vegetation. Moreover, worldwide, the increasing populations and economies are inevitably accompanied by rapid urbanization, and the corresponding requirements of ecological resources must be sustainable. Thus, through policies and agriculture techniques, it is necessary to improve agriculture production and the utilization efficiency of NPP per area, as well as properly control urban expansion.

#### *4.3. HANPP Responses to Coordinated Regional Urban–Rural Development*

As two major systems indispensable to regional development, urban areas and rural areas play an important role in promoting sustainable use of ecological resources. Therefore, coordinated urban–rural development and the impact on the HANPP in 31 provinces in China was quantified. The quantification of the urban development index (*f*(*u*)) indicated that highly urbanized areas (Figure 5a) were mostly concentrated in the coastal provinces but high values of the rural development index (*f*(*r*)) mostly appeared in southern China (Figure 5b). Specifically, Shanghai (*f*(*u*) = 0.90), Beijing (*f*(*u*) = 0.43), Henan (*f*(*r*) = 0.73), and Shandong (*f*(*r*) = 0.68) showed the peak values over the 95th percentiles of *f*(*u*) and *f*(*r*). In contrast, the valley values under the fifth percentiles were in Qinghai (*f*(*u*) = 0.04), Tibet (*f*(*u*) = 0.02), Tianjin (*f*(*r*) = 0.09), and Beijing (*f*(*r*) = 0.08). Furthermore, the degree of coordinated urban–rural development (*D*) among provinces was evaluated by the comprehensive urban– rural development index (T) and the coupling degree (C) (Table S4) and showed an obvious spatial differentiation between coastal and inland areas (Figure 5c). In this study, the *D* value was divided into nine levels (Table 2) to map coordinated urban–rural development

types in China (Figure 5d). Only a few coastal provinces had achieved primary coordination (Jiangsu and Guangdong) and barely coordination (Shandong and Zhejiang). Inversely, most regions showed a universal imbalance and urban–rural development gap. Notably, there was severe imbalance in Heilongjiang Province, Tibet, Gansu Province, and Xinjiang Province. In short, the ideal urban–rural coordination did not exist in China in 2015. achieved primary coordination (Jiangsu and Guangdong) and barely coordination (Shandong and Zhejiang). Inversely, most regions showed a universal imbalance and urban– rural development gap. Notably, there was severe imbalance in Heilongjiang Province, Tibet, Gansu Province, and Xinjiang Province. In short, the ideal urban–rural coordination did not exist in China in 2015.

5c). In this study, the *D* value was divided into nine levels (Table 2) to map coordinated urban–rural development types in China (Figure 5d). Only a few coastal provinces had

*Land* **2023**, *12*, 1062 16 of 22

**Figure 5.** Urban and rural development index (**a**,**b**) and coordinated urban–rural development degree (**c**) and types (**d**) in China. HANPP responses to coordinated regional urban–rural development (**e**). (PC, primary coordination; BC, barely coordination; NI, near imbalance; II, inchoate imbalance; MI, moderate imbalance; SI, severe imbalance; \*\*: *p* < 0.05, and \*\*\*: *p* < 0.01). **Figure 5.** Urban and rural development index (**a**,**b**) and coordinated urban–rural development degree (**c**) and types (**d**) in China. HANPP responses to coordinated regional urban–rural development (**e**). (PC, primary coordination; BC, barely coordination; NI, near imbalance; II, inchoate imbalance; MI, moderate imbalance; SI, severe imbalance; \*\*: *p* < 0.05, and \*\*\*: *p* < 0.01).

The relationship between HANPP indexes and coordinated urban–rural development (Figure 5e) helped create verifications and supplements for Figure 4. It is worth noting that *f(u)* and *f(r)* showed opposite driving effects on HANPP (gC), HANPPharv (gC), and HANPPluc (gC) and the regional ecological pressure characterized by HANPP (%NPPpot) and HANPP (gC/m<sup>2</sup> ) mostly depended on *f(r)*. The highly urbanized regions with high *f(u)* usually appeared in provinces with limited territorial areas (i.e., Shanghai, Beijing, and Tianjin), in which the total amount of occupied biomass was at a relatively low level. On the contrary, provinces with high levels of *f(r)* had relatively large territorial areas in China and developed production capacity. As a component of HANPP, HANPPluc was mainly controlled by the urban development index (*f(u)*). Conversely, *f(r)* made significant positive impacts on HANPPharv. All the occupied biomass indexes per unit area showed the most sensitive response to coordinated urban–rural development degree (*D*) and good urban–rural coordination obviously improved HANPP (gC/m<sup>2</sup> ) and HANPPharv (gC/m<sup>2</sup> ) but relatively weaker urban–rural coordination influenced HANPPluc (gC/m<sup>2</sup> variation. From the perspective of six coordination types among provinces in China, the sample points of HANPPharv (gC/m<sup>2</sup> ) tended to be discrete, along with a higher degree of coordination (Figure 6), and the samples with moderate imbalance (MI) urban–rural coordination showed the highest degree of dispersion of HANPP (gC/m<sup>2</sup> ) and HANPPluc (gC/m<sup>2</sup> ). The relationship between HANPP indexes and coordinated urban–rural development (Figure 5e) helped create verifications and supplements for Figure 4. It is worth noting that *f*(*u*) and *f*(*r*) showed opposite driving effects on HANPP (gC), HANPPharv (gC), and HANPPluc (gC) and the regional ecological pressure characterized by HANPP (%NPPpot) and HANPP (gC/m<sup>2</sup> ) mostly depended on *f*(*r*). The highly urbanized regions with high *f*(*u*) usually appeared in provinces with limited territorial areas (i.e., Shanghai, Beijing, and Tianjin), in which the total amount of occupied biomass was at a relatively low level. On the contrary, provinces with high levels of *f*(*r*) had relatively large territorial areas in China and developed production capacity. As a component of HANPP, HANPPluc was mainly controlled by the urban development index (*f*(*u*)). Conversely, *f*(*r*) made significant positive impacts on HANPPharv. All the occupied biomass indexes per unit area showed the most sensitive response to coordinated urban–rural development degree (*D*) and good urban–rural coordination obviously improved HANPP (gC/m<sup>2</sup> ) and HANPPharv (gC/m<sup>2</sup> ) but relatively weaker urban–rural coordination influenced HANPPluc (gC/m<sup>2</sup> ) variation. From the perspective of six coordination types among provinces in China, the sample points of HANPPharv (gC/m<sup>2</sup> ) tended to be discrete, along with a higher degree of coordination (Figure 6), and the samples with moderate imbalance (MI) urban–rural coordination showed the highest degree of dispersion of HANPP (gC/m<sup>2</sup> ) and HANPPluc (gC/m<sup>2</sup> ).

)

**Figure 6.** Correlations between HANPP indexes and the degree of coordinated urban–rural development. **Figure 6.** Correlations between HANPP indexes and the degree of coordinated urban–rural development.

#### *4.4. Dominant Driving Factors of Each HANPP Index 4.4. Dominant Driving Factors of Each HANPP Index*

Bidirectional effects were found among urban–rural development indictors of biomass appropriation, and the dominant driving factors of each HANPP index were detected by stepwise regression to eliminate the indirect correlation. Table 4 shows the fitting models and correlation coefficients from three perspectives (urban, rural, and urban–rural development), in which stronger decisive effects were displayed by the driving factor with larger standardized coefficients, and a higher R referred to a better fitting model. Among all the driving factors in the 27 fitting equations in Table 4, the rural development indicators (34 times) appeared more frequently than the urban development indicators (21 times). At the same time, in different dimensions of urban–rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE), more likely to be affected by rural population (RP) than urban population (UP), and more likely to be affected by rural agricultural technology (RA) than urban spatial expansion (US). From the analog effect of each equation, driven by 20 urban and rural development indicators, the analog effect of HANPPharv indexes was the best, followed by that of HANPP indexes, and HANPPluc indexes were the weakest; the selected urban and rural development indicators could explain their variation characteristics only to a certain extent. From the perspectives of urban, rural, and comprehensive urban–rural development, Bidirectional effects were found among urban–rural development indictors of biomass appropriation, and the dominant driving factors of each HANPP index were detected by stepwise regression to eliminate the indirect correlation. Table 4 shows the fitting models and correlation coefficients from three perspectives (urban, rural, and urban–rural development), in which stronger decisive effects were displayed by the driving factor with larger standardized coefficients, and a higher R referred to a better fitting model. Among all the driving factors in the 27 fitting equations in Table 4, the rural development indicators (34 times) appeared more frequently than the urban development indicators (21 times). At the same time, in different dimensions of urban–rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE), more likely to be affected by rural population (RP) than urban population (UP), and more likely to be affected by rural agricultural technology (RA) than urban spatial expansion (US). From the analog effect of each equation, driven by 20 urban and rural development indicators, the analog effect of HANPPharv indexes was the best, followed by that of HANPP indexes, and HANPPluc indexes were the weakest; the selected urban and rural development indicators could explain their variation characteristics only to a certain extent. From the perspectives of urban, rural, and comprehensive urban–rural development,

the analog effects of all the nine HANPP indexes were the best when only considering rural development indicators, followed by the regression models of comprehensive urban– rural development indicators, and the ability to explain HANPP indexes was relatively limited when only considering urban development indicators. From the perspective of the explanatory ability of each independent variable, the impact of urban economy (UE) on occupied biomass was much stronger than that of urban population (UP) and urban spatial expansion (US). On the contrary, rural population (RP) and rural agricultural technology (RA) had more significant explanatory effects on each HANPP index than rural economy (RE). Among all the independent variables, UP2, US1, US2, US3, and RA3 did not appear in any of the fitting models.


**Table 4.** Dominant driving factors of each HANPP index.

Comprehensively considering the dominant driving factors and standardized coefficients of each HANPP index, it was found that UE3 had the strongest explanation for the differentiation characteristics of HANPP (gC/m<sup>2</sup> ) and that the restriction of RP and RA indicators on HANPP (%NPPpot) was much stronger than that of urban development indicators. The total amount of HANPP (gC) was affected by the opposite effects of RA4 and US4, and RA4 showed a stronger control of it. In addition, RP and RA indicators showed absolute driving advantages in HANPPharv (%NPPpot) and HANPPharv (gC) and RA indicators had a stronger positive driving effect on them. From the perspective of comprehensive urban–rural development, urban development indicators had no significant driving effect on the above two indexes. The influence of each independent variable on HANPPharv (gC/m<sup>2</sup> ) was relatively complicated, in consideration of RA, UE, and US indicators, and a better-fitting model could obtain HANPPharv (gC/m<sup>2</sup> ); the standardized coefficient showed that the urban development indicators played a more important role. Table 4 shows that the fitting equations between the independent variables and HANPPluc were relatively simple. US4, RE2, and UE1 were the dominant driving factors for HANPPluc

(gC), HANPPluc (%NPPpot), and HANPPluc (gC/m<sup>2</sup> ), respectively, and none of the population factors were correlated. It should be noted that the ability of selected independent variables to interpret all HANPPluc indexes was significantly weaker than their ability to interpret HANPP and HANPPharv, indicating that there are driving factors besides those selected in this study. Further analysis of the positive or negative driving differences of each independent variable showed that UE3, RA1, RA2, and RA4 had a significant positive influence in each fitting model, which could be verified and compared with the results in Figure 4. In summary, HANPP (%NPPpot), HANPPharv (%NPPpot), and HANPPluc (%NPPpot) were more dependent on rural development, while urban development was the decisive factor for all HANPP indexes per unit area, and the total amount of occupied biomass of each province was more affected by rural development than urban.

### **5. Conclusions**

In this study, we estimated the HANPP in China in 2015 based on multi-source datasets and quantified the contribution characteristics of HANPP components as well as the biomass production–consumption types in 31 provinces. The regional differentiations of HANPP indexes among typical regionalization zones were also identified. In addition, an evaluation framework of urban–rural development was built to investigate how HANPP indexes respond to the regional unbalanced development and we identified the dominant driving factors of HANPP indexes from the perspectives of economy, population, urban expansion, and agricultural technology. The results showed that the total amount of HANPP was 2.68 PgC and gradually decreased from the southeast to the northwest of China, which represented 60.33% of the NPPpot. The grain-producing provinces had high HANPP (%NPPpot) values, while the southwestern provinces with complex terrain had low HANPP (%NPPpot) values. HANPP (gC/m<sup>2</sup> ), HANPPharv (gC/m<sup>2</sup> ), and HANPPluc (gC/m<sup>2</sup> ) were significantly different (*p* < 0.05) between the Huanyong Hu Line and the western–eastern part of China, while HANPPharv did not show significant differentiations in the coastal land and inland areas.

From the perspective of HANPP components, in China, the total amount and the average value of HANPPluc and HANPPharv were found to be 1.63 PgC and 216.97 gC/m<sup>2</sup> and 1.05 PgC and 126.53 gC/m<sup>2</sup> , respectively, and the harvest through cropland, livestock grazing, and forestry contributed 29.86%, 8.53%, and 0.91% to the total HANPP, respectively. Land use change had a bidirectional influence on HANPPluc (gC/m<sup>2</sup> ); the negative values were mostly concentrated in a small number of cultivated land patches in northwest China, which indicated that the biomass production under intensive artificial management could be higher than the potential production. Moreover, the North China Plain and northeast provinces had relatively high HANPPharv (gC/m<sup>2</sup> ) values, mainly provided by cultivated land. Different urban functions resulted in various ecological occupation structures. HANPPluc played the dominant role among 21 provinces (over 50% of HANPP), especially in the cities with rapid urban expansion and some coastal ones. With the production–consumption ratio of biomass based on HANPP, 17 provinces in the middle and southeastern coastal areas in China were defined as the consumption type. Specifically, the ratios in Shanghai (0.06) and Beijing (0.09) were under the fifth percentiles, and a great number of ecological resources were needed from surrounding regions. In addition, a universal positive correlation (*p* < 0.05) was found between the production–consumption ratio and the HANPPharv (% of HANPP) through agriculture, forestry, and grazing in 31 provinces in China.

In different dimensions of urban–rural development, each HANPP index was more likely to be affected by urban economy (UE) than rural economy (RE); the impact of rural population (RP) was much stronger than that of urban population (UP); and rural agricultural technology (RA) had more significant explanatory effects on each HANPP index than urban spatial expansion (US). Specifically, in the dimension of economic development, the added value of transportation, warehousing, and postal services (UE3) had a significant influence on most HANPP indexes. From the demographic perspective, rural

population (RP1) had significant and higher positive correlations with all the HANPPharv indexes than UP1 but no statistical correlation with HANPPluc. Meanwhile, agricultural technology indicators showed significant positive effects (*p* < 0.01) on and high correlations with HANPPharv. The higher regional average nighttime light intensity (US1), the proportion of the built-up area (US2), and the proportion of urban road area (US4) always corresponded with large HANPPluc values. Conversely, most HANPP indexes decreased as the proportion of urban green spaces increased (US3). Furthermore, the regional ecological pressure characterized by HANPP (%NPPpot) and HANPP (gC/m<sup>2</sup> ) mostly depended on the rural development index (*f*(*r*)) and HANPPluc and HANPPharv were mainly controlled by *f*(*u*) and *f*(*r*), respectively. A higher degree of urban–rural coordination (*D*) obviously improved HANPP (gC/m<sup>2</sup> ) and HANPPharv (gC/m<sup>2</sup> ) but had a relatively weaker effect on the variation in HANPPluc (gC/m<sup>2</sup> ).

**Supplementary Materials:** The following supporting information can be downloaded at https: //www.mdpi.com/article/10.3390/land12051062/s1: Figure S1: Conceptual model of human appropriation of net primary production (HANPP); Figure S2: Spatial distribution of NPPpot in China; Figure S3: Typical regionalization of China to detect the HANPP differentiation; Table S1: Moisture contents (%), harvest factors and recovery rates of typical crops; Table S2: Conversion coefficients from crop harvest to straw; Table S3: Regulation for assessing NPPpot by vegetation-approach model; Table S4: Urban–rural coordinated development indexes; Supplementary File S1: Estimation of HANPPharv; Supplementary File S2: Evaluation of urban–rural coordinated development.

**Author Contributions:** Conceptualization, T.Z. and J.P.; methodology, T.Z.; software, T.Z.; validation, T.Z. and J.P.; formal analysis, T.Z.; resources, J.P.; data curation, X.C.; writing—original draft preparation, T.Z.; writing—review and editing, J.P.; visualization, T.Z.; supervision, X.C.; funding acquisition, T.Z. and X.C. 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. 42001097) and the National Natural Science Foundation of China (Grant No. 41831284).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We gratefully acknowledge Helmut Haberl from the University of Natural Resources and Life Sciences for his advice on the research methods and the HANPP estimation model used in this manuscript.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
