The Guanzhong Plain Urban Agglomeration (104°34′–112°34′E, 33°34′–36°56′N), centered around Xi’an, is located in the inland area of northwest China and is the second largest urban agglomeration in western China. This urban agglomeration covers an area of 107,000 km² and belongs to the warm temperate continental monsoon climate zone, bordered by the Qinling Mountains in the south and the Yellow River in the east, situated between the Qinba Mountains and the Loess Plateau. In 2020, the GPUA had a resident population of 38.87 million and a regional GDP of CNY 2.2 trillion.

The Guanzhong Plain Urban Agglomeration Development Plan officially delineated the planning area of the GPUA. This area includes Xi’an, Baoji, Xianyang, Tongchuan, Weinan, and Shangluo City (Shangzhou District, Luonan County, Danfeng County, and Zhashui County) in Shaanxi Province; Yuncheng City (except Pinglu County and Quanqu County) and Linfen City (Yaodu District, Houma City, Xianfen County, Huozhou City, Quwo County, Yicheng County, Hongdong County, and Fushan County) in Shanxi Province; as well as Tianshui City and Pingliang City (Kongdong District, Huating County, Jingchuan County, Chongxin County, and Lingtai County) and downtown Qingyang (Xifeng District) in Gansu Province. Because of data availability, this study uses the municipal area as the boundary, including the city-wide areas of Shangluo, Yuncheng, Linfen, Pingliang, and Qingyang. The study area is shown in Figure 1.

1) Vector boundaries of the study area were obtained from the 1:1 million National Basic Geographic Database of the National Geographic Information Resource Catalog Service System (https:// www.webmap.cn/); (2)Land use data at a spatial resolution of 30 m were obtained from the Resource Environment and Science Data Center of the Chinese Academy of Sciences. These data are based on the Landsat series of satellite images and were generated through human-computer interaction and manual visual interpretation (https://www.resdc.cn/). According to the land use classification standard of the Chinese Academy of Sciences, the study area is divided into six primary types, namely, farmland, woodland, grassland, water body, construction land, and unused land. Social and economic data as well as energy consumption data ranging from 2001 to 2021 were obtained from the “Gansu Province Statistical Yearbook” (2001–2021), “Shaanxi Province Statistical Yearbook” (2001–2021), “Shanxi Province Statistical Yearbook” (2001–2021), “China City Statistical Yearbook” (2001–2021), and statistical bulletins of cities and towns from 2000 to 2020. Food price data were obtained from the “National Compilation of Information on Costs and Benefits of Agricultural Products” (2001–2021).

The land use transfer matrix reflects the quantitative transfer status of each land use type in the form of a two-dimensional matrix over a 20-year period. This matrix is used to visualize the quantity and direction of changes between time periods and land use types. S i j = S 11 S 12 ⋯ S 1 n S 21 ⋯ S n 1 S 22 ⋯ S n 2 ⋯ ⋯ ⋯ S 2 n ⋯ S n n (1) where, S i j is the land area that is transferred from land use i to land use j ; n is the number of land classes; i is the land use type at the beginning of the study period; j is the land use type at the end of the study period.

The dynamism model is expressed by calculating the average annual rate of change in the quantity of a particular land use type in the region at the end of the study period and at the beginning of the study period. The degree of single land use dynamic is mainly used to reflect changes in specific land use types in the study area, indicating the rate of change and the magnitude of change in different land use types over time. The calculation formula is: K = U b − U a U a × 1 T (2) where, K > 0 means that the land scale of the type is expanding, K < 0 means that the land scale is shrinking, U a 、 U b are the number of single land use/cover types at the beginning and the end of the study period, respectively; T is the length of the study time period, and when the time period is set to 1 year, the attitude of the single land use dynamics is the annual rate of change of the utilization/cover type of the land.

Comprehensive land use dynamics responds to the rate of change in the overall type of land use in the study area over a given period of time, responding to comprehensive regional land use changes. The calculation formula is: L C = ∑ i = 1 n ∆ L U i − j 2 ∑ i = 1 n L U i × 1 T × 100 % (3) Where, L C is the comprehensive land use dynamics, L U i is the area of type i land use type at the beginning of the study, ∆ L U i − j is the absolute value of the area of type i land use type converted to non-type i land use type during the study period; T is the study period, and n is the number of land use types.

The carbon emissions in the study area were calculated by multiplying the area and carbon emission factors of each land type and then summing the results. The main land types that can be directly measured for carbon emissions are types of carbon sink land, including woodland, grassland, water body, and unused land; these are measured based on the carbon emission factors obtained by previous studies. The direct carbon emission calculation formula is: E l = ∑ e i = ∑ S i × δ i (4) where, E l denotes direct carbon emissions; e i denotes carbon emissions from different site types; S i and δ i denote the area of the site type and the carbon emission factors, respectively. In this study, the carbon emission coefficients of each land type were determined based on the research results of Cai et al. (2005), Fang et al. (2007), Duan et al. (2008), and Lai (2010). Table 1

Carbon emissions from construction land originate from fossil energy consumed in the production process. Because energy consumption varies greatly among different regions, indirect measurement is used to calculate it. The total carbon emissions from construction land are estimated indirectly by accounting for energy consumption of oil, coal, natural gas, and electricity, combined with data obtained from the China Energy Statistics Yearbook and the 2006 IPCC Guidelines for National Greenhouse Gas Emissions Inventories. The formula for calculating carbon emissions from construction land is as follows: E φ = ∑ e i = ∑ i = 1 n m i × β i × θ i (5) where, E φ represents the total carbon emissions from construction land; e i represents the carbon emissions generated by energy source i; m i is the consumption of various fossil energy fuels; β i is the standard coal conversion factor of energy sources; θ i is the carbon emission factor of each energy source. Table 2

The carbon emission intensity is determined by the area of each land use type and the carbon emission factor, where a larger carbon emission intensity means a larger carbon emission from each grid. The calculation formula is: C R I = ∑ i n S i P i S (6) where, C R I is the carbon emission intensity; S i and P i are the area of gridded land type and carbon emission coefficient, respectively; S is the total area of the gridded network.

In reference to the results of Xie et al. (2015), in this study, the ESV equivalent factor of one standard unit of ESV was defined as one-seventh of the economic value of natural food production per unit area of farmland per year nationwide. ESV correction and accounting were performed. Xie et al. (2015) corrected the ESV equivalent factor table adapted to the actual situation in China, which is more suitable for the actual application in China, as shown in Table 3. Based on Table 3, this study utilizes the following methods to make corrections to match the actual social development and regional characteristics of the GPUA.

The calculation formula is as follows: D = 1 7 × ∑ i = 1 n P i G i S i (7) where, D is the ESV for a standard equivalent factor; P i is the unit price of crop category i; G i is the yield of crop category i; S i is the planted area of crop category i.

The ESV coefficient is calculated by the formula: V C j = ∑ E j i × D (8) where, V C j denotes the coefficient of ESV of the jth land; E j i is the equivalent factor of the ESV of the ith land use ecological service function of the jth land; D represents the value of a standard equivalent factor.

The spatial clustering and dispersion between carbon emission intensity and ESV intensity was examined using the bivariate Moran’s index. The calculation formula is as follows: I = ∑ i = 1 n ∑ j = 1 n W i j x i − x x j − x ¯ S 2 ∑ i = 1 n ∑ j = 1 n W i j (9) I i = x i − x ¯ S 2 ∑ j = 1 n W i j x j − x ¯ (10) where, I and I i denote the global Moran’s I and local Moran’s I of carbon emissions and ESV intensity, respectively. The Moran’s I index should be in the range of [−1,1]. When Moran’s I index is greater than 0, it indicates the existence of spatial positive correlation, when Moran’s I index is less than 0, it indicates the existence of spatial negative correlation of agglomeration effect, and when Moran’s I is equal to 0, it indicates the non-existence of spatial relationship. Where, n represents the number of spatial grid units; x i and x j are the attribute values of cell grids i and j , respectively; S 2 is the variance. The bivariate LISA clustering diagram can illustrate the degree of association between the two at a certain significance level. This is mainly used to examine ESV and carbon emissions in the study area in terms of their spatial location, as well as to find the correlation and difference between both variables.

Figure 2 shows the land use types of the GPUA from 2000 to 2020. According to the land use classification standard, the GPUA is dominated by farmland, grassland, and woodland, which together account for about 94.7% of the total land area in 2020. In terms of the 2020 land use structure, for example, farmland accounts for about 39.9% of the total land area of the region, grassland accounts for about 33% of the total land area of the region, and woodland accounts for about 21.8% of the total land area of the region. The GPUA has a beneficial geographical location. The northern part of the urban agglomeration is a loess hilly and ravine area with relatively high vegetation coverage and rich grassland and forest resources. The central part consists of the Guanzhong Plain and the Fenwei Plain and is adjacent to the Weihe River water system, which is suitable for agricultural development. The southern part depends on the Qinling Mountains, which have good ecological diversity. During the period of 2000–2020, farmland decreased by 3,872.546 km² in total, which is the main category of land transferal. Construction land is the land type that increased the most, with an increase of 2002.442 km², whereas grassland increased by 892.98 km², woodland increased by 890.77 km², and water body and unused land increased by 67.351 km² and 19.003 km², respectively. To more intuitively understand the structural characteristics of land use change in the GPUA during this period, a land use transfer matrix was used to calculate land use dynamics. The land use transfer matrix was constructed and the land use dynamic attitude was calculated. The specific results are shown in Tables 4, 5. Overall, unused land has the highest land use motivation, followed by water body and construction land. Unused land and water body are smaller because of their land use area, while construction land increased more because of areas that were transferred in. Although the overall situation from 2000 to 2020 showed an increase in the areas of all land types (except for the largest area of transferred out farmland), the transfer-in attitude of woodland showed a slight decrease after 2005. The change in the transfer attitude of grassland was clearly slower than that of construction land. From 2000 to 2020, the largest amount of land that was transferred out was farmland, and the main transfer directions were grassland, construction land, and woodland. With the continuous promotion of the project of returning farmland to forest and grass and the project of soil erosion control in the northern loess hills and ravines, the degree of vegetation cover increased year by year. Moreover, wasteland and bare land were gradually transformed into grasslands and forests. According to the China Forestry and Grassland Statistical Yearbook, during the period of 2000–2020, the forest coverage rate of Shanxi Province increased from 13.3% to 20.5%; the forest coverage rate of Shaanxi Province increased from 32.6% to 43.1%; the forest coverage rate of Gansu Province increased from 6.7% to 11.3%. The prevention of soil erosion and the return of cultivated land to forest and grassland made great achievements. During this period, the GPUA expanded rapidly: construction land increased, population increased, urbanization expanded rapidly, regional arable land was annexed, and the intensity of land development and utilization also increased. The expansion met the demand for rapid economic expansion and development of this region.

In this study, the carbon emissions of urban land use in the GPUA are measured via carbon emissions, carbon sources/sinks, and net carbon emissions of each category. The dynamic change characteristics of each indicator in 2000–2020 are analyzed, and the results are shown in Table 6; Figure 3. As shown in Table 6, the overall net carbon emissions from land use in the GPUA followed a significant growth trend from 40,815.272 × 10⁴ t in 2000 to 352,655.178 × 10⁴ t in 2020. The growth trend of carbon emissions from construction land is basically consistent with the growth trend of net carbon emissions from land use, and carbon emissions from construction land have become the main source of the growth of net carbon emissions. The farmland area gradually decreased over the study period, and accordingly, the carbon emissions from cultivated land followed a decreasing trend from 289.77 × 10⁴ t in 2000 to 273.407 × 10⁴ t in 2020. Carbon sinks followed a first decreasing and then increasing trend during the study period. Woodland was the main source of carbon sequestration, and the carbon sequestration capacities of grassland and water body were relatively low. The carbon sequestration capacity of woodland showed a small increase over 2000–2010 and was more stable from 2010 to 2020. The carbon sequestration effect of grassland and water body remained basically unchanged during the study period. Unused land showed a decreasing trend from 2000 to 2015 and a significant increase from 2015 to 2020, which matches the trend of carbon sink changes.

Analyzing the data in terms of carbon source/sink composition showed that carbon emissions from construction land were the most important carbon source, accounting for 99.30%–99.92% of the total carbon source and increasing significantly year by year. As one of these carbon sources, the contribution rate of carbon emissions from farmland was significantly smaller than that of construction land, but the area share was high. Carbon emissions per hectare of construction land represent the carbon intensity and energy consumption structure of construction land as a land use category in each study year. The average carbon emission intensity of construction land in the study area in 2000, 2005, 2010, 2015, and 2020 were 843.725 t/hm², 1586.923 t/hm², 2898.595 t/hm², 4211.434 t/hm², and 5153.703 t/hm². The carbon emissions of construction land in 2020 were about six times of those in 2000. This result highlights the great potential for carbon emission reduction on construction land. Regarding the carbon emissions composition of carbon sinks, woodland and grassland are the main carbon sinks, together accounting for about 98.5% of the total carbon sinks. However, woodland and grassland, which cover 54.8% of the land area of the study area, have an extremely important carbon sequestration capacity, while water body and unused land have lower carbon sequestration capacities. Because of the great gap between the amount of carbon sources and carbon sinks, under the integrated regional eco-economic development strategy, the GPUA has begun to focus on controlling carbon sources and increasing carbon sinks to maintain ecosystem services and functions.

Because of differences in the class structure and development level in this urban agglomeration, to analyze the change of carbon emissions more intuitively and enhance comparability, the carbon emission intensity index was analyzed in depth. This index was divided into five intervals of −1700 to −900, −900 to 0, 0 to 300, 300 to 600, and 600 to 1100, and the results are shown in Figure 3. The maximum value of carbon emission intensity in the GPUA decreased from 1053.24 (2000) to 1048.66 (2020). This indicates that ecological and environmental quality issues have attracted attention during the development of urban agglomerations and the uncontrolled brutal development pattern of relying on energy and expanding land has been restrained. High values of carbon emission intensity were mostly distributed in the central, western, and northwestern parts of the urban agglomeration. In addition, the trend of spatial differentiation shows a distribution pattern of high intensity in the central region and low intensity in the surrounding regions, which is roughly similar to the distribution of land use structure of farmland and grassland. High altitude and slope areas are rich in ecological resources, and these areas also feature low land development, well-protected ecosystems, and low carbon emission intensity. High intensity carbon emission areas are mainly located on the banks of the Weihe and Yellow River as well as in the Guanzhong Plain. These areas are mainly influenced by the expansion of urban and rural settlements and changes in the scale of agricultural production land. The development of agriculture and the expansion of the scope of urban settlements have led to a high degree of land development. The high intensity of human activities has directly led to the problem of high carbon emissions from land use, which is difficult to circumvent.

To apply coefficient correction to the actual production capacity of the Guanzhong Plain Urban Agglomeration, the grain yield correction method was used. According to the yield of major grain crops, sown area, and grain purchase price in the study area from 2000 to 2020, the coefficient of ESV of the GPUA was obtained by combining the land use type and the corresponding ecosystem. The results are shown in Table 7.

The total ESV of the GPUA from 2000 to 2020 showed fluctuating changes and an overall increasing trend. The total value increased from 215,263.7 million yuan to 216,776.2 million yuan, which represents an increase of 0.700% (Table 8). Judging from the ESV of each land type, woodland, grassland, and farmland are the main contributors to the total ESV, accounting for 44.674%, 26.160%, and 21.144%, respectively. Among them, ESV contributed by farmland decreased year by year, from 21.836% to 20.459%. The remaining land types showed an overall increasing trend. The functional structure of land use indicates that during the period of 2000–2020, farmland decreased by a total of 3872.546 km², while construction land was the main category of increase, increasing by 2002.442 km². At 67.351 km², the increase in the water area was less. However, the increase in ESV added during the study period reached 4.403%, indicating that water bodies have high ESV intensity and are also the main contributors to the total ESV. Water bodies are also the main contributing land type and hold a very important ecological status. The interconversion among land classes is also an important reason for the fluctuation of total ESV during the study period.

To visually analyze the spatial distribution characteristics of ESV in the study area, while avoiding the influence of the edge area of the study area, a 5 km × 5 km grid was created for an equally spaced sampling of the study area based on 30 m resolution land use data. The study area was divided into 6934 units, and the ESV intensity results of each study unit were identified. The results are shown in Figure 4. In this study, the calculated ESV intensity was divided into four zones, I, II, III, and IV, and the range of values defined for the four zones corresponded to: 0–20, 20–35, 35–55, and >50. The ESV intensity of the GPUA showed a distribution trend of low in the central regions and high in the surrounding regions. The grid in the low ESV intensity range is mainly caused by the relatively uniform spatial distribution and high proportion of farmland and construction land, accounting for about 22.45% of the total grid area. The Guanzhong Plain area is suitable for agricultural development, the scope of urban settlements expands, the ecological and environmental effects are low, and negative effects are difficult to avoid. Areas at high altitude, slopes, and Yellow River and Weihe River basin areas are rich in ecological resources, their ecosystems are more complete, and the positive ecological environment effect is significant; moreover, the value of services available is high, accounting for about 44.71% and 23.61% in zones II and III, respectively, and for about 9.22% in the high ESV intensity zone. Overall, the distribution range of low-intensity areas has been gradually shrinking and transforming into II intervals, and the overall ESV intensity of the GPUA has increased.

To examine the interaction between carbon emission intensity and ESV intensity, the spatial correlation pattern was identified by bivariate spatial autocorrelation. The Moran’s I index was negative in all years with a p-value of less than 0.001 (shown in Table 9). This result indicates the existence of a negative spatial correlation between land-use carbon emission intensity and ESV intensity in the study area, i.e., the ESV intensity followed a decreasing trend with increasing land-use carbon emission intensity.

The bivariate local autocorrelation LISA agglomeration map shows that between 2000 and 2020, the high and high agglomerations are concentrated in the Yellow River and Weihe River basins; low and high agglomerations are mainly concentrated in the surrounding areas with higher slopes and topography, such as the Qinling Mountains in the south and the Loess Plateau area in the west. These areas have low land opening intensities, low carbon emission intensities, better ecological backgrounds, and higher values of ecosystem services available. These factors are consistent with the regional distribution pattern of low carbon emission intensity and high ESV intensity, and present roughly the same spatial distribution trend as woodland and grassland. Figure 5 shows that over the study period, the number of high and high aggregation areas in the Guanzhong Plain increased and the area expanded, while the area of high and low aggregation areas in the northwestern part of the study area decreased. This is probably the result of the increasing intensity of ecosystem service provisioning caused by the structural adjustment of land use functions and the gradual implementation of the policy of returning farmland to forests and grassland in loess hills and gullies. The construction land covers a small area, has a high carbon emission coefficient, and is a carbon source, resulting in high carbon emission intensity in local areas. The coefficient is much higher than that of the surrounding farmland, and thus, the grid area with a high proportion of cultivated land area together with the construction land forms a high–low agglomeration area. In the marginal areas of urban agglomeration zoning, low–low agglomeration areas show a scattered distribution.

Land use is the most direct expression of the interaction between human activities and natural ecosystems. The structure of land use is not only an important part of basic land use units, quantitative characteristics, and spatial layout, but also the core basis of regional ecosystem macro control. Land use change affects the carbon cycle process of terrestrial ecosystems by changing the original land cover type and the socio-economic activities it carries. With the acceleration of the urbanization rate of the GPUA, the development and utilization of land are gradually enhanced, and construction land is the most important application site for fossil fuels. The significant expansion of the areas of construction land causes a shrinkage of the carbon sink area and an expansion of the carbon source area, which directly or indirectly causes changes in carbon emissions and carbon sequestration and affects the supply capacity of ecosystem services. The total carbon emissions from land use in the study area increased significantly from 2000 to 2020, and the carbon emissions from construction land exceeded net carbon emissions. The main reason for this phenomenon is that the development of rapid urbanization has led to significant changes in land use structure. The GPUA, with Xi’an at the core, is rapidly promoting urban construction and social development, thus facilitating the process of industrial construction and service construction. The output value of secondary and tertiary industries accounts for about 90% of the regional GDP. The carbon sequestration capacities of woodland and grassland play an important role in ecological improvement and protection, and these land use types account for a large proportion of the land use area. The maximum carbon emission intensity within the study area decreased from 2000 to 2020, which is generally consistent with the findings of previous related studies (Cao et al., 2021; Zhang et al., 2022b; Wu et al., 2022; Zhao et al., 2023).

From 2000 to 2020, the total value of ESV shows a fluctuating growth trend of “rapid increase–rapid decrease–slow decrease–rapid increase”, including two significant change periods from 2000 to 2010 and from 2015 to 2020. This differs from the results of previous research related to the Guanzhong Plain Urban Agglomeration (Yang et al., 2022b). The higher ESV values were distributed in the southern and western and northwestern edges of the study area, in the Qinba Mountains and along the Yellow River, enabled by the positive effect of forest cover and water conservation on ESV intensity. The rapid increase in the total ESV value in 2005 may have been due to the transfer of forest land area as the largest increase during the study period. The implementation of the project of returning farmland to forest and grassland and the erosion control project provide a good start for the restoration of ecosystem service functions. In 2018 China promulgated The Guanzhong Plain Urban Agglomeration Development Plan, which specifies ecological environmental protection as a task and prerequisite for the construction of city clusters. This provides policy support for the rapid growth of ecosystem service values and the implementation of ecological protection and restoration projects from 2015 to 2020 (Guo et al., 2022; Wang et al., 2022).

The spatial agglomeration of carbon emission intensity and ESV intensity shows a significant negative correlation effect. On production land and living land such as arable land and construction land, the central region of the GPUA generates a large amount of greenhouse gas pollution and emissions due to carrying a large number of human production and living behaviors, and clustering the development and utilization of resources in agriculture, industry and services. Human behavior is also a consumption of ecosystem functions, so it makes the value provided by terrestrial ecosystems impaired, showing a spatial agglomeration distribution pattern of high carbon emission intensity–low ESV intensity. In the southwestern part of the GPUA, the southern hinterland, and part of the northeastern part of the GPUA, i.e., the area with a large distribution of ecological land such as forests and grasslands, a low carbon emission intensity-high ESV intensity distribution pattern is shown. This is in line with the functions of natural ecology such as woodland and grassland in performing soil and water conservation, climate regulation, or other functions. At the same time, a series of man-made ecological projects, such as the government’s establishment of the Qinling Protected Area and the introduction of the Outline of the Plan for the Ecological Protection and High-Quality Development of the Yellow River Basin, have made an important contribution to the ecological environment of the region. How to reduce the total amount of carbon emissions from land use and carbon sources and improve the level of urban ESV is the breakthrough direction for the future green development of the city and the realization of long-term synergistic development of ecology and economy.

The above changes in land use structure and function yield corresponding changes in total carbon emissions and ESV, leading to ecosystem instability. However, it is worth noting that farmland has a weak ecological carrying capacity, and when making changes to utilization strategies, it is important to pay attention to regional food security and farmland red line restrictions, to not develop and transfer out uncontrollably. The GPUA should prioritize the ecological security pattern in the development process, and the governments of Gansu, Shaanxi, and Shanxi provinces should do a good job of coordinating the division of labor based on the consideration of ecological value coordinated with social and economic value. The three provinces should work together to protect the ecological barrier areas, ensure the smooth flow of ecological corridors, and consolidate the green foundation for the high-quality development of the urban agglomeration (Chen et al., 2022; Yang et al., 2022c). This study complements the relevant research results on carbon emissions and ESV under land use change, and using the GPUA as a case study, proposes strategies to promote integrated regional ecological-economic development, which helps to improve the practicality and completeness of the research results.