Determining the Dominant Contributions between Direct and Indirect Impacts of Long-Term Urbanization on Plant Net Primary Productivity in Beijing

Abstract

Rapid urbanization exerts noteworthy impacts on the terrestrial ecosystem carbon budget, with pronounced effects in a metropolis such as Beijing, the capital city of China. These impacts include both Direct and Indirect Impacts. For instance, direct impacts influence regional Net Primary Productivity (NPP) by directly altering the vegetation coverage area. Concurrently, indirect impacts primarily affect regional NPP indirectly through climate change and urban vegetation management. How direct and indirect impacts contribute to the NPP is the core content of our research. Owing to that, we need to precisely assess the spatial and seasonal characteristics of the impact of urbanization in Beijing from 2000 to 2020. Firstly, a novel framework was proposed to analyze the impact of urbanization on NPP at the pixel level. Meanwhile, we employ the Proximity Expansion Index (PEI) to analyze urban expansion patterns. Results reveal that the direct impacts led to a cumulative NPP loss of 0.98 TgC, with the largest loss stemming from cropland conversion to construction land. During the last two decades, there has been a 56.87% increase in the area used for urban development in Beijing, a clear sign of swift urban expansion. Concurrently, this urban growth has had favorable indirect effects on NPP, with an average annual increase of 9.76 gC·m⁻²·year⁻¹, mainly observed in urbanized regions. Moreover, the seasonal analysis underscored that indirect impacts were primarily temperature-related, exhibiting higher values during autumn and winter within urban areas, indicating enhanced vegetation growth suitability in urban areas during these seasons. Our findings quantitatively examine the numerical relationship between direct and indirect impacts at a magnitude level. The carbon gain brought about by indirect impacts surpassed the carbon loss induced by direct impacts, with indirect impacts offsetting 29.41% of the carbon loss due to direct impacts. Ultimately, we advocate for enhanced greening initiatives in areas of Beijing with higher indirect impacts to achieve optimal carbon gain. This strategy might effectively reduce the negative impact of rapid urbanization on the carbon budget of terrestrial ecosystems.

1. Introduction

Urbanization is a dynamic and complex process involving changes in the economy, population, society, and space [1]. This process usually takes many years or even decades to complete, during which time the long-term effects, both direct and indirect, on local ecosystems can be both positive and negative [2]. Each city has its own unique pattern of urbanization [3]. In contrast to well-developed cities in Europe and North America, the rate of urbanization in the peri-urban regions of East Asian metropolitan cities is relatively high, indicating a distinct development rhythm [4,5]. For example, large-scale development of single-use land in urban construction has significant negative impacts on the ecological environment, including the squandering of natural resources and the loss of biodiversity [6]. Urbanization has produced irreversible changes to the productivity of terrestrial ecological vegetation [7,8], including the massive expropriation of natural land such as cropland and grassland and the massive increase in impervious surfaces. These changes have had a significant impact on climate and atmospheric conditions, with the urban heat island effect being one of the most typical examples [9,10]. The aforementioned changes brought about by urbanization have also affected the terrestrial carbon stocks [11]. Accordingly, it is necessary to analyze the impact of intra-metropolitan urbanization on the regional carbon balance of ecosystems. What are the dominant contributions of the direct and indirect impacts of long-term urbanization on a typical metropolitan city in East Asia? This is our key research question.

Net Primary Productivity (NPP) is the quantity of carbon fixed through photosynthesis by plant communities per unit of time and area. Not only is NPP an essential indicator of the dynamic exchange of net carbon between plants and terrestrial ecosystems and the atmosphere, but it also plays a vital role in the global carbon cycle [12,13]. Therefore, NPP was employed as the key ecological indicator in this study. Moreover, as a sensitive indicator of climate and environmental changes, NPP is an efficient tool for assessing the carbon budget of terrestrial ecosystems [14]. The analysis of the long-term impacts of urbanization on NPP necessitates remote sensing datasets with adequate spatial resolution and extended time series. However, common Moderate Resolution Imaging Spectroradiometer (MODIS) NPP products do not estimate NPP in urban building land areas due to limitations [15]. Therefore, we employed the Carnegie Ames Stanford Approach (CASA) model to estimate NPP in our study area to surmount the technical challenges of inadequate spatial resolution and the absence of NPP in urban areas. Finally, an NPP dataset with a 250-m resolution was obtained for the study area, which satisfies the temporal and spatial resolution requirements of this research [16,17].

Research on the impacts of urbanization on NPP has been conducted at different scales. For example, Zhao et al. [18] classified the impact as direct and indirect. Indirect impacts are defined as natural climate impacts and human factors, while direct impacts are defined as changes in land use types that contribute to a reduction in vegetation. At the city scale, Urumqi’s urban expansion directly caused a total NPP loss of 51.45 GgC, offsetting 20.86% of the carbon gain resulting from indirect impacts [19]. At the scale of urban agglomerations, urban expansion in the key development area in China’s central plains urban agglomeration has resulted in regional carbon losses of 4.69 TgC, with indirect impacts offsetting the carbon sink loss of 0.04 TgC [20]. During the period 2000–2010, at the national level, the decline in China’s NPP due to urbanization offset 5.42 percent of climate-driven growth [21]. It is important to note that urban expansion does not diminish the carbon sequestration capacity of vegetation, but it does reduce the growth rate of that capacity [22]. Globally, metropolitan cities have more significant positive indirect impacts than small and medium-sized developing cities [1]. However, few studies have investigated the spatial and seasonal heterogeneity of urbanization’s effects on NPP, and there is a need to analyze the differential impacts of urbanization on terrestrial vegetation NPP at the spatial and seasonal levels. Moreover, many previous studies have not focused on the long-term changes in carbon budgets during urbanization, only analyzing the differences in NPP at two selected time points [23,24,25]. Therefore, it is necessary to comprehensively explore the long-term impacts of urbanization on NPP. This research focuses on Beijing, a global metropolis. Metropolis has a high pace of urbanization and significant ecological impacts, making it an ideal region to investigate the mechanisms of urbanization’s impacts on ecosystems [26]. Beijing is still transitioning from accelerated development to maturity, and this trend will persist until 2030. Throughout this process, urbanization will have an ongoing impact on the terrestrial ecosystem carbon cycle. To sum up, the goals of this paper are as follows:

(1)

Develop an improved, innovative framework to analyze the direct and indirect impacts of urbanization on NPP, specifically at the pixel level;

(2)

Analyze the temporal and spatial changes of direct and indirect impacts and their driving mechanisms. Determine which impacts of long-term urbanization contribute the most to the impacts of urbanization on the NPP of Beijing;

(3)

Estimate the carbon loss attributable to Beijing’s urbanization over the past two decades and propose a green space planning strategy that will effectively reduce this loss.

2. Materials and Methods

2.1. Study Area and Data Sources

Beijing (28–41°05′N, 115°25–117°30′E) is situated on the northwest border of the North China Plain. From northwest to southeast, the altitude decreases, and the Taihang Mountains border the region to the west, the Yanshan Mountains and the Inner Mongolia Plateau to the north, the Huanghuaihai Plain to the south, and the Songliao Plain to the northeast, with a total area of 16,410 km². Beijing has a typical continental monsoon climate that is mild, temperate, and semi-humid. Spring is dry and breezy, whereas summer is humid and sweltering. 70% of annual precipitation is concentrated during the summer months. The weather in autumn is sunny, while winter is chilly and dry. Since the late 1980s, Beijing has undergone an urbanization process of international renown. Currently, Beijing is divided into four functional areas (Figure 1), namely the Capital Core (CC) area (Dongcheng District, Xicheng District); the Function Expansion (FE) area (Haidian District, Chaoyang District, Shijingshan District, Fengtai District); the New Developing (ND) area (Changping District, Shunyi District, Tongzhou District, Fangshan District, Daxing District), and the Ecological Conservation (EC) area (Yanqing District, Huairou District, Miyun District, Pinggu District, Mentougou District).

In this study, the datasets used include the MODIS NPP data set, the Normalized Difference Vegetation Index (NDVI) data, land use data, and meteorological data. The meteorological data, sourced from the China Surface Climate Data Daily Value Dataset (V3.0) provided by the National Meteorological Administration (https://data.cma.cn/, accessed on 15 June 2023.), includes monthly total precipitation, mean temperature, and solar radiation data. The land use data were obtained from the Land Use and Land Cover (LULC) dataset of Wuhan University, which covers a specific time series from 1990 to 2020 and is utilized to quantify the urban expansion process [27]. NPP and NDVI data are derived from the Google Earth Engine (GEE) platform (https://earthengine.google.com, accessed on 3 October 2023.), a geographical big data analysis, and the visual cloud platform. The NDVI data comes from the MODIS Collection6 MOD09A1 product, and monthly NDVI data were synthesized using the Maximum Value Composite (MVC) method on the GEE platform. The NPP data source is the MOD17A3HGF product, with a time series from 2001 to 2020 for subsequent accuracy verification. The Digital Elevation Model (DEM) data were downloaded from the Resource and Environmental Science Data Center (https://www.resdc.cn, accessed on 6 November 2023). The DEM was used as an input value for the calculation of the total solar radiation area.

During the data processing phase, the Ångström function was used to convert effective sunlight duration into solar radiation data to address the limited solar radiation stations in Beijing [28]. All meteorological data were processed using the ANUSPLIN specialized climatic interpolation software. Compared to spatial interpolation methods such as Kriging and Inverse Distance Weighting (IDW), ANUSPLIN yields the most accurate and applicable interpolation results with the lowest interpolation error [29]. To ensure data compatibility, the previously specified data were resampled to a resolution of 250 m using the WGS1984 coordinate system.

2.2. Methods

2.2.1. Estimating NPP Using the CASA Model

The CASA model is a light-use efficiency model that simulates the process by which vegetation absorbs photosynthetically active radiation and converts it into its own organic carbon. It estimates vegetation NPP using remote sensing data, temperature, precipitation, solar radiation, vegetation type, and land use type data [30,31]. In this model, NPP is estimated predominantly based on Absorbed Photosynthetically Active Radiation (APAR) and actual light use efficiency (ε). The equation is as follows:

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

(1)

where $\mathrm{t}$ represents time, $\mathrm{x}$ represents spatial position, $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},\mathrm{t})$ represents the NPP value of pixel $\mathrm{x}$ at time $\mathrm{t}$, and $\mathrm{A}\mathrm{P}\mathrm{A}\mathrm{R}$ and $\mathsf{\epsilon }$ value calculation equation utilizes (2)–(4).

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

(2)

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

(3)

$$\mathsf{\epsilon }(\mathrm{x},\mathrm{t})=\mathrm{T}\mathsf{\epsilon }1(\mathrm{x},\mathrm{t})\times \mathrm{T}\mathsf{\epsilon }2(\mathrm{x},\mathrm{t})\times \mathrm{W}\mathsf{\epsilon }(\mathrm{x},\mathrm{t})\times \mathsf{\epsilon }\ast$$

(4)

where $\mathrm{L}(\mathrm{x},\mathrm{t})$ represents the total solar radiation (MJ·m⁻²) at pixel $\mathrm{x}$ during month $\mathrm{t}$. $\mathrm{F}\mathrm{P}\mathrm{A}\mathrm{R}(\mathrm{x},\mathrm{t})$ denotes the fraction of incident Photosynthetically Active Radiation (PAR) absorbed by the vegetation canopy, reflecting the degree of $\mathrm{P}\mathrm{A}\mathrm{R}$ absorption by vegetation. The constant 0.5 is the proportion of solar radiation useable by vegetation (with a wavelength range of 0.4~0.7 μm) in relation to total solar radiation. $\mathrm{T}\mathsf{\epsilon }1$ and $\mathrm{T}\mathsf{\epsilon }2$ represent the impact of temperature on the light conversion rate, while $\mathrm{W}\mathsf{\epsilon }$ represents the water stress impact coefficient, which illustrates the effect of water conditions. $\mathsf{\epsilon }\ast$ is an abbreviation for the highest light conversion rate under ideal conditions. Estimating vegetation NPP requires the maximum and minimum values of NDVI and the Simple Ratio (SR) vegetation index for distinct plant types, in addition to the maximum light-use efficiency. This study cites NPP estimation parameters computed by Zhu Wenquan and others using measured data for diverse vegetation types in China [32].

2.2.2. Analysis of Urban Expansion Patterns

Urban expansion is a dynamic process of geographical spatial alteration that is typically characterized by the introduction of new urban land patches and the replacement of existing ones. The geographical link between these emerging urban patches and the existing urban areas is a critical question in studying the dynamic evolution of urban expansion. We use the Proximity Expansion Index (PEI) in this study to assess the degree of spread of new urban land patches. This offers a morphological viewpoint on the process of urban expansion [33]. The equation is as follows (5):

$$\mathrm{P}\mathrm{E}\mathrm{I}=\frac{1}{\mathrm{N}+(1-\frac{{\mathrm{p}}_{\mathrm{i}}}{{\mathrm{p}}_{\mathrm{n}}})}$$

(5)

where PEI is an abbreviation for Proximity Expansion Index, $\mathrm{N}$ is the number of multi-buffers for each new patch, indicating the approximate distance between old and new patches. ${\mathrm{p}}_{\mathrm{n}}$ is the area of the outermost buffer zone; ${\mathrm{p}}_{\mathrm{i}}$ is the intersection area between the buffer zone and the old patch; $\frac{{\mathrm{p}}_{\mathrm{i}}}{{\mathrm{p}}_{\mathrm{n}}}$ represents the boundary sharing rate between the old and new patches.

2.2.3. The Dynamic Influence Framework of Urbanization on NPP

The process of urban expansion has various impacts on regional ecosystems. Zhao et al. [18] proposed a conceptual framework to quantify these impacts, breaking them down into direct impacts (conversion of vegetation land into other land uses) and indirect impacts (natural climate variations or those caused by human factors). However, this framework cannot capture the continuous changes in urban influences over a specific time period. Based on Zhao’s research, Guan et al. [34] devised an analytical framework that discerns the ongoing effects of urbanization on regional NPP within specific durations. Based on previous research, this study improved the urban impact framework to express the intensity of urban expansion at a pixel level at a specific time. The revised framework for the impact of urbanization on NPP is as follows:

If only considering the impact of changes in urban land cover, this paper postulates that the NPP of a pixel entirely covered by vegetation before urbanization is set as $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{p}}$ (8):

$$\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{p}}(\mathrm{x},{\mathrm{t}}_0)=\frac{{\mathrm{M}}_{\mathrm{x}}}{{\mathrm{M}}_{\mathrm{x}}-{\mathrm{N}}_{(\mathrm{x},{\mathrm{t}}_0)}}\times \mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)$$

(6)

$$\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)=\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)+\frac{{\mathrm{N}}_{{\mathrm{t}}_0}-{\mathrm{N}}_{{\mathrm{t}}_1}}{{\mathrm{M}}_{\mathrm{x}}}\times \mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{p}}(\mathrm{x},{\mathrm{t}}_0)$$

(7)

where $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)$ represents the NPP estimated value for the $\mathrm{x}$ pixel during the ${\mathrm{t}}_0$ period. The $\mathrm{x}$ pixel is divided on average into $\mathrm{M}$pixel points with the total represented by ${\mathrm{M}}_{\mathrm{x}}$. ${\mathrm{N}}_{(\mathrm{x},{\mathrm{t}}_0)}$ signifies the number of non-vegetation pixels within the $\mathrm{x}$ pixel during the ${\mathrm{t}}_0$ period (${\mathrm{N}}_{(\mathrm{x},{\mathrm{t}}_0)}$ ≤ ${\mathrm{M}}_{\mathrm{x}}$), $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{p}}(\mathrm{x},{\mathrm{t}}_0)$ means the NPP value of the $\mathrm{x}$ pixel with full vegetation coverage before urbanization during the ${\mathrm{t}}_0$ period, and $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$ denotes the NPP value of the $\mathrm{x}$ pixel after urbanization in the ${\mathrm{t}}_1$ period (7). Based on this, according to $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$ during the ${\mathrm{t}}_1$ period, the direct and indirect impacts of urbanization on NPP can be derived. The equation is as follows (8) and (9):

$$\mathrm{D}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{I}\mathrm{m}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{t}(\mathrm{x},{\mathrm{t}}_1)=\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)-\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)$$

(8)

$$\mathrm{I}\mathrm{n}\mathrm{d}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t} \mathrm{I}\mathrm{m}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{t}(\mathrm{x},{\mathrm{t}}_1)=\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_1)-\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$$

(9)

where $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_1)$ represents the actual NPP value of pixel $\mathrm{x}$ after urbanization in the period ${\mathrm{t}}_1$, $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)$ is the estimated NPP value of pixel $\mathrm{x}$ in the period ${\mathrm{t}}_0$, and $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$ is the NPP value of pixel $\mathrm{x}$ in the period ${\mathrm{t}}_1$ post-urbanization, considering only the impacts of land cover change (direct impacts). Therefore, the difference between $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$ and $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_0)$ represents the direct impact of urbanization on NPP in pixel $\mathrm{x}$ during the period ${\mathrm{t}}_{0-1}$. The value of the indirect impacts is determined by subtracting $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{u}}(\mathrm{x},{\mathrm{t}}_1)$ from the actual NPP value $\mathrm{N}\mathrm{P}\mathrm{P}(\mathrm{x},{\mathrm{t}}_1)$ of pixel $\mathrm{x}$ during the urbanization period ${\mathrm{t}}_1$, representing the indirect impact of urbanization on NPP.

2.2.4. Analysis of Dynamic Changes in NPP Trend

To demonstrate the dynamic change trend in NPP in Beijing from 2000 to 2020, this study adopted the method of Sen’s slope in conjunction with the Mann-Kendall trend significance test to analyze every pixel in the research area. The computation formula for Sen’s slope is (10):

$$\mathsf{\beta }=\frac{\mathrm{n}\times {\sum }_{\mathrm{m}=1}^{\mathrm{n}}\mathrm{m}\times \mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{m}}-{\sum }_{\mathrm{m}=1}^{\mathrm{n}}\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{m}}{\sum }_{\mathrm{m}=1}^{\mathrm{n}}\mathrm{m}}{\mathrm{n}\times {\sum }_{\mathrm{m}=1}^{\mathrm{n}}{\mathrm{m}}^2-({\sum }_{\mathrm{m}=1}^{\mathrm{n}}\mathrm{m}{)}^2}$$

(10)

where $\mathsf{\beta }$ represents the trend value of NPP changes, where 1 < $\mathrm{m}$ < $\mathrm{n}$, $\mathrm{n}$ stands for the study duration ($\mathrm{n}$ = 20), $\mathrm{m}$ represents the year’s serial number, and $\mathrm{N}\mathrm{P}{\mathrm{P}}_{\mathrm{m}}$ represents the NPP value of year $\mathrm{m}$. When $\mathsf{\beta }$ > 0, it indicates an increasing trend of NPP, and the larger the value, the more NPP increases. Conversely, it signifies a decreasing trend of NPP if the value is negative.

To further analyze the significance of the simulated trends for each pixel point, the Z value of the Mann-Kendall test was used to assess the significance of the NPP change trend in Beijing [35,36]. When |Z| > 1.65, |Z| > 1.96, and |Z| > 2.58, the results pass the significance test with confidence levels of 0.9, 0.95, and 0.99, respectively.

2.2.5. Geoda Spatial Autocorrelation Analysis

To study the spatial distribution characteristics of the impact of urbanization on NPP, this article employs spatial autocorrelation analysis. Spatial autocorrelation measures whether spatial elements are clustered and includes both global and local spatial autocorrelation [37]. Global autocorrelation reflects the similarity in attribute values between a spatial element and its neighbors, represented by the Global Moran’s Index [38]. The calculation formula is as follows (11)–(14):

$$\mathrm{I}=\frac{{\sum }_{\mathrm{a}=1}^{\mathrm{n}}{\sum }_{\mathrm{b}\ne \mathrm{a}}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{a}\mathrm{b}}({\mathrm{x}}_{\mathrm{a}}-\bar{\mathrm{x}})({\mathrm{x}}_{\mathrm{b}}-\bar{\mathrm{x}})}{{\mathrm{S}}^2\left({\sum }_{\mathrm{a}=1}^{\mathrm{n}}{\sum }_{\mathrm{b}\ne \mathrm{a}}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{a}\mathrm{b}}\right)}$$

(11)

$${\mathrm{S}}^2=\frac{1}{\mathrm{n}}{\sum }_{\mathrm{a}=1}^{\mathrm{n}}({\mathrm{x}}_{\mathrm{a}}-\bar{\mathrm{x}}{)}^2$$

(12)

$${\mathrm{I}}_{\mathrm{i}}=\frac{1}{{{\mathrm{S}}_{\mathrm{i}}}^2}({\mathrm{x}}_{\mathrm{i}}-\bar{\mathrm{x}}){\sum }_{\mathrm{b}=1,\mathrm{b}\ne \mathrm{a}}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{a}\mathrm{b}}({\mathrm{x}}_{\mathrm{a}}-\bar{\mathrm{x}})$$

(13)

$${{\mathrm{S}}_{\mathrm{i}}}^2=\frac{1}{\mathrm{n}-1}{\sum }_{\mathrm{b}=1,\mathrm{b}\ne \mathrm{a}}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{a}\mathrm{b}}-{\bar{\mathrm{x}}}^2$$

(14)

where $\mathrm{I}$ represents the Global Moran’s Index; ${\mathrm{I}}_{\mathrm{i}}$ stands for the Local Moran’s Index; n is the total pixel count in Beijing; ${\mathrm{x}}_{\mathrm{a}}$ and ${\mathrm{x}}_{\mathrm{b}}$ are the values by which pixel a and b NPP are changed due to urbanization, respectively; $\bar{\mathrm{x}}$ is the average value of NPP changed by urbanization; ${\mathrm{w}}_{\mathrm{a}\mathrm{b}}$ is the sum of spatial weights for all pixels. If pixels a and b are adjacent, the corresponding element in ${\mathrm{w}}_{\mathrm{a}\mathrm{b}}$ is 1; otherwise, it is 0.

The Global Moran’s Index can determine whether there is spatial clustering among the elements within the study area. The value range for I is from −1 to 1. If I > 0, it indicates a positive spatial correlation among the spatial elements in that area; the larger the value, the more pronounced the spatial correlation [39]. When I < 0, it indicates a negative spatial correlation, and the smaller the value, the greater the spatial disparity. Local spatial autocorrelation analysis, on the other hand, can pinpoint the locations of spatial clustering among elements and is crucial for analyzing the distribution characteristics of specific spatial attributes in the study area.

3. Results

3.1. The Urban Expansion Process of Beijing over the Last 20 Years

Over the past two decades, Beijing has undergone rapid urbanization. Based on the PEI values (

Table 1. Statistics of patches under different urban expansion modes.

Expansion Type	PEI Value Range	2000–2005		2005–2010		2010–2015		2015–2020
		Number	Proportion	Number	Proportion	Number	Proportion	Number	Proportion
Infilling	(2/3,1]	1044	23.79	27	0.63	5	0.13	2	0.12
Edge-expansion	(1/2,2/3]	3157	71.93	3143	73.62	2701	70.43	1074	65.89
Leap	(0,1/2]	188	4.28	1099	25.74	1129	29.44	554	33.99

), this study categorizes Beijing’s urban expansion into three distinct patterns: infill expansion, edge expansion, and leapfrog expansion. Detailed analyses were conducted for four distinct time periods. The results indicate that edge expansion dominated in all four periods. Infill expansion was particularly evident between 2000 and 2005, where new urban patches frequently emerged between existing ones, suggesting a trend toward denser urban development. The leapfrog expansion pattern was most prominent from 2005 to 2015, with new patches extending in various directions, indicating an increasingly dispersed urban morphology. Between 2015 and 2020, the number of new patches resulting from all three expansion modes exhibited varying degrees of decline, indicating a gradual slowdown in Beijing’s urban expansion rate during this period.

During the urban expansion period, the primary changes in land use types are illustrated in Figure 2a. In 2000, 2010, and 2020, the urban land area was 2290.91 km², 3135.51 km², and 3593.77 km², respectively, reflecting an overall growth of 56.87%. During this period, the cropland area decreased by 1381.91 km², accounting for 24.76% of the total cropland area. Furthermore, the forest area increased during the urbanization process. To better display the afforested areas over the past twenty years, six primary land use categories (cropland, forest, grassland, water bodies, impervious, and barren land) were generated based on the GEE platform using Landsat images and human-computer interaction interpretation, with an accuracy exceeding 94.3% [27]. Comparisons of forest areas with Global Forest Change (GFC) data reveal that the spatial distribution of forests in this study is broadly similar to GFC data. These six primary land use categories were then reclassified into forest and non-forest types, and a forest/non-forest map of the study area was drawn. Figure 2b shows significant afforestation in the EC region, with an increase of 171.11 km², mainly concentrated in Huairou, Miyun, and other areas. During the research period, the newly afforested area and the areas of forest logging in Beijing are shown in Figure A1, Figure A2 and Figure A3 (Appendix A). The forest area increased from 47.48% in 2000 to 48.74% in 2020. This demonstrates that the afforestation plan initiated by the Beijing municipal government in 2000 has achieved noticeable success.

Based on the changes in different land use types, this study further analyzed the characteristics of land use type transformation during the urbanization process. The results show that from 2000 to 2020, approximately 12.73% of the land area underwent land use type changes, highlighting the urbanization process in Beijing. As illustrated in Figure 3a, the orange color represents the type of transformation from cropland to urban land, mainly concentrated in the FE and ND areas. Notably, the orange area is the largest, indicating that most of the newly added urban land is transformed from cropland. In the FE and ND areas, the type of transformation from cropland to urban land accounts for 86.65% and 79.45% of their total transformations, respectively, as seen in Figure 3b. This reveals that a large amount of highly productive cropland has been occupied during the urbanization process, leading to significant losses in vegetation productivity in many areas.

3.2. Spatiotemporal Changes in Net Primary Productivity

3.2.1. Temporal Characteristics of Net Primary Productivity

From 2000 to 2020, the overall annual average NPP of Beijing exhibited a wavy growth trend, with an annual growth rate of 5.16 gC·m⁻²·year⁻¹. In 2020, it reached its peak, where areas with an annual average NPP exceeding 800 gC·m⁻²·year⁻¹ accounted for 17% of the total study area. This is double the percentage from 2000, indicating that the vegetation’s carbon sequestration capacity within the study area is gradually strengthening year by year. Furthermore, as can be seen from Figure 4, the annual average NPP showed significant fluctuations during the period from 2000 to 2007. This is because, during this time, there was a higher prevalence of infill urban expansion, where new urban land patches often emerged between old ones, making urban development more compact. Compact cities tend to have a more significant negative impact on NPP. For example, decentralized cities, with their abundance of open spaces, facilitate the dispersion of internal pollutants [40]. However, compact cities, often associated with the heat island effect, affect the dispersion of pollutants, leading to higher concentrations of pollutants within these cities. This can easily exceed the ecosystem’s regulatory capacity, resulting in an unstable trend in the regional ecosystem state [41,42]. During this period, the disordered expansion of urban areas at the edges also severely impacted the NPP, which exhibited a fluctuating downward trend [43]. After 2007, the annual average NPP gradually stabilized, with areas having an annual average NPP of less than 400 gC·m⁻²·year⁻¹ accounting for 37%. Overall, urbanization significantly reduced the regional NPP. However, with the continuous promotion of ecological construction and the constant increase of green spaces within cities, the direct negative impacts brought about by urbanization have significantly been alleviated.

3.2.2. Spatial Characteristics of Net Primary Productivity

The spatial distribution of NPP in Beijing generally increases from east to west and from south to north. High NPP values are mainly concentrated in the EC area. To study the impact of urbanization on NPP, the spatial trend changes of NPP in Beijing from 2000 to 2020 were analyzed. The results show that 84.27% of the study area exhibits an increasing trend in NPP, of which 13.01% of the study area showed a highly significant increase in NPP. As depicted in Figure 5, the CC and FE areas have a larger area where NPP has increased significantly. For example, in the ecological management project of the Beijing-Tianjin Sand Source Control, measures such as converting farmland back to forest and afforestation of wastelands were taken. These ecological management projects have led to an increase in vegetation in urban areas, and human management has enhanced the vegetation’s carbon sequestration capacity. For instance, artificial irrigation during summer can mitigate the restrictive effects on vegetation growth due to water shortages.

Due to the influence of urbanization, the NPP levels in the four areas all showed growth trends at varying rates (Figure 6). Among them, the CC area had the fastest growth rate (6.65 gC·m⁻²·year⁻¹). This is because the CC area is already highly urbanized; the urban heat island effect generated internally promotes vegetation growth in winter, while the direct negative impacts of urbanization are relatively minor. In contrast, the ND area exhibited the slowest growth rate (4.34 gC·m⁻²·year⁻¹), indicating that urbanization has had a significant negative impact on the regional ecosystems of the ND area. Overall, urban expansion in Beijing has not reduced the carbon sequestration capability of vegetation, but it has slowed the growth rate of this capability.

3.3. Direct and Indirect Impacts of Urbanization on NPP

3.3.1. Inter-Annual Variation

The impact of urban expansion on NPP is complex and varies with the patterns of urban expansion. While urban land replacing vegetation cover directly damages urban NPP during urban expansion, natural climatic conditions, and human activities also positively influence the growth state of vegetation. Hence, this study categorizes the impact of urbanization on NPP into direct and indirect impacts. The mechanisms of these two impacts differ in Figure 7; the direct impact brought about by urbanization over the past two decades has consistently been negative and has been intensifying, with an annual direct impact of −2.87 gC·m⁻²·year⁻¹. This means that when vegetation cover is converted to urban land due to urbanization, the regional ecosystem incurs a carbon loss of 2.87 gC·m⁻²·year⁻¹. In contrast, excluding the indirect impacts caused by changes in land use types, despite fluctuations, the impacts have consistently been positive, with an annual average of 9.76 gC·m⁻²·year⁻¹, resulting in a carbon gain of 9.76 gC·m⁻²·year⁻¹. This indicates that the climatic changes and human activities triggered during the urbanization process have had a positive impact. In this context, the carbon loss due to urban expansion counterbalances 29.41% of the carbon gain resulting from the indirect influences.

3.3.2. Regional Differences

Due to the vast spatial extent of Beijing, there are often significant spatial variations within the metropolis. This study divided the analysis into four periods (the past twenty years have been divided into four stages on average) to calculate the direct and indirect impacts in different areas (Figure 8). The results indicate that the FE and ND areas experienced the most significant direct impacts due to changes in land use types, and these impacts were negative. Between 2000 and 2005, the direct impact peaked in the FE area (49.21 gC·m⁻²·year⁻¹), suggesting that the FE area underwent rapid urbanization during this period. From 2005 to 2020, the direct impact in the CC area was positive and continually increasing, indicating an increasing focus on greening efforts in the CC area and a gradual shift in land use type toward vegetation cover. In contrast, the indirect impacts were positive in most areas. Between 2010 and 2015, the indirect impact peaked in the ND area (109.87 gC·m⁻²·year⁻¹). During this period, the negative direct impact offset 23.92% of the NPP gain caused by the indirect impact. The annual averages for the ND and FE areas were both higher than those of the ecological conservation areas, indicating that the positive indirect impact on vegetation within urban areas exceeds that of non-urban forested regions.

3.4. Analysis of Indirect Impact Driving Factors of Urbanization

3.4.1. Correlation Analysis

The indirect impact of urbanization refers to its influence on plant growth, that is, the response of plant growth to changes in environmental conditions [18]. The direct impact of urbanization on regional NPP is theoretically a linear decline due to vegetation loss caused by land use changes. The indirect impacts of urbanization are more complex, often affecting vegetation NPP indirectly by altering the growth environment crucial for plant growth. Studies show that climate change indirectly affects regional vegetation NPP by influencing the Start of Photosynthesis Period (SOP) [44,45]. This study quantifies the indirect impacts using a modified urban impact framework, specifically the difference between actual vegetation NPP and theoretical vegetation NPP (considering only the impact of land use changes during urbanization) in a given area after urbanization [46]. Climatic factors are the dominant natural elements controlling vegetation growth, with temperature, precipitation, and solar radiation jointly determining the input of water and energy during the vegetation growth process [47]. Studies have shown that temperature and precipitation have significant effects on the growth state of plants, with these two factors making the most significant contribution to changes in the NPP [48,49,50]. Therefore, this paper selects temperature and precipitation for correlation analysis with the indirect impacts of urbanization, with the aim of determining which of these two factors is the main driving force behind the indirect impacts of urbanization. The results showed that, between 2000 and 2020, over 80% of the area exhibited a positive correlation between the indirect impact and temperature (Figure 9a). In contrast, only around 20% of the area showed such a correlation with precipitation (Figure 9b), suggesting that temperature is the primary driver of the indirect impacts. The regions showing a significant positive correlation between temperature and indirect impact were mainly concentrated in the EC and ND areas. In the ND area, vegetation growth was enhanced due to the rise in temperature brought about by rapid urbanization (including both natural climate change and the urban heat island impact) [51]. On the other hand, the EC area, being less impacted by urbanization, saw natural climate variations as the main driving force behind the indirect impacts. Notably, in some urban areas, temperature and the indirect impact showed a negative correlation, as vegetation growth can be inhibited when temperatures, exacerbated by the urban heat island impact, exceed the optimal range for plant growth.

3.4.2. Seasonal Heterogeneity

Direct impacts influence NPP by altering land use types and are, therefore, not affected by seasonal variations. Although fluctuations exist, these impacts are consistently negative. However, indirect impacts mainly affect the growth state of vegetation, where NDVI is an effective indicator reflecting the growth state of plants. As illustrated in Figure 10c, in summer and autumn (July to September), the climate is warm, and precipitation is plentiful, resulting in more significant indirect impacts and better vegetation growth. Conversely, in winter (October to January), the lower temperatures and dry conditions lead to lesser indirect impacts and poorer vegetation growth. The extent of indirect impacts is mainly influenced by seasons, necessitating further exploration of their seasonal heterogeneity. The results reveal that the peak indirect impacts for four areas predominantly occur in the summer (August), with the lowest values manifesting in the winter (February). This pattern aligns with the annual temperature distribution depicted in Figure 10b, reaffirming temperature fluctuations as the primary driver of indirect impacts. Conversely, in the CC area, the maximum indirect impact occurs in winter (January), showing a U-shaped curve with higher values in autumn and winter than in summer, indicating that urbanization’s indirect impact on vegetation growth is more pronounced during colder periods (Figure 10a). This finding is consistent with prior research conclusions, which posit that vegetation in urban areas of colder environments provides positive feedback to warming climates [52,53,54]. As observed from Figure 10b, the temperatures in urban areas (CC, FE, ND) during winter are consistently higher than in suburban areas (EC), with the FE area exhibiting the most pronounced difference, reaching its peak in autumn rather than summer. Thus, it becomes evident that during winter, when temperature is the primary limiting factor for plant growth, the elevated temperatures brought about by the urban heat island impact within the cities favor vegetation growth more than in the less urban-affected suburban areas, resulting in higher values in autumn and winter than in summer.

4. Discussion

4.1. Evaluation of NPP Simulation Results

Validation of NPP can be conducted using two methods: the empirical validation method and the relative validation method. Given the vastness of the study area and the challenges in obtaining empirical data, this study opted for the relative validation method. We compared the NPP data calculated using the CASA model directly with the MODIS NPP data downloaded from NASA’s official website. The MOD17A3 NPP data has been validated and endorsed by many studies, confirming its suitability for global and regional scale NPP research [16,55,56]. Consequently, many studies have juxtaposed model estimates with this product dataset to assess the reliability of their estimates [15]. As this study investigates the impacts of urbanization in five-year intervals and the MOD17A3 product starts from 2001, we chose the years 2001, 2005, 2010, 2015, and 2020 for accuracy verification. The data fit between the two sources has an R² value of 0.75 (p < 0.001), indicating that the NPP data computed in this study is precise and suitable for analytical research.

The average NPP estimates of both CASA and MODIS products were extracted from 2000 random points within the study area. This is because the MODIS product does not estimate the NPP of built-up lands, necessitating the exclusion of central urban areas during accuracy verification. As illustrated in the Figure 11, the NPP absolute values estimated by the CASA model are higher than those of the MODIS product. This discrepancy might be attributed to the MODIS product’s omission of urban area NPP estimates. Moreover, previous research has compared the average MOD17A3 NPP data with observational data, finding that the former is lower than the latter, but they are closely correlated [57]. Ultimately, this study also demonstrates that the NPP growth trend in urban centers is higher than in other areas, and it cannot be overlooked.

4.2. Necessity of Separating Direct and Indirect Impacts

Urban expansion has caused irreversible changes to terrestrial ecosystems, damaging the vegetation’s carbon sequestration capability and leading to carbon loss. However, the impacts of urbanization are complex, encompassing not only urban expansion but also anthropogenic influences and the resultant environmental impacts. Urban expansion has caused irreversible changes to terrestrial ecosystems, damaging the vegetation’s carbon sequestration capability and leading to carbon loss. However, the impacts of urbanization are complex. It does not only involve urban expansion, but anthropogenic impacts and their subsequent environmental consequences are equally pivotal. Through a per-pixel analysis of NPP trends from 2000–2020 in the study area, it was found that both the CC area and FE area exhibited a significantly increasing trend in NPP. In a comparison between the CC and EC areas, there’s a growth trend discrepancy of 1.08 gC·m⁻²·year⁻¹. This indicates that the positive indirect impacts of urbanization ensure an additional 1.08 gC is produced by vegetation per square meter annually. The above conclusion is consistent with findings that vegetation’s carbon sequestration capability is enhanced in urban environments under human intervention [1,58]. For instance, under the context of urban vegetation management, the rise in urban environmental temperature can advance the start of the growing season and, in colder conditions, delay its conclusion, thus resulting in longer vegetation growing season in urban areas compared to suburban regions [59,60,61]. The findings of this research, derived from considering indirect impacts, underscore their importance and indicate they should not be overlooked.

By calculating the linear trend of indirect impacts across different areas, as shown in

Table 2. Annual average growth rate of NPP and NPP_II in different areas.

Unit: gC·m−2·year−1	Period	ND	FE	CC	EC	BJ
NPP	2000–2020	4.34	6.23	6.65	5.57	5.16
NPPII	2000–2020	0.41	0.29	0.48	0.09	0.22

Note: ND represents the New Development area; FE stands for the Function Expansion area; CC represents the Capital Core area; EC represents the Ecological Conservation Zone; BJ represents the entire Beijing city; NPP_II is the value of NPP that changes due to the indirect impact of urbanization.

below, the EC area experienced the smallest increase. This suggests that vegetation in urban areas undergoes a more significant indirect influence than in natural settings. This finding aligns with previous research [58], emphasizing that the positive indirect impact of urbanization on vegetation growth is most evident in urban regions where vegetation is relatively scarce. This further underscores the many enhancing impacts urbanization brings to vegetation growth in urban areas. In summary, if only the direct impact of land use type changes is considered, overlooking the indirect positive impacts on vegetation growth, this study would only conclude that urbanization in Beijing over the past 20 years has reduced the city’s NPP. Therefore, it is necessary to dissect the impacts of urbanization on NPP.

4.3. Spatial Heterogeneity of Direct and Indirect Impacts

The impact of urbanization on NPP might vary due to differences in temporal and spatial scales; previous studies have predominantly focused on temporal scales, often neglecting the spatial impact variations caused by different areas and vegetation types. The direct impacts arise as urban land use replaces plant cover, affecting regional NPP; however, different types of vegetation result in varying carbon storage losses. For instance, Jiang et al. [62] found that urban expansion significantly reduces the NPP of areas with high NPP values, such as cropland and forests, in suburban regions.

This study utilized Geoda software to calculate the Global Moran’s I for direct and indirect impacts over four periods separately, evaluating the spatial characteristics of the direct and indirect impacts. The results show that the Global Moran’s I for both types of impacts is greater than 0 (p < 0.001), illustrating a significant spatial autocorrelation for both impacts in Beijing. For all four periods, the Global Moran’s I for indirect impacts was consistently higher than that for direct impacts (Figure 12), suggesting that indirect impacts have a stronger spatial correlation.

However, Global Moran’s I only indicates whether the two impacts are aggregated but does not reveal the location and intensity of the aggregation. Hence, the local indicators of spatial association (LISA: Local Moran’s I analysis) and significance for the two impacts were computed. The figure below further illustrates the geographic spatial distribution characteristics of the two impacts. Beijing’s urbanization process can be inferred from the LL (Low-Low) clusters of direct impacts. The LL (Low-Low) clusters, having negative values, symbolize the negative impacts brought by urbanization, mainly reflecting the changes in NPP caused by the transformation of LULC. Over the past 20 years, the expansion of Beijing’s outskirts has been dominant, mainly occurring in the ND area, where the NPP is most directly impacted (Figure 13i–l). This aligns with previous studies where the expansion primarily occurs around larger original cities [63]. It is noteworthy that the land occupied by the peripheral expansion is typically fertile and highly productive, exacerbating significant loss in suburban vegetation coverage [64,65]. Between 2015 and 2020, the number of HH (High-High) clusters of direct impacts within urban areas increased, indicating that by then, urban areas had started emphasizing the development of green spaces.

The spatial correlation of indirect impacts is stronger. From 2005 to 2010, the HH (High-High) clusters of indirect impacts were highly aggregated in urban areas. During this period, the growth of urban vegetation improved significantly due to indirect impacts (Figure 13n). From 2010 to 2020, the HH (High-High) clusters of indirect impacts were more dispersed. This is because the urban expansion pattern gradually shifted to a leapfrog mode, with new urban patches sparsely distributed around the city, far from existing patches. The leapfrog expansion may be due to human activities and living needs, such as the development of new farmland and irrigation fertilization, causing the HH (High-High) clusters of indirect impacts to spread rapidly outward (Figure 13o–p). In summary, the HH (High-High) cluster areas of indirect impacts can effectively improve the growth status of vegetation, so establishing new urban green spaces in these areas can achieve higher carbon gains.

4.4. Exploring the Impact of Urbanization on Regional Carbon Budgets

Terrestrial ecosystems, as an essential component of the global carbon cycle, possess the capability to fix and absorb carbon dioxide, playing a crucial role in offsetting a portion of industrial carbon emissions [66,67,68]. Urbanization has significantly impacted ecosystems in recent years; thus, it is necessary to analyze its impact on regional carbon budgets further. NPP_D represents the total NPP directly affected by regional urbanization, indicating carbon loss post-urbanization. Research suggests that over the past two decades, vegetation covering approximately 1321.26 km² in Beijing has been replaced by other land cover types, resulting in a carbon loss of 0.985 TgC (

Table 3. Changes in carbon budget and NPP_LV after urbanization in different areas.

	Period	ND	FE	CC	EC	BJ
NPPD (TgC)	2000–2005	−0.196	−0.077	−8.929 × 10−5	−0.045	−0.319
	2005–2010	−0.191	−0.050	−1.727 × 10−5	−0.043	−0.284
	2010–2015	−0.191	−0.035	1.439 × 10−5	−0.049	−0.275
	2015–2020	−0.073	−0.016	−1.067 × 10−4	−0.018	−0.107
	2000–2020	−0.651	−0.179	−1.988 × 10−4	−0.156	−0.985
NPPLV (gC·m−2)	2000	657.03	800.87	262.37	612.69	642.51
	2020	1160.21	1380.96	444.91	827.23	1002.17

Note: NPP_D represents the total NPP lost in the area due to urbanization; NPP_LV represents the actual NPP when the area is fully covered with vegetation.

). Among the different areas, the ND area contributes the most, with losses reaching 66.09% (0.651 TgC), primarily due to the substantial loss of cultivated land in this area. Except for the CC area, the EC area shows the least variation, highlighting the critical position of the EC area in the sustainable development of the capital. Notably, from 2015 to 2020, carbon losses due to urbanization in Beijing were reduced by half, with the ND area showing the most significant decrease. This underscores the success of Beijing’s efforts in alleviating non-capital functions, mitigating “ metropolis syndrome” [69].

Post-urbanization, the entire area’s NPP_LV rose to 359.66 gC·m⁻², marking a 1.56-fold increase from 2000. This indicates that an additional 847.08 km² of green space is needed to offset the carbon loss during urbanization. Furthermore, the higher the NPP_LV after urbanization, the more apparent the positive effects of the indirect impacts of urbanization. Therefore, the vegetation requirements to offset carbon losses vary. For instance, the ND, FE, CC, and EC areas, respectively, require 748.23 km², 766.25 km², 779.17 km², and 978.59 km² of green space to compensate for carbon loss during urbanization. The ND area requires the least amount of green space. It is the ideal area for converting land into vegetation, followed by the FE area. The CC area, the central urban area, does not have much area converted into vegetation. Based on the analysis of spatial heterogeneity of indirect impacts, the ND area, particularly Tongzhou and Daxing districts, has a higher concentration of HH (high-high) clusters, indicating significant indirect impacts. In summary, urban greening should proceed according to the ND, FE, and EC areas. When constructing green spaces, the Tongzhou and Daxing districts in the ND area might be prioritized for their potential superior carbon sequestration benefits. Rational urban greening is crucial for achieving sustainable urban development. It can effectively counteract the carbon sink loss brought about by urban sprawl.

4.5. Limitations and Future Outlook

This study first uses the CASA model to estimate the NPP of terrestrial vegetation. Then, an improved urbanization impact framework is employed to analyze the impact of urbanization on NPP. This impact is further subdivided into two distinct impacts for in-depth analysis. This research contributes to a comprehensive understanding of the impact of urbanization on carbon storage. Offering insights into how to compensate for carbon storage losses during urban expansion. Urban carbon storage losses can be mitigated through more detailed green space planning. However, there are certain limitations to the study. Firstly, indirect impacts are complex, and human-induced effects can be overlooked when environmental factors dominate. Secondly, the urban heat island effect generated within cities can accelerate soil carbon decomposition while promoting vegetation growth in urban areas during winter, potentially offsetting the carbon gains from indirect impacts. Therefore, to better understand how to mitigate the impact of urbanization on carbon storage, it is necessary to incorporate the changes in soil carbon storage into further research. Lastly, due to the accuracy of the CASA model in estimating NPP being influenced by input data and model parameters, future studies should focus on higher-resolution data and optimized models.

5. Conclusions

This study dynamically analyzes the direct and indirect impacts of urbanization on regional NPP in Beijing from 2000 to 2020. It was determined that indirect impacts are the significant contributors to the impact of urbanization on NPP. The regional NPP was estimated using the CASA model, and the urban expansion process was further analyzed using PEI. We proposed a framework predicated on pixel-level analysis of urban expansion extent, aimed at dissecting the dynamic composite impacts on NPP amid urban expansion processes. Over the past 20 years, Beijing has undergone rapid urbanization, with a metropolitan area growth of 56.87%. This led to a significant loss of cropland, accounting for 24.76% of the total cropland area. During this period, NPP in Beijing showed a fluctuating upward trend with an average annual growth rate of 5.16 gC·m⁻²·year⁻¹. Spatially, 32.34% of the area showed a significant increase in NPP trends. High NPP values were mainly concentrated in the northern and northwestern regions, with lower values in the southern and southeastern areas. The direct impact harms urban NPP by directly replacing vegetation cover, and this increases with the progression of urbanization. The mechanisms causing NPP changes due to indirect impacts are more complex. The study found that indirect impacts are significantly related to temperature, with evident seasonal differences within urban areas; the results are higher in autumn and winter than in summer. After urbanization, the total carbon loss in the region reached 0.985 TgC, with the highest contribution from the ND area. To further compensate for carbon sink losses, increasing the green space in development zones with significant indirect impacts will enhance the offsetting effect against direct impacts. The research results are of considerable importance in maintaining the carbon balance of urban ecosystems against the backdrop of accelerated urban development and climate change in the future.

Future urban development should avoid becoming too compact, adopting a more dispersed pattern to expand into surrounding areas and reserve space for urban greenery construction. The construction of urban green spaces should be prioritized in areas where the indirect effects of urbanization have a higher positive impact and between new and old urban land patches by creating more green spaces within the city to mitigate the pollution caused by rapid urbanization.