**Assessment of Urban Ecological Resilience and Its Influencing Factors: A Case Study of the Beijing-Tianjin-Hebei Urban Agglomeration of China**

**Chenchen Shi 1,2,3, Xiaoping Zhu 4, Haowei Wu <sup>5</sup> and Zhihui Li 2,6,\***


**Abstract:** Climate change and rapid urbanization bring natural and anthropogenetic disturbance to the urban ecosystem, damaging the sustainability and resilience of cities. Evaluation of urban ecological resilience and an investigation of its impact mechanisms are of great importance to sustainable urban management. Therefore, taking the Beijing-Tianjin-Hebei Urban Agglomeration (BTHUA) region in China as a study area, this study builds an evaluation index to assess urban ecological resilience and its spatial patterns with the resilience surrogate of net primary production during 2000–2020. The evaluation index is constructed from two dimensions, including the sensitivity and adaptability of urban ecosystems, to capture the two key mechanisms of resilience, namely resistance and recovery. Resilience-influencing factors including biophysical and socio-economic variables are analyzed with the multiple linear regression model. The results show that during 2000–2020, the spatial pattern of urban ecological resilience in the BTHUA is characterized by high resilience in the northwest and relatively low resilience in the southeast. High resilience areas account for 40% of the whole region, mainly contributed by Zhangjiakou and Chengde city in Hebei Province, which is consistent with the function orientation of the BTH region in its coordinated development. Along with urbanization in this region, ecological resilience decreases with increased population and increases with GDP growth; this indicates that, although population expansion uses resources, causes pollution and reduces vegetation coverage, with economic growth and technological progress, the negative ecological impact could be mitigated, and the coordinated development of social economy and ecological environment could eventually be reached. Our findings are consistent with mainstream theories examining the ecological impact of socio-economic development such as the Environmental Kuznets Curve, Porter Hypothesis, and Ecological Modernization theories, and provide significant references for future urbanization, carbon neutrality, resilience building, and urban ecological management in China.

**Keywords:** urban resilience; land use/cover change; urbanization; carbon neutrality; Beijing-Tianjin-Hebei urban agglomeration

#### **1. Introduction**

Land use change induced by urbanization has greatly changed the physical conditions of urban ecosystems. Meanwhile, external risks faced by urban systems, such as climate

**Citation:** Shi, C.; Zhu, X.; Wu, H.; Li, Z. Assessment of Urban Ecological Resilience and Its Influencing Factors: A Case Study of the Beijing-Tianjin-Hebei Urban Agglomeration of China. *Land* **2022**, *11*, 921. https://doi.org/ 10.3390/land11060921

Academic Editor: Richard C. Smardon

Received: 16 May 2022 Accepted: 14 June 2022 Published: 16 June 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/).

change, require resilient urban ecosystems to provide a solid urban foundation to support urban production and livelihood [1]. Measures such as emission reduction are taken in urban systems to adapt to environmental change along with urbanization, for example, China is taking the lead in committing to achieving carbon neutrality in 2060 [2]. These measures would further change urban ecosystems with the goal of enhancing resilience. Unlike the sustainability concept, which is too comprehensive and hard to quantify in social-ecological systems, the resilience concept is derived from ecosystems and proves to be a good analytical tool for urban systems [3]. To quantitively assess urban ecological resilience and elucidate its driving mechanisms is of great importance to guide sustainable urbanization, mitigate climate change impact and achieve high-quality development.

The concept of ecological resilience is a measure of the persistence of an ecosystem and refers to the ability of a system to absorb natural and anthropogenic disturbances while still maintaining the interplay between population or ecosystem state variables and maintaining system function [4]. Ecological resilience in urban systems is a complex of socio-economic human activities and biophysical habitats [5]. An urban ecosystem is a specific type of ecosystem that integrates human activities into the biophysical sphere. In recent years, there has been much research on urban resilience that originates from socialecological resilience. The urban system itself is a more complex system than an ecosystem, as an ecosystem is only one of the subsystems within the complex urban system, which includes society, economy, culture, ecology and environment, etc. [6,7] Urban resilience refers to the capability of urban systems to remain intact or quickly recover to desired levels in the face of shocks or stresses [8,9]. As resilience in different scales embeds different connotations and mechanisms, this research focuses on urban resilience to theoretically conceptualize and empirically explore resilience in urban systems.

Utilization of the resilience concept must be based on the characterization and measurement of resilience. According to its concept, resilience can be measured by the rate of return of ecosystem state after change or disturbance, and the measurement of resilience is, in effect, the measurement of threshold-crossing [10]. Therefore, detecting an ecological threshold along disturbance gradients is essential to protecting the threshold from being crossed [11]. Previous research selects a key indicator or index to evaluate ecosystem resilience [12–15] or quantifies the economic value of resilience [16]. For urban resilience, a popular measurement is to build a comprehensive evaluation index with indicators representing the resilience of urban elements [17–19] or resilience process [20–22]. However, there are also scholars who believe that instead of seeking accurate metrics to measure resilience or trying to develop a general resilience index, it might be better to use surrogates or proxies [23,24]. As a matter of fact, evaluation of the previously mentioned resilience indicators is resilience surrogates in the sense that they demonstrate the impact factor of resilience rather than resilience (threshold-crossing) itself.

However, many resilience assessments tools, especially comprehensive evaluation indexes, tend to use state variables or cross-section data to measure resilience [17,25,26], while resilience is a process concept that encompasses two processes, namely resistance and recovery [27], making its measurement difficult as threshold-crossing often does not happen [24]. Therefore, in this study, drawing upon the assessment of economic resilience, where the change rate of GDP (also known as the sensitivity index) is used to measure resilience of the economic system [28], a resilience evaluation index is built to measure urban ecological resilience. Specifically, we use the change of net primary productivity (NPP) as the resilience surrogate. A concrete calculation method will be presented in the following method section. The surrogate selection of NPP is referenced from research that uses normalized difference vegetation index (NDVI) and GPP to quantify ecosystem resilience [29,30].

In terms of the impact factors of ecological resilience, as mentioned before, some studies use the impact factor of resilience to quantify system resilience and identified factors such as climate (e.g., precipitation, sunshine hour, temperature), hydrology (e.g., surface and groundwater resources), and land cover factors [31]. Among them, land cover is

an important and widely studied impact factor of ecosystem services and resilience [5,32,33]. Research shows urbanization increases land cover fragmentation, thus decreasing urban ecological resilience. In addition to biophysical factors, social-economic factors, such as GDP, population, land use, industrial structure, infrastructure, institutional arrangements, etc. also influence the resilience in urban systems [3,34–36]. An urban system, as a complex adaptive system, has a socio-economic subsystem that directly contributes to system resilience with its own urban functions, such as engaging in social production, providing employment, increasing productivity, and improving the livelihood of urban residents. Meanwhile, it influences system resilience indirectly through the interactions with the ecological-environment subsystem, which also directly impacts urban ecological resilience.

Therefore, this study sets out to measure urban ecological resilience with the empirical case of the Beijing-Tianjin-Hebei Urban Agglomeration and explore the impact factors including physical geographical factors (temperature, precipitation, elevation, slope, land cover type) and socio-economic factors (GDP, population, industrial structure, urbanization rate and carbon emission). Rationales for impact factor selection will be presented in the empirical results section. In doing this, we aim to empirically assess the ecological impacts of urbanization through the quantification of urban ecological resilience and explore the influencing mechanism of urban change on ecosystems. The following sections of this paper are structured as follows. The next section presents the overall methodology of this study with elaborations on the case study, research data and method. The empirical results are illustrated in Sections 3 and 4, with Section 3 presenting the resilience assessment result and analyzing the spatial difference of ecological resilience in the study region. Section 4 demonstrates the impact factor analysis and the correlation between ecological resilience and various impact factors. Section 5 concludes this paper by discussing the empirical, methodological, and theoretical implications and contributions of this study.

#### **2. Methodology**

#### *2.1. Study Area*

Previous empirical studies indicate that from 2006 to 2013, there was an N-shaped relationship between urbanization and ecological efficiency, and other studies show that from 2008 to 2017 [37], ecological efficiency slightly declined along with China's urbanization [38], while the current new-type urbanization in China (Refer to "The National New-type Urbanization Plan (2014–2020)" from National Development and Reform Commission of China: https://www.ndrc.gov.cn/xwdt/ztzl/xxczhjs/ghzc/201605/t20160 505\_971882.html?code=&state=123 (accessed on 30 April 2022)) that emphasizes the improvement of ecological environment, has yielded good ecological outcomes in pollution reduction and energy efficiency increase [39]. This study explores the ecological consequence of urbanization in China. Specifically, we selected the Beijing-Tianjin-Hebei Urban Agglomeration (BTHUA) region in China as the case study area to empirically assess ecological resilience and its influencing factors (Figure 1). The research was conducted on a 1 km × 1 km grid scale to observe fine-scale variation of urban ecological resilience change, and the research period is set from 2000–2020 to observe the rapid urbanization in China in the last twenty years and its ecological impacts.

The BTHUA region, also known as the Capital Economic Circle, is one of the major urban agglomerations in China, with high levels of urbanization and industrialization, though it suffers from ecological and environmental problems such as water resource shortage, air pollution, extensive land use and forest degradation, urban floods, and heatwave, etc. This region is semi-arid and semi-humid, has a temperate and warm temperate continental monsoon climate with an average temperature of 1 ◦C to 15 ◦C, abundant light, uneven spatial distribution of annual precipitation (about 300 mm to 750 mm from west to east), and the average annual evaporation is generally 900 mm to 1000 mm. The region is divided into plateau, mountain hills, basin and plain from northwest to southeast.

**Figure 1.** Location, land use and administrative division of the Beijing-Tianjin-Hebei Urban Agglomeration.

Currently, under the national policy of "Coordinated Development of the Beijing-Tianjin-Hebei Region", the overall positioning of the BTHUA region is "a world-class city cluster with the capital as the core, a new engine for economic growth driven by innovation and coordinated development, and a demonstration area for ecological restoration and environmental improvement" (Source: Outline of Coordinated Development of the Beijing-Tianjin-Hebei Region issued by State Council in 2015). Now, this region is still in the middle of industrialization and undergoing rapid urbanization and is also the main area of population inflow. In 2021, the permanent resident population rate in Beijing, Tianjin and Hebei were 87.5%, 84.88% and 61.14%, respectively (Data source: https://data.cnki.net/, accessed on 15 May 2022). In the future, this region will continue to host population inflow, with the ecological environment pressure, the contradiction of space resources utilization and the lack of regional infrastructure caused by urbanization expected to be more prominent [39].

#### *2.2. Data Source*

The data used for evaluating ecological resilience and driving factors of ecological resilience changes mainly include the Net Primary Production (NPP) data, meteorological data, topographic condition information, land use/land cover data, gridded GDP and population data, and carbon emission data. The NPP data for the years 2000–2020 with a spatial resolution of 500 m was derived from the MODIS MOD17A3HGFV06 product of vegetation NPP. The daily temperature and precipitation station monitoring data for the years 2000–2020 were obtained from the National Meteorological Information Center, and then the annual average temperature and precipitation were calculated and interpolated into a 1 km grid based on Kriging interpolation. The digital elevation model (DEM) and land use/land cover, population density and GDP per unit area with a spatial resolution of 1 km for the years 2000, 2005, 2010, 2015 and 2019 (using the population and GDP of 2019 to approximately represent that of 2020 due to data availability) were obtained from the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences. The CO2 emission data for the years 2000–2020 with a spatial resolution of 1 km were derived from the Open-source Data Inventory for Anthropogenic CO2 (ODIAC) datasets. All the data were interpolated or resampled into 1 km × 1 km resolution.

#### *2.3. Method*

In the quantification of ecological resilience, an evaluation index is built in this study. The resilience assessment results are classified with the natural breaks method as an objective classification method and in the impact factor analysis, a multiple linear regression model is adopted. For the resilience evaluation index, NPP is selected as the resilience surrogate. The evaluation model is constructed from two dimensions, including sensitivity and adaptability of urban ecosystems, to capture the two key mechanisms of resilience, namely resistance and recovery [27].

Sensitivity is the system's responsiveness to disturbance during normal operation. For a particular ecosystem, sensitivity is defined as the degree to which an ecosystem responds to disturbances such as climate change [40–42]. In this study, NPP is used to characterize ecosystem function. System sensitivity is represented by the interannual fluctuations of NPP from 2000 to 2020 to reflect the degree of dispersion of NPP from the average. The calculation formula is as follows:

$$S = \frac{\sum\_{i=1}^{n} |P\_i - \overline{P}|}{\overline{P}} \tag{1}$$

where *i* is the year (*n* = 21), *Pi* is the NPP value in year *i*. *P* is the mean value of NPP. *S* is the sensitivity index, which reflects the dispersion of NPP relative to the mean value over a specific time period.

Adaptability is the ability of a system to maintain and restore its structure in the face of disturbances [43]. Ecosystem adaptation, which is the self-regulation mechanism of an ecosystem, can be regarded as a measure to keep the system in a relatively stable state. In a certain period, the trend of variability of an ecosystem is used to measure its departure from homeostasis, which can be called ecosystem adaptation. If the variability trend decreases or does not change, the system tends to be relatively stable. Increased variability indicates that unstable systems adapt to changes and may indicate increased vulnerability [3,44]. Over a certain period, ecosystem adaptation, the self-regulating capacity of the system, can be expressed as the slope of a linear trend line fitting the interannual variability of NPP. In this study, adaptation is represented by the slope of a linear fitting trend line for the interannual variability of NPP from 2000 to 2020:

$$y = A\mathfrak{x} + B \tag{2}$$

$$A = \frac{n\sum xy - \left(\sum x\right)\left(\sum y\right)}{n\sum x^2 - \left(\sum x\right)^2} \tag{3}$$

where *x* is the time series from 2000 to 2020, *y* is the annual change rate of NPP, which is the annual absolute value change of NPP. It is the annual value of NPP minus the average value of NPP from 2000 to 2020. *A* is the changing trend of NPP variability, which is the regression slope of data sets *y* and *x* and indicates ecosystem adaptability; *B* is the intercept.

Resilience is a function of the characteristics, amplitude, and range of change rate, as well as the sensitivity and adaptability of the ecosystem [27,45]. Resilience is negatively correlated with sensitivity and positively correlated with adaptability. The levels of adaptability, sensitivity and resilience in each region are relative, and the sensitivity and adaptability calculated according to the preceding formula may not be in the same dimension. To analyze regional differences in resilience, the calculation results of sensitivity and adaptability should be standardized respectively before resilience calculations [46]. The resilience formula is as follows:

$$R = A - S \tag{4}$$

where *R* ecosystem resilience, *S* is sensitivity index, *A* is adaptability index.

#### **3. Spatial Pattern of Urban Ecological Resilience in the BTHUA Region**

The urban ecological resilience in the BTHUA region ranges from −0.23 to 0.74, with the mean value being 0.48, and is classified with the natural break method (Figure 2). The results show that during 2000–2020, urban ecological resilience in the BTHUA demonstrates discontinuous regional variance and a relatively small degree of dispersion, with the coefficient of variation being 0.1036. The overall spatial pattern is characterized by high resilience in the northwest and relatively low resilience in the southeast. Most of the regions are at middle-to-high ecological resilience levels, indicating high resilience in the Beijing-Tianjin-Hebei region over the past twenty years, which could be attributed to the improvement of the ecological environment in this region, especially since the coordinated development from 2015 during the 13th Five-Year-Plan period (2016–2020).

**Figure 2.** Spatial pattern of urban ecological resilience in the BTHUA region.

High resilience areas are concentratedly distributed and account for 40% of the whole region, mainly contributed by Zhangjiakou and Chengde city, as well as areas along the Taihang mountains in Baoding, Shijiazhuang, Xingtai and Handan cities in Hebei Province. This is consistent with the function orientation of Hebei Province as the "ecological environment support area" in the BTHUA region. While Beijing is mainly covered by high and middle-high resilience areas, with north and southwest Beijing being more resilient than the southeast. The good resilience performance in Beijing is also consistent with Beijing's development orientation in this region to serve as the "ecological restoration and environmental improvement demonstration area". Along the coastal lines in Qinhuangdao, Cangzhou and Tangshan, there are also areas with high resilience.

The middle-high resilient regions cover 27% of the whole region and are relatively concentrated in Beijing, Tianjin, coastal regions in Hebei Province and areas along the Taihang mountains. Areas with middle resilience account for 27% of the whole region and are relatively concentrated in south Hebei. Middle-low resilience regions account for 5% of the study area and are distributed in Tianjin and south Hebei. Low resilience areas, accounting for less than 1%, are scattered across the region, with continuous distributed

low resilience mainly in the southeast corner of Handan, northwest of Xingtai, southwest of Shijiazhuang, northeast of Cangzhou, central Baoding, north Tangshan in Hebei and middle and south Tianjin. This is attributed to the vegetation reduction along with land-use change induced by urbanization in these areas.

As for land-use type, areas with high resilience are mainly forest and grassland; middlehigh areas are grassland, waterbody and cultivated land; middle resilience areas are mainly grassland, cultivated and built-up land; and middle-low and low resilience areas are mainly cultivated and built-up land. The empirical results in this region show that cultivated and built-up land, disturbed by anthropogenic activities, demonstrates higher ecological risk and lower resilience. These lands with lower resilience are mainly located in the plain area of the BTHUA region between the Taihang mountains in the northwest and coastal lines in the southeast, covering southeast Beijing, central Tianjin, Shijiazhuang, central Tangshan, Qinhuangdao, Handan, Xingtai, Hengshui, Cangzhou, southeast Baoding, and Langfang cities in Hebei. These areas are also the highly urbanized areas in this region with dense populations, and they carry the burden of regional socio-economic development. However, the areas with high forest and grassland coverage are less influenced by human disturbance and exhibit high resilience levels. These areas are mainly in the Taihang mountains, Yanshan mountains, and Bashang grassland in mountain hill regions; forest and grassland can conserve water, clean air, regulate climate, and improve greening rate, forest coverage, air quality and biodiversity, generating positive ecological benefits.

Our findings show that high ecological resilience areas are mainly located in ecological land use-dominated areas, for example, mountain hills with forest landscapes. This indicates that during urbanization, ecological land is better preserved and subjected to less disturbance, and low resilience is mainly in built-up and cultivated land of plain areas. Though cultivated land has certain ecological regulating functions, contributing to ecological resilience and anthropogenic use of cultivated land, for example, pesticide application and mechanized production as well as intensive and irrational use, can cause ecological damage and reduce the resilience of ecosystems. The ecological benefits of built-up land are limited. If the development mode is not reasonable, for example, extensive built-up land use, which will result in the disorderly spread and expansion of built-up land, it will increase the negative impact on the ecosystem. Therefore, during urbanization, ecological protection and restoration projects should be implemented to conserve ecological land by protecting forests, grasslands, wetlands, and other key ecological resources, and to improve ecosystem function and stability.

For now, China is implementing 'Planning for Major Ecological Protection and Restoration Projects in the Northern Sand Control Belt (2021–2035)' in the study region (Source: http://www.gov.cn/zhengce/zhengceku/2022-01/14/content\_5668161.htm, accessed on 15 June 2022), to promote comprehensive ecological management in Zhangjiakou, Chengdu, Yanshan and Taihang mountains and Baiyangdian in Xiongan New District in Baoding. According to this plan, by 2035, the quality and stability of natural ecosystems such as forests, grasslands, rivers and lakes, wetlands and deserts will be significantly improved to enhance ecological service and resilience. This will help achieve carbon peak and carbon neutrality and build an ecological security barrier in northern China, laying a solid ecological foundation for realizing the goal of building a beautiful China. While this study consolidated the necessity of implementing such a plan, its effects on ecological resilience-building deserve further research attention.

#### **4. Impact Factors of Urban Ecological Resilience in the BTHUA Region**

To analyze the correlation between urban ecological resilience and its impact factors (Table 1), a multiple linear regression model is built. Impact factors are selected based on previous studies and data availability. In this study, we try to include both physical geographical and socio-economic factors. For biophysical factors, climate factors including temperature and precipitation play important roles in the formation of ecosystem structure and function. They directly impact the growth of vegetation and thus tend to impact ecological resilience. In the study region, the annual mean temperature is decreasing from the southeast to the northwest (Figure 3a), while precipitation is concentrated in east Hebei plain and south Hebei (Figure 3b). Vegetation coverage is also significantly affected by topography and correlated with elevation and slope. The topography of the study region is high in the northwest and low in the southeast, tilting from the northwest to the southeast (Figure 3c). The slope rises from the southeast to the Taihang Mountains and drops to the northwest after the Taihang Mountains (Figure 3d).


**Table 1.** Descriptive statistics of resilience impact factors.

Socio-economic factors, including population, GDP, carbon emission, and built-up land, mark the impact of human activities on the urban ecosystem. As an important approach to regional sustainable development, urban resilience is a coordinated and optimized combination form of multiple factors such as economy, population, land and market. Among them, the built-up area is the spatial carrier of regional resource convergence and diffusion, which influences urban resilience through multiple forms of resource combination. Reasonable population spatial pattern and economic development pattern are the main ways for cities to deal with natural disasters and urban diseases. In the study region, the selected socio-economic indicators all show the characteristics of clustering towards the urban center (Figure 3e–h).

The multiple linear regression results show that from the biophysical side, ecological resilience in the BTHUA region decreases with the increase of temperature and precipitation, showing that ecosystems in the high temperature and heavy rainfall areas in this region are less stable. Resilience decreases by 0.0389 for every 10-degree Celsius increase in temperature, and by 0.0089 for every 1000 mm precipitation increase, and the results are significant with a P value less than 0.001 (Table 2). Elevation and slope determine the spatial distribution of vegetation types to a large extent, and the degree of ecological resilience increases with the increase of average elevation in the region, with a 0.0093 increase in ecological resilience for every 1000 m elevation increase. In addition, ecological resilience in the study region increases by 0.0303 with every 10-degree slope increase. The system resilience is significantly correlated with topography (with P values both less than 0.001).

However, from the socio-economic side, urban ecological resilience is positively related to GDP and negative related to population and built-up land. For every extra 100 people per unit km2 area, resilience decreases by 0.0187%, while when built-up land expands, resilience tends to decrease significantly. For every 1 million increases in GDP per unit km2 area, resilience increases by 0.00187%. It shows that the level of economic development has a significant positive impact on urban resilience. This can be attributed to the high overall economic development level in the BTHUA region; resources could be allocated to technological (new technologies that increase resources efficiency and reduce pollution) and institutional (stringent environmental protection policies and enforcements) advancement, generating benign ecological and environmental effects.

**Figure 3.** Spatial pattern of resilience impact factors. (**a**) temperature; (**b**) precipitation; (**c**) elevation; (**d**) slope; (**e**) population; (**f**) GDP; (**g**) carbon emission; (**h**) built-up land.


**Table 2.** Regression results of the resilience impact factor analysis.

Along with urbanization in this region, ecological resilience decreases with population boost and built-up land expansion and increases with GDP growth; this indicates that, though urban expansion takes up resources, causes pollution and reduces ecological resilience, with technological and institutional progress generated by economic growth, the negative environmental impact could be mitigated and we could eventually reach the coordinated development of social economy and ecological environment. In addition, as ecological resilience is calculated from NPP change, it shows the destruction of vegetation by population expansion and urbanization. However, such damage could be mitigated by technological advances and institutional change along with economic development. The results in this study are consistent with mainstream environmental theories such as the Porter Hypothesis (the Porter Hypothesis holds that strict environmental standards stimulate innovation and improve environmental quality to offset the conflict between economic development and environmental protection [47,48]. Interpretations of the Porter Hypothesis show that stringent environmental regulations boost the environmental services sector, induce technological innovation, and provide some firms with an early mover advantage; large firms benefit more than small firms because of their lower compliance costs [49]), the Environmental Kuznets Curve (the Environmental Kuznets Curve (EKC) holds that economic modernization will reconcile the conflict between economic and environmental interests [50]. EKC is about the relationship between environmental degradation and growth at different levels of economic development. At low levels of income, people pay attention to growth and worry little about the environment, so growth produces degradation; at higher incomes, people pay more attention to the environment compared to growth, therefore degradation declines with increasing incomes, making the relationship between growth and degradation takes the shape of an upside-down U) and Ecological Modernization (Ecological Modernization is a model of environmental governance originating from Europe in the 1980s, which holds that the capitalist economy could solve this problem through technological advances that mitigate the trade-off between economy and environment, with the participation of civil society actors in addition to state and market actors in a democratic setting [51–53]. According to ecological modernization theory, the conflict between industrial development and environmental protection could be resolved through environmental innovation without fundamental changes to production and consumption) theories.

#### **5. Discussion and Conclusions**

In this study, with the BTHUA as the case study area, urban ecological resilience and its spatial pattern are assessed to quantify the ecological impacts of urbanization. The results of the ecological resilience assessment show high ecological resilience areas are mainly in ecological land with forest and grassland as the main landscape. Therefore, it is important to protect key ecological resources to improve resilience through the provision of ecosystem services [54]. The current "Planning for major ecological protection and restoration projects in the northern sand control belt" in this region should be duly implemented. In addition, as the results indicate the difference in spatial resilience pattern, with the overall spatial pattern being characterized by high resilience in the northwest and relatively low resilience in the southeast, future development in the BTHUA region should adhere to the spatial orientation of regional coordinated development, which position Hebei as the ecological support area and Beijing as the ecological restoration demonstration area. This also indicates that the use of geospatial information on resilience could provide effective management references [55].

The resilience impact factor analysis integrates both natural-physical and socio-economic indicators. From the physical side, climate and topographical factors are found to have a significant impact on urban ecological resilience. From the human side, the economy, population and land use are also correlated with urban resilience. The empirical case in the BTHUA region indicates that along with urbanization, ecological resilience is negatively related to population and built-up land expansion and positively related to GDP growth. To explore the reasons for such an impact, although urban expansion utilizes resources, causes pollution and reduces resilience, economic growth and technological and institutional advancement could mitigate negative ecological impact. The findings are consistent with the Environmental Kuznets Curve, the Porter Hypothesis, and Ecological Modernization theory and the like, which also examine the environmental/ecological impact of socioeconomic development. With further economic development in this area, the environmental and ecological burden could be eased with advanced technology and stringent regulation. The carbon neutrality scheme, for example, is a case in point, as it promotes technological (e.g., energy efficiency increase, electrification, renewable energy technology, etc.) and institutional changes (e.g., carbon pricing, carbon tax, carbon market, green finance, etc.) to conserve resources and reduce pollution.

The empirical results in this study provide important policy references for future urbanization, carbon neutrality, resilience building and urban ecological management in this region in specific and in China in general. With the findings on the ecological impact of urbanization in this study, future urbanization in China should integrate the concept and principles of ecological civilization into the whole process of urbanization and pursue a new type of urbanization characterized by an intensive, smart, green and low-carbon growth pattern (both ecological civilization and new-type urbanization are political discourses in China. Ecological civilization was first introduced into Chinese ideology in 2007 at the 17th Congress of the Communist Party, endorsed by President Xi in 2013 in environmental law and policy-making and written into the constitution in 2018. It aims at solving ecological and environmental problems with technological innovation as well as improved governance institutions. New-type urbanization is a guideline put forward in the report of the 18th National Congress of the Communist Party. New-type urbanization is characterized by urban-rural integration, city-industry interaction, conservation and intensification, ecological livability and harmonious development. It is characterized by coordinated development and mutual promotion of large, medium and small cities, small towns and new-type rural communities). With the carbon neutrality program expecting to drastically change the industrial structure as well as urban-rural relations, it is very important for the future urban development mode to find a realistic path of carbon neutrality from the coordination of urban and rural ecology. Carbon neutrality itself would contribute to the strengthening of urban resilience, and urban resilience governance should focus on all urban subsystems, including economic, social, natural and built-up environment, etc., as system elements interacting with each other and contributing to the emergence of the complex urban system.

To explore the change mechanism of urban ecological resilience from the perspective of resilience characteristics and the temporal and spatial differences of urban ability to prevent and defuse ecological risks, we clarify that the important task of ecological governance in urban zones is a favorable way to realize ecological risk prevention and control in resilient cities. In previous studies, resilience assessment can be categorized into two groups: compound index systems that capture as many resilience characteristics as possible and the use of a single indicator as a resilience surrogate. Though compound index systems excel in capturing complex resilience characteristics, they are generally based on static

indicators and describe system states. However, for single indicator evaluation, thresholdcrossing can be measured. In this study, the change rate of NPP is measured to characterize threshold-crossing. In addition, the resilience surrogate methods that use a single indicator require fewer data and are more operable in empirical studies. The resilience evaluation index in this study can easily be applied in other cases with available data.

The urban ecological resilience index we built is a methodological contribution as well as a new conceptualization of ecological resilience, which uses the NPP change to quantify the sensitivity and adaptability of urban ecosystems to shocks. With the urban system as a special type of ecosystem that integrates human activities with natural habitat, and with its complex adaptive system nature, the analysis in this study integrates urban elements in different dimensions and analyzes their impact on urban ecosystem resilience. The application of complex adaptive system theory in ecosystem resilience research is also one of the theoretical contributions of this study to further bridge social–physical complex networks, as previously done by Cavallaro et. al. [56]. One empirical advantage of this study is that it was conducted on a grid-scale, while most other urban resilience research is done on a prefectural city level. Further, as we were constrained by the availability of grid-scale socio-economic data, which generally depend on administrative statistical units, some key factors (for example, industrial structure) might be omitted in this research. This constitutes our future research agenda, provided that fine-scale data are available or when the research is conducted upscale.

**Author Contributions:** Conceptualization, C.S. and Z.L.; methodology, Z.L. and C.S.; formal analysis, Z.L. and C.S.; data curation, X.Z. and H.W.; writing—original draft preparation, C.S. and X.Z.; writing—review and editing, C.S. and Z.L.; visualization, Z.L. and H.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 number 72104149; Academy of Metropolis Economic and Social Development at Capital University of Economics and Business, grant number ZSXM2021003; Capital University of Economics and Business: The Fundamental Research Funds for Beijing Universities, grant number XRZ2021048 and QNTD202009.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All data and materials are available upon request.

**Acknowledgments:** We would like to thank the reviewers for their thoughtful comments that helped improve the quality of this work.

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

#### **References**


### *Article* **Spatial Distribution of Optimal Plant Cover and Its Influencing Factors for** *Populus simonii* **Carr. on the Bashang Plateau, China**

**Yu Zhang 1,2,3, Wei Li 4,5,\*, Shaodan Li 1,2,3, Baoni Xie 4,5, Fangzhong Shi <sup>6</sup> and Jianxia Zhao <sup>7</sup>**

	- Shijiazhuang 050024, China <sup>4</sup> School of Land Science and Space Planning, Hebei GEO University, Shijiazhuang 050031, China; xbn-feya@nwafu.edu.cn
	- <sup>5</sup> International Science and Technology Cooperation Base of Hebei Province: Hebei International Joint Research Center for Remote Sensing of Agricultural Drought Monitoring, Hebei GEO University, Shijiazhuang 050031, China
	- <sup>6</sup> School of Natural Resources, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China; fzshi@mail.bnu.edu.cn
	- <sup>7</sup> College of Geography and Land Engineering, Yuxi Normal University, Yuxi 653100, China; zjx@yxnu.edu.cn
	- **\*** Correspondence: weil87land@hgu.edu.cn

**Abstract:** The Bashang Plateau is the core zone of the agro-pastoral ecotone in northern China and represents an ecological barrier for preventing the invasion of wind-blown sand in the Beijing–Tianjin– Hebei region. Increasing plant cover to control soil erosion is an effective measure to address land degradation; however, plant cover is different from climatic conditions. In this study, we determined the optimal spatial distribution of *Populus simonii* Carr., which is a widely planted species used for revegetation on the Bashang Plateau. A modified Biome-BGC model was used to simulate the dynamics of the net primary productivity (NPP), actual evapotranspiration (AET), and leaf-area index (LAI). The model was validated using field-observed tree-ring and MODIS AET and NPP data. The dynamics of AET, NPP and LAI for *P. simonii* at 122 representative sites in the study area for the period 1980–2019 were simulated by the validated model. The results showed that the spatial distributions of mean AET, NPP, and LAI generally decreased from southeast to northwest. The ranges of optimal plant cover in terms of maximum LAI for *P. simonii* were 3.3 in the Fengning– Weichang area, 1.9 in the Shangyi–Zhangbei–Guyuan area and 1.3 in the Kangbao area. Mean annual precipitation (MAP), elevation, soil texture and mean annual temperature were the main factors influencing the distribution of AET, NPP and LAI. As the MAP decreased, the correlations between AET, NPP, LAI and precipitation gradually decreased. In different subregions, the factors influencing optimal-plant-cover distribution varied significantly. These quantitative findings provide the optimal plant cover for the dominant tree in different subregions and provide useful information for land degradation management on the Bashang Plateau.

**Keywords:** Bashang Plateau; *Populus simonii* forest; optimal plant cover; leaf-area index; Biome-BGC model

#### **1. Introduction**

Afforestation of degraded land is one of the major strategies for improving the ecological conditions at the regional levels by controlling soil erosion [1], enhancing soil-vegetation carbon sinks, maintaining biodiversity [2,3] and regulating hydroclimatic conditions [4]. However, immense afforestation may consume excess water and cause a severe soil water deficit [5,6], which in turn will negatively affect the ecological functions of vegetation in semiarid and arid regions. Attaining an optimal plant cover is vital to afforestation

**Citation:** Zhang, Y.; Li, W.; Li, S.; Xie, B.; Shi, F.; Zhao, J. Spatial Distribution of Optimal Plant Cover and Its Influencing Factors for *Populus simonii* Carr. on the Bashang Plateau, China. *Land* **2022**, *11*, 890. https://doi.org/10.3390/ land11060890

Academic Editor: Manuel Pulido Fernádez

Received: 18 May 2022 Accepted: 9 June 2022 Published: 11 June 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/).

projects for balancing soil water consumption, fostering plant growth and maintaining the ecosystem health of restored ecosystems [7,8].

An optimal plant cover is one that reaches the maximum leaf-area index (LAI) in the given climate–soil conditions and vegetation types without causing a permanent soil water deficit [9]. Several studies on optimal plant cover for different vegetation types have been conducted at varying scales ranging from plot and field scales to small watershed and regional scales using field experiments and numerical simulations in semiarid and subhumid areas [9–12]. Fu et al. combined field experimental data with the SHAW model to calculate the optimal cover corresponding to a maximum LAI of 1.27 for *C. korshinkii* and 0.70 for *S. psammophila* in the northern region of China's Loess Plateau [9]. Mo et al. used Eagleson's ecohydrological optimality method to determine the optimal canopy cover in Horqin Sands of China [8]. In the Northeast China Transect, Cong et al. determined that the optimal canopy cover for forests equates to 0.822 [13]. However, most field experiments used in research on optimal plant cover last only 1–3 years and cannot represent variations in long-term climate characteristics. An eco-physiological process model can estimate the dynamics of plant coverage using long-term climate records and has proven to be an effective tool to estimate optimal plant cover. Jia et al. estimated the spatial distributions of optimal plant cover for three exotic plant species across China's Loess Plateau using the Biome-BGC model [12]. The Biome-BGC model has been proven to simulate the optimal plant cover for typical species in arid and semi-arid regions.

The Bashang Plateau is one of the most vulnerable and sensitive portions of northern China and is a major source of dust and sand storms that affect Beijing and its surrounding areas (the Beijing–Tianjin–Hebei region) [14,15]. To prevent further land degradation and restore degraded lands is especially important in this area. A series of afforestation projects have been implemented by China's government over the last four decades [16], e.g., the Grain-for-Green Program, the Three-North Shelterbelt Program and the Beijing–Tianjin Sand Source Restoration Project, to address land degradation and improve the regional ecological environment. Since then, the vegetation cover has dramatically increased, and wind and sand erosion have been effectively controlled [17,18]. The excessive introduction of fast-growing plant species with high planting density caused a very large consumption of local water resources and further aggravated land degradation. *Populus simonii* Carr. is the dominant revegetation species used [19]; however, in recent decades, there has been a great deal of degeneration and death among *P. simonii* trees. Nearly 80% of the planted forest is degraded, whereas approximately one-third of it is dead [20]. Drought is associated with tree mortality at the regional level [21]. Thus, determining the optimal plant cover for *P. simonii* under Bashang Plateau conditions is crucial for accelerating recovery and combating land degradation.

The optimal plant cover for a given vegetation type depends on many factors and is sensitive to changes in the local environment and climate [12]. There were obvious subregional differences in climatic and soil conditions in the Bashang Plateau. However, previous studies have only determined the optimal plant cover and influencing factors at the regional scale. Few studies analyze the main influencing factors in different subregions. We hypothesized that the spatial distribution of optimal plant cover for *P. simonii* showed heterogeneity and the main influencing factors varied in different subregions. Thus, the objectives of this research were to (1) determine the spatial distribution of optimal plant cover for *P. simonii* in the study area based on long-term variations in net primary productivity (NPP), actual evapotranspiration (AET), and LAI for *P. simonii* using a modified Biome-BGC model and (2) determine the main factors that influence the spatial distribution of the optimal plant cover in three subregions with different climate, soil and topography conditions. This information should be useful in policy decision support for addressing land degradation through spatially explicit optimal plant cover for the major afforestation species on the Bashang Plateau of China.

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

#### *2.1. Study Area*

The Bashang Plateau, which is located between 114◦35 ~116◦45 E and 41◦00 ~42◦20 N, encompasses six counties (Guyuan, Zhangbei, Kangbao, Shangyi, Weichang, and Fengning) within the northern part of Hebei Province, China (Figure 1). It covers a total area of 18,202 km<sup>2</sup> and has an elevation range of 831~2215 m. The Bashang Plateau is characterized by a continental monsoon climate that is cold, windy, and prone to drought. The mean annual precipitation (MAP) varies from 330 mm in the northwest to 460 mm in the southeast. The mean annual temperature (MAT) ranges from 1.5 to 5 ◦C, and the annual potential evaporation is approximately 880~1000 mm. The soil textures are mainly composed of clayey and sandy soils. The vegetation types in this area are warm temperate deciduous broadleaf forest and temperate grassland [22]. The Bashang Plateau is characterized by vulnerable ecosystems, and various non-native plant species (including *Populus* L., *Ulmus pumila* L., *Pinus tabuliformis* Carr., *Pinus sylvestris var. mongolica* Litv., *Larix principis-rupprechtii* Mayr., *Betula platyphylla* Suk., *Caragana microphylla* Lam. and *Hippophae rhamnoides* L.) have been introduced in this area in the efforts to restore vegetation. *P. simonii*, a fast-growing tree that was cultivated in the 1980s to promote ecological protection, was investigated in this study. Based on the geographical and ecological characteristics, we divided the Bashang Plateau into three subregions (Table 1).

**Figure 1.** The location of the study area and the distributions of the 34 model evaluation sites, 122 simulation sites and precipitation contours in the region.



#### *2.2. Biome-BGC Model Description*

ĉ

The Biome-BGC model is a typical process-based model designed to simulate biogeochemical element cycling. It can predict the dynamics of water, carbon, nitrogen and energy in terrestrial ecosystems with general stand information and daily meteorological data [23,24]. The model provides complete parameter settings for different biomes (i.e., evergreen needle leaf, evergreen broadleaf, deciduous broadleaf, shrub, and grass) [23]. Leaf area determines the potential forest productivity as well as the rainfall partitioning and evaporation from soil. The LAI is a vegetation state variable that is updated every day on the basis of the estimated leaf carbon [25]. NPP is partitioned into biomass compartments according to dynamic allocation patterns based on nitrogen limitations [26]. The effect of water stress on the rate of photosynthesis and respiration is regulated through stomatal adjustments to simulate plant growth (Figure 2). A physically based equation (i.e., the one-dimensional Richards equation with a root water extraction term) was implemented in water-cycling routines to simulate soil water movement under limited water conditions by Huang et al. [27]. The inner process of the Biome-BGC model and the rationale have been described in depth by White et al. and Thornton et al. [23,24]. The original model is not suitable for coarse soil, and this study has therefore used the modified Biome-BGC model.

**Figure 2.** Detailed flow chart of the Biome-BGC model [28].

#### *2.3. Parameterization of the modified Biome-BGC Model*

#### 2.3.1. Meteorological Parameters

Measured values of the daily maximum and minimum air temperature (◦C), precipitation (mm) and relative humidity (%) for 12 meteorological stations inside and around the study area were provided by the national ground stations of the China Meteorological Administration. The data were interpolated using thin-plate smoothing splines to produce maps at a spatial resolution of 1 km for 1980–2019. The MT-CLIM model was used to generate other necessary meteorological parameters based on the temperature, precipitation, and humidity data [24].

#### 2.3.2. Soil Hydraulic Parameters

We investigated and measured the soil samples in each subregion under different topography, soil and vegetation conditions for the model evaluation. The latitude, longitude and elevation of each site were recorded by a GPS receiver, while the site slope and aspect were determined using a geological compass. At each site, several soil points were randomly selected to collect soil samples. The bulk density (*BD*), saturated SWC (*θs*) and saturated hydraulic conductivity (*Ks*) at each site were measured using undisturbed soil cores. Disturbed soil samples were collected using a soil-drilling sampler at depths of 0–0.2 m at each point. Various soil hydraulic parameters are necessary for the modified Biome-BGC model, and they include the residual soil water content (*θr*), *θs*, *Ks* and the van Genuchten model parameters (α, *n* and *m*). *Ks* was determined using the constant-head method, and *BD* was determined from the volume–dry-mass relationship [29] for each core sample. The hydraulic parameters *θr*, *α* and *n* were determined using the estimated soil–water-retention curves with the RETC program for each site [30]. A standard soil suction of 33 kPa was used to determine the soil water content at field capacity. The soil characteristic parameters at the simulation sites were based on previously published data from studies in this area [31–33] and the "National 1:1 million Digital Soil Map" developed by the Nanjing Soil Research Institute of the Chinese Academy of Sciences and Soil Environment Department of the Ministry of Agriculture of China. Digital-elevationmodel (DEM) data, ASTER GDEM 30 m, were acquired from the Geospatial Data Cloud (http://www.gscloud.cn (accessed on 16 May 2022)).

#### 2.3.3. Ecophysiological Parameters

The ecological parameters for *P. simonii* are summarized in Table 2. These parameters were derived from field observations for *P. simonii* by Kang et al. [34] and the published literature of the genera. Based on the study by Liu et al. [35], the distribution of roots for *P. simonii* was determined. The maximum root depth was assumed to be constant at 2 m.

**Table 2.** The ecophysiological parameters for *P. simonii* forest used in the Biome-BGC model.




Note: \* is the default parameter value of the model, \*\* is the optimized parameter value, and the remaining data are from the literature [34]. FC is the field water capacity, and WP is the wilting point.

#### *2.4. Validation of the Biome-BGC Model*

Tree-ring chronology can record long-term tree growth under natural conditions and can thus be used as an indicator of NPP for evaluating forest productivity in response to climate change [36]. A model simulation of the interannual variations in the NPP of *P. simonii* forests was tested using the tree-ring-width index (RWI) during 1982–2017 for *P. simonii* at the Zhangbei Ertai site derived from the published data of Liu et al. [37]. To test the performance of the modified Biome-BGC model at the spatial scales of NPP and AET, the simulated annual mean values of NPP and AET from 2001 to 2019 were compared with remote-sensing observation data in this study.

MODIS satellite images (MOD 17A3 H v006 and MOD16 A3) with a 500 m resolution were used to obtain the AET and NPP of the study area for the period from 2001 to 2019. Image processing was conducted using ENVI 5.1 software. The modified Biome-BGC model was used to simulate the past NPP and AET and to compare them with the MODIS AET and NPP. The accuracy of the modified Biome-BGC model was assessed by the coefficient of determination (*R2*), mean difference (*MD*) and root mean square error (*RMSE*).

#### *2.5. Determination of Optimal Plant Cover*

The modified Biome-BGC model was used to simulate annual variations in AET, NPP, and LAI for *P. simonii* across 122 representative sites on the Bashang Plateau from 1980 to 2019. Due to seasonal changes in LAI, the optimal plant cover equals the average of the maximum LAI during the simulation.

#### *2.6. Statistical Analyses*

Basic statistical analyses were used to analyze the simulations of AET, NPP and LAI for each subregion and the Bashang Plateau. Pearson correlation analyses were used to assess the relationships between AET, NPP and LAI and climate (precipitation and temperature), soil texture (sand, silt and clay contents) and elevation for *P. simonii*. A stepwise regression analysis was used to identify the major variables that accurately predicted the maximum LAI for *P. simonii*. All statistical analyses were conducted using SPSS 22.0. Geostatistical analyses and maps were produced by ArcGIS 10.1.

#### **3. Results**

#### *3.1. Model Evaluation*

The simulated NPP values for the sample plot in Zhangbei County were consistent with the interannual variation trend in the RWI during the period of 1982–2017, except in 1998 (r = 0.40, *p <* 0.05) (Figure 3a), and agreed with the MODIS data during the period of 2001–2017 (r = 0.90, *p <* 0.01). Regardless, the values simulated by the modified Biome-BGC model were consistent with the MODIS NPP and ET data during the period of 2001–2019, and the correlations were significant (r > 0.77, *p <* 0.01) (Figure 3b,c). The results indicated

that the temporal dynamics of NPP and AET simulated by the modified Biome-BGC model were acceptable for *P. simonii*.

**Figure 3.** (**a**) Comparisons between the *P. simonii* forest NPP based on the Biome-BGC model and MODIS data for NPP and the tree-ring-width index (RWI) at the Zhangbei Ertai site during the period of 1982–2017; (**b**) comparisons between the *P. simonii* forest NPP based on the Biome-BGC model and MODIS data and (**c**) AET based on the Biome-BGC model and MODIS data at the Zhangbei Ertai sites during the period of 2001–2019.

The comparisons of the annual average MODIS NPP and AET data with those simulated by the modified Biome-BGC model during the period of 2001–2019 for *P. simonii* at the regional scale are shown in Figure 4. In the whole Bashang Plateau, the correlations between the simulated and observed NPP and AET were >0.80, and the *MD* and *RMSE* values were low. For the three subregions, the correlations between the simulated and observed NPP and AET ranged from 0.72 to 0.94, and the *MD* and *RMSE* values were 22.53~35.51 g C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> and 26.29~57.34 g C m−<sup>2</sup> <sup>a</sup>−<sup>1</sup> for NPP, and −31.35~33.18 mm and 15.28~60.07 mm for AET, respectively (Table 3). Each subregion showed a low degree of simulated NPP and AET dispersion. The results thus indicated that the spatial dynamics of NPP and AET simulated by the modified Biome-BGC model were acceptable for *P. simonii*.

**Figure 4.** Comparisons of the Biome-BGC modeled (**a**) NPP, (**b**) AET and MODIS data for NPP and AET from 2001 to 2019 on the Bashang Plateau.


**Table 3.** Assessment index between the modified Biome-BGC model and MODIS data for NPP and AET in the study area.

#### *3.2. Spatial Distributions of AET and NPP*

The AET and NPP were simulated by the modified Biome-BGC model using 122 representative sites across the Bashang Plateau from 1980 to 2019. The spatial distributions of the mean values of AET and NPP of *P. simonii* were then mapped by kriging interpolation (Figure 5a,b). The AET generally increased from the northwest to the southeast. The annual average value of the estimated AET was 398.7 mm, with a range of 283.6~532.7 mm and a spatial variation coefficient of 13.7% (Figure 5a). There were significant differences in the AET for *P. simonii* in different subregions, of which zone III had the highest AET (447.8 mm), and zone I had the lowest AET (323.1 mm). The AET for zone II was 408.7 mm (Table 4). The spatial pattern of the mean AET was consistent with that of the MAP.

**Figure 5.** Spatial distribution of the mean (**a**) AET and (**b**) NPP and (**c**) maximum LAI for *P. simonii* forest on the Bashang Plateau during 1980–2019.

Consistent with the mean AET, the mean NPP increased from northwest to southeast. The annual average value of the estimated NPP was 360.5gCm−<sup>2</sup> yr−1, with a range of 191.3~542.7 g C m−<sup>2</sup> yr−<sup>1</sup> and a spatial variation coefficient of 23.7% (Figure 5b). There were significant differences in NPP for *P. simonii* in different subregions. The highest mean NPP was 435.7gCm−<sup>2</sup> yr−<sup>1</sup> in zone III, the mean NPP for zone II was 373.0gCm−<sup>2</sup> yr<sup>−</sup>1, and the lowest mean NPP was 251.3gCm−<sup>2</sup> yr−<sup>1</sup> for zone I (Table 4).

**Table 4.** Mean NPP, AET and maximum LAI values simulated by the modified Biome-BGC model for *P. simonii* forest in three zones of the Bashang Plateau.


Note: The values of each property in the table are the average values ± standard deviation for the study sites in three zones of the Bashang Plateau.

#### *3.3. Spatial Distribution of the LAI and Optimal Plant Cover*

The spatial patterns of the mean maximum LAI of *P. simonii* were also mapped by kriging interpolation (Figure 5c). The optimal plant cover was calculated as the mean maximum LAI of *P. simonii* during the study period. In accordance with the MAP distribution, the optimal plant cover generally increased from northwest to southeast, with a range of 1.0~4.1 m2 m−2. The overall mean optimal-plant-cover value was 2.08 m2 m−<sup>2</sup> and the spatial variation coefficient was 38.7%. The optimal plant cover varied by subregion (Table 4). It was less than 1.5 m<sup>2</sup> m−<sup>2</sup> in zone I with a MAP < 350 mm, and the mean value was LAI = 1.3 m<sup>2</sup> m−2. The optimal plant cover ranged from 1.5 to 3.0 m2 m−<sup>2</sup> in zone II with a MAP = 350~450 mm, and the mean value was LAI = 1.9 m2 m<sup>−</sup>2. The optimal plant cover was higher than 3.0 m<sup>2</sup> m−<sup>2</sup> in zone III with a MAP > 420 mm, and the mean value was LAI = 3.3 m2 m<sup>−</sup>2.

#### *3.4. Mean AET, NPP and LAI Distribution Factors*

Highly significant positive relationships were found between the mean AET and MAP for *P. simonii* on the whole Bashang Plateau and in its three subregions (r > 0.965, *p* < 0.01) (Table 5), suggesting that the MAP was a main influencing factor of the spatial pattern of the mean AET on the Bashang Plateau. Furthermore, the mean AET was negatively correlated with elevation (r = −0.748, *p* < 0.01), sand content (r = −0.714, *p* < 0.01) and mean annual temperature (MAT, r = −0.417, *p* < 0.01) on the Bashang Plateau. There were significant differences in distribution factors in different subregions. The mean AET was significantly positively correlated with the clay (r = 0.635, *p* < 0.01) and silt contents (r = 0.689, *p* < 0.01) and negatively with the sand content (r = −0.659, *p* < 0.01) and elevation (r = −0.390, *p* < 0.01) in zone I (n = 28). The mean AET was significantly correlated with elevation (r = 0.493, *p* < 0.01) but negatively correlated with MAT (r = −0.606, *p* < 0.01) in zone II (n = 68). A negative correlation was found between the mean AET and sand content (r = −0.427, *p* < 0.01) in zone III (n = 26). The influencing factors of the spatial variation in the NPP for *P. simonii* on the Bashang Plateau and its subregions were similar to those of AET except in zone III. NPP was significantly positively correlated with MAT (r = 0.413, *p* < 0.01) in zone III (n = 26). The influencing factors of the spatial variation in the maximum LAI in the three subregions were similar to those of the NPP. However, the maximum LAI was significantly positively correlated with silt content (r = 0.230, *p* < 0.01) and was negatively correlated with sand content (r = −0.550, *p* < 0.01), elevation (r = −0.531, *p* < 0.01) and MAT (r = −0.491, *p* < 0.01) on the Bashang Plateau (n = 122).

**Table 5.** Pearson correlation coefficients between mean ET, mean NPP, maximum LAI, climate variables, and soil texture on the Bashang Plateau.


Note: \* = Significantly correlated at *p* < 0.05, \*\* = Significantly correlated at *p* < 0.01.

A stepwise regression analysis was further used to determine the relationship between the maximum LAI for *P. simonii* and the dominant climatic, soil and topography variables (Table 6). For the whole Bashang Plateau, the MAP, sand and silt contents explained 79.7% of the spatial variation in the maximum LAI (*p* < 0.01). However, for zone I with a MAP < 350 mm, the MAP and elevation explained 98.5% of the spatial variation in the maximum LAI (*p* < 0.01). For zone III with a MAP > 420 mm, 77.5% of the spatial variation in the maximum LAI was explained by the MAP and MAT (*p* < 0.01). Additionally, 78.1% of the spatial variation in the maximum LAI was explained by the MAP (*p* < 0.01) in zone II with a MAP = 350~450 mm.

**Table 6.** Stepwise regression for the major variables influencing the spatial distribution of the mean maximum LAI for *P. simonii* forests on the Bashang Plateau.


Note: MAP = mean annual precipitation; ELV = elevation; MAT = mean annual temperature; Sand = sand content; Silt = silt content; \*\* = Significantly correlated at *p* < 0.01.

#### **4. Discussion**

#### *4.1. Model Validation*

The performance of the modified Biome-BGC model was evaluated by comparing it with the MODIS AET and NPP on both a local scale and a regional scale. Furthermore, the RWI was used to evaluate the simulated annual NPP based on the modified Biome-BGC model. It was demonstrated that the temporal and spatial dynamics of AET and NPP for *P. simonii* could be rationally simulated on the Bashang Plateau. The modified Biome-BGC model could be used to calculate the LAI for *P. simonii* considering the strong correlations among the AET, NPP, and LAI [38–41].

The MODIS data and simulated values differed for some study sites for *P. simonii*, which may have been due to limitations in the observed data. For zone I and zone II, some of the study sites were shelter forests planted at equal intervals with low coverage, which may have resulted in overestimations of the plant coverage for the study sites. For zone III, relatively few sites were sampled, and the soil data required for the modified Biome-BGC model were acquired from the Harmonized World Soil Database. The RWI of *P. simonii* was consistent with the NPP simulated by the modified Biome-BGC model in the Zhangbei area except for some dry years. This result may be caused by the number of tree ring samplings and the differences in tree growth of *P. simonii*. Another limitation of the simulated values was related to the uncertainty of the eco-physiological parameters used as model inputs. Some of these parameters were taken from Kang et al. [34] rather than from measurements made under the different climate and soil conditions at the study sites, which might have led to inaccurate estimations of AET and NPP [23,42]. Regardless, based on the three metrics (*MD*, *RMSE* and *R2*) for the observed and simulated AET and NPP, the modified Biome-BGC model has been proven to be an appropriate tool for simulating the relationships among climate, soil, and vegetation on the Bashang Plateau.

#### *4.2. Spatial Variations in AET and NPP and the Influencing Factors*

The spatial pattern of the mean AET on the Bashang Plateau was governed by a decrease in the MAP values from the southeast toward the northwest. The Pearson correlation coefficients between the AET and MAP ranged from 0.965 to 0.998, indicating that MAP was the major factor controlling the AET in the study area. Zhang et al. [11] and Jia et al. [12] found that the long-term mean AET and MAP were well-correlated for dominant species on the Loess Plateau. Liu et al. [43] reported that precipitation was the main driver of spatiotemporal changes in AET in semiarid and arid areas of China. These results are consistent with the findings of our study. Because groundwater levels in the semiarid plateau are generally 10–11 m below the land surface, precipitation is the main source of

soil water on the Bashang Plateau. A negative correlation was found between the mean AET and MAT, and the result is consistent with the study of Ma et al. [44] and Zhu et al. [45]. The reason for this result is that the AET in the study area is mainly controlled by the MAP and may also be affected by the actual sunshine duration and wind speed; therefore, the MAT has a limited influence on the AET. Additionally, the narrow range of the MAT may reduce the possibility of detecting significant variations in the AET in the study area. A positive correlation was found between the mean AET and clay content, and a negative correlation was found with the sand content. This finding is consistent with the study of Jia et al. [12]. Soil texture strongly influenced the AET [46]. This result was mainly due to the effect of soil texture on soil water storage. Khakural et al. [47] found that soil water storage was positively correlated with silt and clay contents but was negatively correlated with sand content. Fine soil texture has high water holding capacity, high available soil water and low drainage loss, and hence, a higher AET [12]. A negative correlation was found between the mean AET and elevation, indicating that elevation had a negative effect on the AET in the Bashang Plateau. This effect may be due to the fact that high elevation was correlated with a low MAT and hence a low AET [48,49].

In different subregions, the factors influencing AET distribution varied significantly. As the MAP decreased, the correlations between AET and precipitation gradually decreased. This result indicated that the influence of precipitation on the dynamic change of evapotranspiration decreased as precipitation increased. In drier areas, rainfall directly impacts water availability, thus affecting the AET. In zone I with a MAP < 350 mm, in addition to MAP, the AET was mainly affected by soil properties and elevation. This result was due to soil properties affecting soil water availability and hence affecting the AET in this area with limited precipitation [50,51]. There was no correlation between the AET and MAT because the region is drier than other regions and the water is the main limitation. For zone II with a MAP = 350~450 mm, the AET was mainly affected by the MAP, MAT and elevation but had no significant correlation with soil properties. This result may be because there was no obvious difference in soil properties in this region, and the simulation results failed to reflect the effect of soil texture on the AET. In zone III, with a MAP > 420 mm, in addition to MAP, the AET was negatively correlated with sand content. Sand content affected the AET by affecting the soil water holding capacity [52]. There was no correlation between the AET and MAT because the elevation in this region is relatively higher and no obvious differences exist. The influencing factors of the spatial variation in the NPP on the Bashang Plateau and in in zone I and zone II were similar to those of the AET. There was a strong correlation between the NPP and AET [53,54]. However, in zone III, a semi-humid region, the NPP was positively correlated with MAP and MAT but not significantly correlated with soil properties. The increase in temperature was conducive to *P. simonii* growth in zone III due to the low temperature. This finding is consistent with the studies of Cui et al. [55] and Chen et al. [56]. The temperature produced a greater effect on the NPP in humid regions than in arid regions.

#### *4.3. Optimal Plant Cover and the Influencing Factors*

The LAI for *P. simonii* was simulated by considering variations in the long-term climatic conditions throughout the period of 1980–2019. The optimal maximum LAI in the study area indicated that the plant-cover ranges were LAI = 1.3 for Kangbao, LAI = 1.9 for Shangyi–Zhangbei–Guyuan and LAI = 3.3 for Fengning–Weichang. The simulated values for *P. simonii* were within the LAI range of 0–6.2 with a mean of 2.14 given by Zhu [57] for vegetation on the Bashang Plateau. The optimal plant cover (expressed by the maximum LAI) was similar to that given by Jia et al. [12], with an LAI = 1.1–3.5 for deciduous broadleaf forests on China's Loess Plateau. This is because parts of the Loess Plateau and Bashang region of Hebei Province are located in China's agro-pastoral transitional zone and have similar climate and soil conditions.

The optimal-plant-cover variation showed significant regional differences in the study area. The maximum LAI for *P. simonii* indicated that the optimal plant cover with the lowest values was distributed in the west and higher values were distributed in the east of the Bashang Plateau. Precipitation was the key factor controlling the annual maximum LAI for *P. simonii* in the Bashang Plateau. Zhao et al. [58] reported that the precipitation gradient is the most important driver of LAI. As the MAP decreased, the correlations between the LAI and precipitation gradually decreased. Precipitation has a greater impact on soil water availability in drier regions, which indirectly affects plant growth [59]. In addition to precipitation, sand and silt contents were the main factors affecting the maximum LAI spatial distribution in the study area. This finding is consistent with the study of Jia et al. [12]. Coarse-textured soils have low water-holding capacity, high drainage loss and low soil water storage for plant growth. However, there is a lower moisture limit for plant growth on a medium-textured soil due to its water-retention characteristics [46,60].

The major factors influencing spatial distribution of optimal plant cover showed significant regional differences in the study area. In addition to the MAP, the mean maximum LAI of *P. simonii* in zone I was highly influenced by elevation, as it modulated the climate and/or water availability in this area with low precipitation and temperature levels. This result is consistent with the report of Jia et al. [12]. The maximum LAI was strongly affected by precipitation in zone II. The spatial variation in the mean maximum LAI was strongly correlated with temperature for *P. simonii* in zone III since the precipitation resources could satisfy *P. simonii* growth and the temperature was relatively low in this area. Similar to our findings, Kong et al. [61] found that the influence of temperature on vegetation variation in high latitude or mountainous regions was strong. According to the above results, the MAP, MAT, soil texture, and elevation are reliable indicators of the optimal plant cover of *P. simonii* in the Bashang Plateau.

#### *4.4. Implications for Afforestation and Forest Management*

Increasing vegetation cover is one of the most effective methods of countering land degradation. Introduced species such as *P. simonii* were widely planted mainly because they grow faster than native tree species and provide sand-fixation functions within a short period of time. Zhang et al. [19] found that at a *P. simonii* stand age of 40 years, the near-surface wind speed is reduced by more than 80%, while the surface coverage exceeds 70% and the wind erosion volume is nearly zero. Previous studies also found that planting *P. simonii* could improve soil conditions and the microclimate [62,63]. The *P. simonii*, however, is a species with high water consumption [64], and excessive planting of *P. simonii* increases the soil water deficit, consequently resulting in poor plant growth. As a result, an optimal plant cover is vital to maintain water balance and vegetation sustainability. In this study, the optimal plant cover for *P. simonii* was determined by the modified Biome-BGC model in the Bashang Plateau. Recent studies have shown that artificial mixed forests reasonably utilize forest resources and play a certain role in promoting nutrient accumulation compared with artificial pure forests. In this study, however, we considered the optimal plant cover of *P. simonii* pure forests because of the limitations of the model mechanism. Our research on the optimal plant cover of mixed forests will be strengthened in the future.

#### **5. Conclusions**

In this study, the AET, NPP and maximum LAI dynamics of *P. simonii* were simulated by the modified Biome-BGC model in the Bashang Plateau. The spatial distributions of optimal plant cover were determined by the maximum LAI. Model simulations confirmed the following scientific hypotheses of the present study: (1) The optimal plant cover for *P. simonii* was spatial heterogeneity on the Bashang Plateau. In general, the spatial distributions of the mean AET, NPP and maximum LAI increased from northwest to southeast. The optimal plant cover for *P. simonii* was 3.3 in the Fengning–Weichang area (MAP > 420 mm), 1.9 in the Shangyi–Zhangbei–Guyuan area (350 mm < MAP < 450 mm) and 1.3 in the Kangbao area (MAP < 350 mm). The MAP, elevation, soil texture and MAT were the main factors influencing the AET, NPP and maximum LAI. (2) The factors influencing the spatial distributions of the mean AET, NPP and maximum LAI were different for each subregion. For all subregions, the spatial variations in AET, NPP, and maximum LAI were primarily affected by precipitation. As the MAP decreased, the correlations between the AET, NPP, LAI and precipitation gradually decreased. In addition to precipitation, soil texture and elevation were the main influencing factors for the AET, NPP and LAI of *P. simonii* in zone I. Elevation was the main influencing factor for the AET, NPP and LAI of *P. simonii* in zone II. MAT was the main influencing factor for the NPP and LAI of *P. simonii*, while sand content was the main influencing factor for the AET of *P. simonii* in zone III. These quantitative results indicating optimal plant cover and influencing factors for different subregions on the Bashang Plateau could provide important guidance for non-native vegetation restoration.

**Author Contributions:** Conceptualization, Y.Z.; methodology, Y.Z. and W.L.; data curation, F.S., S.L. and B.X.; Supervision, J.Z.; Writing—original draft, Y.Z.; Writing—review & editing, Y.Z. and W.L. 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 (No.42101019, 42001027), Natural Science Foundation of Hebei Province (D2019205274, D2019403115), Science and Technology Project of Hebei Education Department (QN2019152), Key Laboratory of Agricultural Water Resources & Hebei Key Laboratory of Agricultural Water-Saving, Center for Agricultural Resources Research, Institute of Genetics and Developmental Biology (KFKT201903), Science and Technology Project of Hebei Education Department (BJK2022022, BJK2022031).

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

#### **References**

