Carbon Metabolism in Urban “Production–Living–Ecological” Space Based on Ecological Network Analysis

: To understand the changing pattern of urban carbon metabolism from the perspective of urban “production–living–ecological” (PLE) space, taking Suzhou City as an example, this study constructed a carbon metabolic network model in urban PLE space, analyzed the changes of horizontal carbon ﬂow, and evaluated the comprehensive effect of the PLE space changes using the ecological network analysis method. The results showed that the total carbon sequestration showed a ﬂuctuating change of increasing and then decreasing, while the total carbon emissions grew dramatically. Production spaces were the key nodes for the generation of horizontal carbon ﬂow. The exploitation relationship was the dominant ecological relationship in the network, the mutualism relationship was abundant from 2005 to 2010 and gradually decreased from 2010 to 2018, and the frequency of competition relationship appeared gradually increased. The ecological network hierarchy evolved from an irregular shape dominated by primary consumers in 2000–2005 to a pyramidal shape dominated by producers in 2010–2018 at the driving weight end, and the pull weight showed a declining trend, with pull weight of producers increasing from 1.72% to 24.33%. The results can provide a theoretical basis for planning adjustments to the city’s PLE space structure to achieve low-carbon goals.


Introduction
Urban areas, which account for about 2% of the global land area, contribute nearly 75% of global carbon emissions and are one of the major causes of global warming [1]. Land use change in the urban "production-living-ecological" (PLE) space is an important source of urban carbon emissions [2]. Land use change, as the dual body of both carbon emissions and sequestration, directly transforms or indirectly affects the carbon emission and sequestration processes between the terrestrial ecosystem and the atmosphere by changing the land use type, function, and structure [3,4]. Therefore, it is necessary to understand the linkage between urban carbon and the change of PLE space to identify the critical path of carbon emissions reduction through land use adjustment in spatial planning [5]. By constructing a carbon metabolic network model in urban PLE space based on land use change, we can simulate and account for the carbon metabolism flow in cities and analyze the ecological relationship and hierarchical structure in the process of carbon metabolism, which will firstly help to reflect the interrelationship and coordination between the urban economy, environment, and resources [6]. Secondly, it can measure the changes in carbon metabolism flow, metabolism anomalies, and the rationality of urban PLE space structure, to provide a scientific reference for regulating the process of urban different carbon metabolic subjects due to the transition of PLE space. It is more conducive to policymakers carrying out land planning from a socio-economic macro perspective to promote urban carbon neutrality, taking into account the overall characteristics of the region at the same time.
With the release of the outline of the integrated regional development in the Yangtze River Delta, the integration of the Yangtze River Delta Economic Zone is elevated to a national strategy, emphasizing the importance of realizing sustainable urban development by means of territorial space optimization [26,27]. As an important central city in the Yangtze River Delta region, Suzhou is a key node for promoting the integrated ecological and green development of the Yangtze River Delta [28]. Therefore, this study focused on the complex impacts of PLE space changes on urban carbon metabolism, taking Suzhou City as an example. We constructed a carbon metabolic network model based on the classification of PLE space, accounted for the carbon emissions and sequestration of each node in the network, and analyzed the changes in horizontal carbon flow among the nodes in the network with the land use transition matrix. Then we evaluated the influences of the change in urban PLE space on the carbon metabolic network by ENA. The results of the study can provide a scientific reference for optimizing urban spatial layout and promoting low-carbon urban land use.

Study Area
Suzhou City (119 • 55 -121 • 20 E, 30 • 47 -32 • 02 N) is located in the southeast of Jiangsu Province and on the east bank of Taihu Lake (Figure 1), with a total area of 8657.32 km 2 . The topography in Suzhou is flat and low, with plains accounting for 54.8% of the total area, numerous lakes, and dense rivers. As of 2018, Suzhou had a resident population of 10.72 million people, an urbanization rate of 76.05%, and a GDP of CNY 1.86 trillion (USD 0.27 trillion), accounting for 20.08% of the Jiangsu Province. With the rapid development of Suzhou's economy and the gradual expansion of its industrial scale, the total carbon emissions from socio-economic activities have been increasing, and the ecological spaces with carbon sink function such as forest and water have been compressed due to the continuous promotion of the urbanization process, which has posed a great challenge for Suzhou in achieving the "carbon peak and carbon neutral" strategic goal. How to optimize the spatial layout of the city and explore long-term and profound emissions reduction paths while ensuring scientific development has become the key to its low-carbon transformation.

Data Sources
The socio-economic data used in this study, including agriculture, industry, population, energy consumption data and so on, were derived from the Suzhou Statistical Yearbook, Jiangsu Statistical Yearbook, China Transport Statistical Yearbook, China Energy Statistical Yearbook, and Suzhou National Economic and Social Development Statistical Bulletin for the corresponding years, etc. (http://tj.jiangsu.gov.cn/, accessed on 22 April 2022). Land use data with a resolution of 30 m in 2000, 2005, 2010, and 2018, including 6 primary land categories such as cultivated land and forest land and 25 secondary land categories, were provided by the Chinese Academy of Sciences Resource and Environmental Science Data Center (http://www.resdc.cn, accessed on 22 April 2022). Based on the PLE dominant function of land use, and referring to the research of Yang et at. [29], we merged the basic land use types taking their dominant function into account and established a classification system of the PLE space (Table 1). The carbon metabolic network model in urban PLE space constructed in this study took into account the carbon exchange of network nodes between the biosphere and the atmosphere, as well as the potential carbon flow in the horizontal direction resulting from the transitions between the PLE space types ( Figure 2). The manner of carbon exchange between each network node and the atmosphere determines its carbon metabolism function (Table 1). In this study, IN, UR, and RU were the socio-economic components of the urban system, which were the main sources of urban carbon emissions and possessed the carbon source function; G, F, and W were the natural compartments of the urban system, which were responsible for the carbon sequestration process in urban carbon metabolism and possessed the carbon sink function; CU had dual attributes and assumed both functions. B was hardly used as the carrier of socio-economic activities, so it assumed no carbon exchange with the atmosphere. Note: IN, industrial production space; UR, urban living space; RU, rural living space; G, grassland ecological space; F, forest ecological space; W, water ecological space; CU, agricultural production space; B, other ecological space; Red arrow, carbon source function; Green arrow, carbon sink function.

Carbon Emissions Accounting
Carbon emissions from CU were mainly from agricultural production activities, livestock respiration, and crop cultivation [30]. The agricultural production activities mainly included crop sowing, agricultural machinery activities, irrigation process, and fertilizer application; the livestock types mainly included pigs, cattle, and sheep according to the specific characteristics of the study area; the agricultural crop in this study was mainly rice. Its total carbon emissions ECU was specifically calculated as followed: where E CU was the total carbon emissions from CU; E a was the carbon emissions from agricultural production activities; E b was the carbon emissions from livestock respiration; E c was the carbon emissions from crop respiration; K 1 was the total sown area of crops; K 2 was the total power of agricultural machinery; K 3 was the irrigated area; K 4 was the fertilizer application; K 5 , K 6 , and K 7 were the annual pig, cattle, and sheep rearing in the study area, respectively; K 8 was the rice cultivation area; f i was the carbon emissions conversion coefficient corresponding to each agricultural activity (Supplementary Table S1). IN included industrial, mining, and transportation construction land, so its carbon emissions were also divided into industrial energy consumption carbon emissions and transportation carbon emissions. Its total carbon emissions E IN were specifically calculated as follows: E e = K 9 f 9 +K 10 f 10 +K 11 f 11 +K 12 f 12 +K 13 f 13 +K 14 f 14 (7) where E IN was the total carbon emissions from IN; E d was the carbon emissions from industrial energy consumption; E e was the carbon emissions from transportation; C i was the consumption of the i th energy source; d i was the carbon emissions coefficient of the i th energy source; K 9 , K 10 , K 11 , and K 12 were the annual running kilometers of buses, cabs, private cars, and motorcycles, respectively; K 13 and K 14 were the annual passenger and freight volume of roads and railroads, respectively; f i was the carbon emissions coefficient of each transportation type (Supplementary Table S1). Carbon emissions from UR were mainly the carbon emissions of urban residents' respiration and the carbon generated by urban residents' electricity consumption. Similarly, the carbon emissions from RU mainly came from the carbon emissions of rural residents' respiration and the carbon generated by rural residents' electricity consumption. Their total carbon emissions E UR and E RU were specifically calculated as followed: where E UR was total carbon emissions from UR; E RU was total carbon emissions from RU; K 15 was urban living electricity consumption; K 16 was urban resident population; K 17 was rural living electricity consumption; K 18 was rural resident population; f i was carbon emissions coefficient of each index (Supplementary Table S1).

Carbon Sink Accounting
CU, G, F, and W can absorb carbon from the atmosphere and sequester it. In this study, the carbon sink coefficient was used to calculate the amount of carbon sink by combining the area of each space that possessed the carbon sink function. The specific calculation formula was as follows: where V i was the carbon sequestration of the i th type of PLE space; S i was the area of the i th type of PLE space; k i was the carbon sink coefficient of the i th type of PLE space (Supplementary Table S2). In order to improve the accuracy of calculating the carbon sequestration of each space, all secondary land types were considered in this study.

Horizontal Carbon Flow Accounting
In this study, the net carbon emissions or sequestration per unit area (i.e., carbon metabolism density) of each PLE space type was used to characterize the carbon metabolism capacity of the PLE space type. The change of PLE space type led to the change in its carbon metabolism density, and the horizontal carbon flow among network nodes can be quantified according to the change of area of each PLE space type, which was reflected by the land use transfer matrix [31]. The specific calculation formula was as follows: where f ij was the horizontal carbon flow from PLE space type j to i; ∆W was the difference in carbon metabolism density between PLE space type i and j; ∆S was the area transferred from PLE space type j to i; S i and S j were the areas of PLE space type i and j, respectively; V i and V j were the net carbon exchange between PLE space type i/j and the atmosphere. Because carbon metabolism capacity included carbon emissions and sequestration, there was positive and negative horizontal carbon flow among the PLE spaces in the network. If ∆W > 0, it indicated a positive carbon flow, which represented a decrease in carbon emissions or an increase in carbon sequestration, and helped to alleviate the carbon metabolism imbalance of the urban system; if ∆W < 0, it indicated a negative carbon flow, which represented an increase in carbon emissions or a decrease in carbon sequestration, aggravating the carbon metabolism imbalance in the urban system. Based on the above equation, the carbon flow among PLE spaces was calculated to constitute the direct carbon flow matrix F in the horizontal direction.

Ecological Network Analysis (ENA) Method
The ENA method originated from the analysis of currency circulation in economic problems. In 1973, Hannon first applied the input-output model in economics to analyze the energy flow in an ecosystem and the hierarchical structure of the entire ecosystem, forming the prototype of the ENA method [32]. As research progressed, the ENA method was widely used to analyze the structure of ecosystems and various ecological attributes, and the quantitative models used were also enriched [33][34][35]. In this study, we used the throughflow, utility, and structure analysis of the ENA method to quantitatively assess the carbon throughflow of each node, the ecological relationships between each PLE space type, and the hierarchical structure of the network, to reflect the impact of urban PLE space changes on urban carbon metabolism.

1.
Throughflow analysis Each PLE space type had the potential for absorbing carbon or releasing carbon from the environment. In this study, the state variable x k was defined to reflect the change in the carbon stock of PLE space type k to balance the total outflow and inflow of carbon in the network [36]. If x k > 0, then the carbon throughflow T k of the space was equal to the sum of all carbon flow out of the space plus the state variable; if x k < 0, then the carbon throughflow T k of the space was equal to the sum of all carbon flow into the space minus the state variable [19,37]. The specific calculation of throughflow T k was as follows: where f ik was the horizontal carbon flow from PLE space type k to i; y k was the carbon flow released from PLE space type k to the environment; (x k ) + represented the increase of carbon stock; f kj was the horizontal carbon flow from PLE space type j to k; z k was the carbon flow absorbed from the environment by PLE space type k; (x k ) − represented the decrease of carbon stock; n was the number of network nodes, and n was 8 in this study.

Utility analysis
Based on the direct carbon flow matrix F calculated above and the throughflow analysis results, with reference to the results obtained in previous research [38,39], the utility analysis method was used to determine the ecological relationships between each PLE space type. The calculation steps were as follows: 1 Based on the carbon throughflow of each PLE space type, the direct carbon flow matrix F was normalized to compose the direct utility matrix D, which reflected the direct utility of the horizontal carbon flow between each node. 2 Based on the direct utility matrix D, the overall utility matrix U was calculated to fully reflect the direct or indirect utility relationships of the horizontal carbon flow. 3 Based on the properties of positive and negative of the elements in the overall utility matrix U, the ecological relationships in the network were determined; 4 the mutualism index MI and the synergy index SI were defined to quantify the comprehensive impact of the PLE space changes on the urban carbon metabolic network. The specific calculation was as follows: where d ij was the net utility of carbon flow from PLE space type j to i; u ij represented the elements in the overall utility matrix U; I represented the unit matrix; the superscript m of the direct utility matrix D represented the m paths that the carbon flow had to go through to reach the final PLE space type; N +/− represented the number of positive/negative elements in the overall utility matrix U. If MI > 1, it indicated that the PLE space changes had a positive effect on the urban carbon metabolic network, and the larger the MI, the stronger the positive effect. If MI < 1, it indicated that the PLE space changes had a negative effect on the urban carbon metabolic network, and the larger the MI, the stronger the negative effect. SI was the sum of the elements in the overall utility matrix U. A larger SI meant that the positive effect of horizontal carbon flow was more obvious, and the urban carbon metabolic network was more coordinated.
Based on the properties of positive and negative of elements in the overall utility matrix U, the ecological relationships between each PLE space type in the network can be determined ( Table 2). Although there were theoretically nine ecological relationships, four of them were more common (mutualism, competition, exploitation, and control) [40,41]. The control and exploitation relationships were similar in nature, so we combined them into a single category (exploitation relationship), both indicating that one PLE space type obtained more carbon storage in the process of carbon transfer, while the other space type lost carbon storage. The mutualism relationship indicated that both PLE space types obtained more carbon storage as a result of horizontal carbon flow in the network, while the competition relationship indicated the contrary. Both were the outcome of indirect effects of carbon flow on each PLE space type in the network. Table 2. Ecological relationships between PLE space types in the network.

Structure Analysis
The utility analysis can reflect the ecological relationships in the network, but it was difficult to reflect the overall structure of the urban carbon metabolic network, and the level and specific role of each node in the network. Therefore, the biological concepts of producer, consumer, and decomposer were introduced [42], referring to previous research [43,44], to analyze the hierarchical structure of the carbon metabolic network model in urban PLE space (Figure 3), and quantify the capacity of each node to absorb and export horizontal carbon flow to reflect its importance in the network.
Grassland (G), forest (F), and water (W) ecological space, as producers in the network structure, were located at the bottom of network, providing natural resources to the upper layer, and maintaining the carbon balance of the whole carbon metabolic network through carbon sinks. Agricultural production space (CU) and rural living space (RU), as primary consumers, were at the junction of natural and socio-economic compartments, and used the natural resources provided by producers for agricultural production and primary agricultural products processing in order to provide products to more advanced consumers, and their carbon metabolic capacity was also relatively low. Industrial production space (IN), as secondary consumer, carried out secondary processing of agricultural products and transported them to senior consumers, which was an important node in the generation of horizontal carbon flow. Urban living space (UR) acted as a senior consumer and used these products and natural resources to improve human living standards and develop the economy. Other ecological space (B) acted as a decomposer and was highly susceptible to human activities, often serving as a "transfer station" for PLE space type change, e.g., bare space may be formed by the conversion of abandoned buildings in living and production space, but it can be transformed into ecological space through a series of treatment measures. Thus, its ability to release or absorb carbon flow was very low. Based on the results of the throughflow analysis above, we can obtain the structural utility matrix Y. By calculating the elements in Y, we can derive the driving weight and pull weight of each PLE space type in the network. The driving weight indicated the ability of the space type to influence other space types by exporting horizontal carbon flow, while the pull weight represented the ability of the space type to receive carbon flow [45]. The specific calculation formula was as follows [46]: where I was the unit matrix; diag(T) was the diagonal matrix composed of throughflow T i ; ; f ij was the horizontal carbon flow from PLE space type j to i; T j was the carbon throughflow of PLE space type j; W j was the driving weight of PLE space type j; W i was the pull weight of PLE space type i.

Carbon Emissions and Carbon Sequestration Changes
The results of carbon emissions and sequestration accounting in Suzhou from 2000 to 2018 are shown in Table 3. The total carbon sequestration in Suzhou showed a fluctuating change from rising to falling during the 18 years. From 2000 to 2005, the total carbon sequestration increased from 12.40 × 10 4 t C to 12.51 × 10 4 t C. The carbon sequestration in CU decreased by 9.17%, while the carbon sequestration in W increased by 12%, accounting for 90.05% of the total. From 2005 to 2010, the total carbon sequestration rapidly declined, with an average annual decrease of 4.23%, and the carbon sequestration of the PLE space which possessed the carbon sink function experienced different degrees of decline. From 2010 to 2018, the decline in total carbon sequestration tended to level off, with a rate of 4.93%. The carbon sequestration in CU and F had not changed much, while the carbon sequestration of W decreased 6.43%, which is the most important space type for carbon sink.

Horizontal Carbon Flow Changes
The results of the horizontal carbon flow accounting are shown in Supplementary Tables S3-S5. The total amount of negative carbon flow in Suzhou shows an increase over the three study periods. From 2000 to 2005, the total negative carbon flow was 1.22 × 10 6 t C, mainly from the transition of CU→IN and CU→UR, accounting for 78.22% and 10.07% of the total, respectively. From 2005 to 2010, the total negative carbon flow was 3.26 × 10 6 t C, which was 2.68 times of the previous study period, and the transition of CU→IN and W→IN contributed 59.25% and 13.76% of the total, respectively. From 2010 to 2018, the total negative carbon flow rose further, to 7.26 × 10 6 t C, mainly from CU→IN, W→IN, and G→IN, which accounted for 45.35%, 40.34%, and 6.90% of the total, respectively.
The spatial distribution of negative carbon flow is shown in Figure 4. In order to reveal the multifaceted characteristics of the spatial changes of carbon flow more comprehensively, Standard deviational ellipse (SDE), a spatial statistical method widely used in many fields such as sociology and geology was introduced (Supplementary Figure S1).  The SDE analysis showed that the positive carbon flow appeared along the north-south direction as a whole, similar to the negative flow, but the distribution range and the mean center shifted more southward than the negative. The short axis of the SDE did not change much in each study period, the area of SDE showed a fluctuation trend of first becoming smaller and then larger, and the azimuth angle changed more from 21 • →26 • →−9 • . It showed that the transition of IN in Kunshan and Gusu had a large traction effect from 2005 to 2010. During the following eight years, the areas of positive carbon flow decreased, and the distribution pattern tended to be more fragmented. Huqiu and Wuzhong had played a positive role in the environmental management of Taihu Lake, which contributed to the reverse deflection of the azimuth of SDE.

The Ecological Relationships between PLE Spaces
From 2000 to 2018, the proportion of each ecological relationship between PLE spaces in Suzhou had been changing (see Figure 6a, Supplementary Tables S6-S8). From 2000 to 2005, the exploitation relationship was absolutely dominant and accounted for 75% of the total, and each PLE space type was fiercely plundering each other's carbon stock. Although that had unexpectedly improved the overall utility of carbon flow in the network, the development mode in which one side gained and one side lost was not conducive to the long-term carbon health of the city. The competition relationship accounted for 17.86% of the total, mainly existing in such artificial compartments as IN, RU, and UR. The proportion of the mutualism relationship was only 7.14%, all of which were related to B, reflecting their neutralizing role as decomposer in the urban ecological process. From 2010 to 2018, the proportion of the exploitation relationship increased slightly to 53.57%, the proportion of the competition relationship increased to 32.14%, while the proportion of the mutualism relationship decreased to 14.29%. During 2000-2018, the spatial distribution of ecological relationships between PLE spaces in Suzhou is shown in Figure 7. From 2000 to 2005, exploitation relationship was distributed in all districts and counties, being was the highest in Xiangcheng, Gusu, and Kunshan, accounting for 12.68%, 11.02%, and 11.02% of their regional areas, respectively. The spatial distribution of the competition relationship was more dispersed, with Kunshan and Xiangcheng still dominating, accounting for 4.15% and 3.18% of their regional areas. The distribution of the mutualism relationship was relatively less dispersed, mainly concentrated in Wuzhong. From 2005 to 2010, the changes in PLE space intensified and the spatial distribution of each ecological relationship became intersecting. The main distribution regions of the exploitation relationship had not changed. The competition relationship increased significantly in all districts and counties, especially in Xiangcheng, accounting for 18.22% of its regional area, mainly from IN and UR. The mutualism relationship was mainly distributed in central Suzhou and along Taihu Lake, occurring mainly from UR and W. From 2010 to 2018, the distribution area of each ecological relationship decreased. The exploitation relationship increased in Taicang, where the ratio of the area increased to 9.53%. The competition relationship was more apparent in Zhangjiagang, Wujiang, and Changshu. The spatial distribution of the mutualism relationship was becoming less scattered and mainly concentrated in Taicang.
The overall utility matrix U was used to calculate the mutualism index MI of the carbon metabolic network in urban PLE space in Suzhou for each study period (see Figure 6b). The MI for 2000-2005 and 2005-2010 were greater than 1, which indicated that the spatial transition of PLE space had a positive impact on the network. From 2005 to 2010, MI attained the highest value of 1.13, which was mainly due to the greater prominence of the mutualism relationship between UR, G, F, and B. From 2010 to 2018, MI decreased to 0.94 and was less than 1, indicating that the spatial transition of PLE space had a negative impact on the network. The synergy index SI was 6.69→5.81→5.17 in the three study periods, showing a decreasing trend, mainly due to the increase in the competition relationship in the network, which generated intense competition for carbon storage. The self-feedback of each PLE space type for horizontal carbon flow was reduced, leading to a certain negative effect on the network.

Hierarchical Structure of Urban Carbon Metabolism Network
The driving and pull weights of each PLE space type in the network during 2000-2018 are shown in Figure 8. At the driving end, the network hierarchy from 2000 to 2005 showed an irregular state, with CU as primary consumer taking the largest driving weight (47.65%), followed by IN as secondary consumer (41.79%), while F, G, and W as producers only contributed 1.93%. The network experienced a distinct production shortage, with yield struggling to meet the demand for urban carbon metabolism. From 2005 to 2010, this shortage was alleviated, and the driving weight of producers increased to 8.25%. As primary consumer, the driving weight of CU decreased to 31.41%. In order to maintain the stability of the whole network, primary consumers assumed part of the ecological functions of producers. From 2010 to 2018, the network hierarchy gradually evolved into a pyramidal structure, with the driving weight of producers further increasing, especially the driving weight of W, rising to 35.71%. The weight of the secondary consumer (IN) decreased to 10.74%, and the weight of the senior consumer (UR) decreased to 1.73%.
At the pull end, from 2000-2005 the network hierarchy showed an unhealthy structure similar to an ellipse. The secondary consumer (IN) contributed the largest pull weight (61.91%) and absorbed a large amount of carbon flow for development, followed by CU as the primary consumer (24.22%). The producers and decomposers bore very low pull weight during this period, which was not conducive to the realization of the urban carbon balance. This phenomenon changed during 2005-2010, with a decrease in the pull weight of CU (24.22%→14.54%), an increase in the weight of the senior consumer (UR) (8.27%→21.46%), and the ability of the secondary consumer (IN) to absorb carbon flow decreased slightly (−7.04%). The problem of smaller carbon flow absorption by producers was slightly alleviated, and the producers' pull weight increased to 4.53%, but it was still difficult to maintain the carbon balance of the whole network. From 2010 to 2018, the pull weight of the secondary consumer (IN) and the senior consumer (UR) further decreased, and the producers' weight increased further (+19.80%).

The Changes in PLE Space Types and Carbon Flow
Urban space is a complex multi-layered spatial system consisting of various elements such as geographical features, human resources, economic conditions, and policies and management [47]. The transition of PLE space reflects the changes in the relationship between the natural environment and socio-economic factors from the perspective of urban land use, which is the main external expression of urban spatial evolution [48]. The transformation of PLE space brought about the transition of PLE dominant functions of land. It is a process of orderly allocation and spatial reconfiguration of limited urban land resources among different dominant functions [49], which greatly affects the carbon cycle of the urban space and the carbon metabolism efficiency of the regional land system. Yang et al. [29] studied the eco-environmental effects and changes in terrestrial carbon stocks caused by land transition in the Beijing-Tianjin-Hebei (BTH) metropolitan region based on the classification of PLE dominant functions, further confirming the practicality of this classification method. Therefore, we adopt a similar classification method based on PLE dominant functions to identify the nodes in the urban carbon metabolic network, considering the multi-functionality of urban land with regard to socio-economic aspects. It was more appropriate to the evolution law of the urban land system and will help to formulate more low-carbon land management policies adapted to each stage of urban socio-economic development.
Suzhou is a representative region of China's economic transformation and development, reflecting the change process of the regional development path since China's reform and opening up [28]. It had undergone phased urbanization and industrialization from 2000 to 2018 [50], and the layout of PLE space had changed dramatically. In this study, we measured carbon emissions and sequestration at four time points of each PLE space type based on statistical data, and calculated the horizontal carbon flow in the carbon metabolic network in three study periods. The negative carbon flow in Suzhou City mainly came from the transition from CU to IN, accounting for more than 40% of the total in each study period. This was similar to the research on Nanjing City by Zhao et al. [51]. And it was also appropriate to the characteristics of Suzhou's industrial structure adjustment, which showed a shift from agriculture to industry and services as the main driving force of the city's economic development. However, through further research, we found that the proportion of negative carbon flow generated from the transition from CU to IN was decreasing over the three study periods, and the spatial pattern of negative carbon flow changed from a large concentration of high-value regions in 2005-2010 to a fragmented distribution in 2010-2018. This was mainly because Suzhou's economic development had entered a new stage, the expansion of UR in core urban regions such as Kunshan and Gusu had gradually stabilized, the transformation of large-scale agricultural to non-agricultural land had come to an end, and the restrictions on the expansion of IN had gradually strengthened. The positive carbon flow was mainly generated from the transition within artificial compartments in the first two study periods, especially from 2005 to 2010 when the conversion of IN to UR generated a large amount of positive carbon flow; this reflected the extremely high carbon emissions density of IN on the one hand, and confirmed the key role of industrial land restrictions in promoting local carbon balance on the other [52]. During 2010-2018, due to the vigorous promotion of policies such as the Suzhou Low Carbon Development Plan (2010-2020), the transition of production and living space to ecological space dominated the positive carbon flow. Therefore, scientific spatial planning was indispensable in low-carbon land management [4].

The ecological Relationships between PLE Space and Urban Expansion
The utility analysis in the ENA method had revealed the ecological relationships formed in the process of carbon transfer in urban PLE space. From 2000 to 2005, the exploitation relationship was absolutely dominant, mainly reflected in the deprivation of ecological space and CU by artificial compartments, especially in Kunshan and Gusu, which caused the rapid increase in carbon emissions in the central part of Suzhou. The reason for this deprivation was the development orientation at the early stage of Suzhou's urbanization, which was dominated by small towns [53]. This orientation led to the expansion of towns in a chaotic way, which reduced the efficiency of spatial utilization and formed an inefficient and disordered spatial pattern of urbanization, which was not conducive to the effective control of urban carbon emissions. Similar conclusions had been reached in some studies from Switzerland [54] and Mexico [55], highlighting the important influence of urban spatial form optimization on carbon emissions reduction.
From 2005 to 2010, Suzhou experienced the most rapid transition of PLE space, gradually getting rid of the "Southern Jiangsu Pattern", which was dominated by the expansion of small towns and industrial development in the countryside [56], and presented a larger scale and more concentrated spatial characteristics of urbanization. Given these characteristics, on the one hand, the orderly expansion of cities further compressed the areas of agricultural and ecological space and forced them to compete with each other in a smaller space, coupling with China's unique policy of protecting farmland, i.e., compensating for the loss of farmland due to urban expansion by occupying other space types to maintain the dynamic balance of farmland quantity [57]. This policy further exacerbated the appropriation of agricultural land for ecological space, and the increase in carbon emissions from this occupation was noteworthy due to the large differences in carbon metabolic capacity. A study from the United States showed that slash-and-burn cultivation resulted in a loss of 436. 8 TgC/year of carbon sequestration in woodlands [58]. Similar phenomena had occurred in Hangzhou [46], which shared many similar urban characteristics with Suzhou. On the other hand, it was also the government-led expansion of cities that enabled some strong environmental protection measures to be effectively implemented from the top down, resulting in some transient appearance of the mutualism relationship, not only within ecological space but also between UR and ecological space. This confirmed the possibility of promoting a healthy carbon cycle by planning urban green landscapes and protecting key ecological space [59].
The rate of urban spatial and functional expansion in Suzhou had slowed down from 2010 to 2018, and the peri-urban areas around the core urban regions had developed into fully urbanized regions with an enhanced stability of the ecological relationship. A noteworthy development was the disappearance of the mutualism relationship between key nodes in network. Although many ecological red line protection plans had been implemented in Suzhou during this time period [60], and the rate of PLE space conversion had slowed down, a certain number of peri-urban or rural areas still existed in certain districts, which were far away from main city areas [61]. How to intensively and efficiently turn these areas into UR by means of spatial planning and moderately increased urban population density, in order to improve energy consumption efficiency and reduce carbon emissions, was the main concern in Suzhou's future low-carbon development.
In our study, the impacts of urban PLE space transition on the carbon metabolic network were quantified by the mutualism index MI and the synergy index SI. The MI reached a maximum of 1.13 from 2005 to 2010, indicating that the transition of PLE space had a positive impact on the urban carbon metabolic network during this period. But it was still far below the 1.91 found by Chen and Chen in the carbon metabolic system of Vienna [39]. This difference may be the result of the difference in the maturity of the two cities: Vienna is an old and relatively mature city with stable development [62], while Suzhou has undergone an unparalleled rapid urbanization process during the past 18 years, and the city lacked stability in this phase of rapid change. Moreover, the two indices showed different trends in each study period. The main reason was that MI reflected primarily the degree of influence of the ecological relationships between each PLE space type on the carbon metabolic network, and the calculation was based on the positive and negative properties of each element in the overall utility matrix [39]; thus, the competition relationship and the mutualism relationship were crucial to MI, and the number of different network nodes also had a greater influence on it. Zhang et al. [18] and Xia et al. [19] used the same method to analyze the mutualism index of urban carbon metabolism in Beijing with 8 and 18 network nodes, and the mean difference was as high as 0.94. However, SI, originating from the sum of the values of elements in the overall utility matrix, could quantify the utility of carbon flow in the network more finely and reduce the errors due to the different ways of node classification. The brief emergence of the partial mutualism relationship in Suzhou from 2005 to 2010 improved the MI of the network, but in the subsequent study period, the mutualism relationship shifted back to exploitation relationship, indicating the fragility of this brief mutualism relationship. However, SI captured this difference well and showed a decreasing trend in the three study periods. Therefore, when analyzing the overall utility of carbon flow in the network, the results of both indices should be combined to reflect the ecological utility of carbon flow more comprehensively.

Policy Recommendations for Low-Carbon Urban Development
The low-carbon city has been a hot concept in China in recent years [63]. The National New Urbanization Plan (2014-2020) released by the Chinese government in 2014 emphasized that urban development should be people-oriented and environment-friendly [64]. Although people have reached a consensus on the concept of low-carbon, its realization path still needs to be explored in practice because of the different PLE space patterns and socio-economic development stages of Chinese cities. As a representative highly urbanized region in China, the rapid urbanization in Suzhou over the past 18 years had caused a major increase in carbon emissions from the land type with carbon source function, further negative ecological relationships in the carbon metabolic network, and the imbalance of the hierarchical structure. The essential reason was that the economy, society, and ecological system interacting with each other in the course of the rapid development [65] resulted in the intense conflict between PLE functions of land, and the formation of a less compatible spatial pattern of PLE types [66]. The effective combination of the ENA method and the carbon metabolic network in urban PLE space analysis framework can provide a more social and multifaceted perspective. By exploring the comprehensive impact of the transition of PLE space on the ecological relationship and hierarchical structure of the network, more targeted low-carbon land management policies can be proposed: (1) For fully urbanized regions such as Kunshan and Gusu, in the central part of the study area, where the senior consumer (UR) was the core, on the one hand, it is necessary to restrain the scale of urban space, promote the integrated development of land and transportation, and protect key ecological regions within the city to increase the amount of urban carbon sink; on the other hand, it is necessary to promote the low-carbon transformation of human activities at the level of social life [67] through the optimization of habitable space and the reasonable layout of public facilities, achieving carbon emissions reduction. (2) For industrial parks such as Kunshan Industrial Park, etc., where the secondary consumer (IN) was the core, the production function was its dominant function. Under the premise of ensuring sufficient employment opportunities and meeting production targets, measures such as optimizing the energy structure, introducing clean renewable energy, and improving resource utilization efficiency should be taken to reduce the energy consumption, carbon emissions density, and the negative carbon flow in the network. (3) For peri-urban and rural areas, e.g., Changshu and Zhangjiagang in the northern part of the study area, where CU and RU as primary consumers were the core, the production function was its dominant function, so the main way to achieve carbon emissions reduction is to improve agricultural production efficiency and land intensification. Measures such as promoting large-scale operations, introducing the smart agricultural production model [68], and promoting high-efficiency, low carbon cycle agriculture, should be taken. With the improvement of production efficiency, the conflict resulting from the supplementing of farmland at the cost of occupying the ecological space will be eased, thus promoting the mutual symbiosis between CU and ecological space to achieve carbon balance in this region. (4) For ecological reserves, e.g., Taihu Lake, Yangcheng Lake, etc., where producers (F, G, and W) were the core, the ecological function was its dominant function. Priority should be given to maintaining the scale of ecological space, improving the regional habitat quality, and strengthening ecological restoration and treatment to ensure the stability of the natural carbon sink function. In addition, it is necessary to control the disorderly spread of urban land and strictly adhere to the ecological protection red line, so that the pull weight of consumers in this region could remain low, thus realizing a healthy pyramid structure at the pull end of the carbon metabolic network, and promoting a positive cycle of carbon metabolism.

Limitations and Future Directions of This Study
This present study was limited by data acquisition and quantitative assessment methods, and there was room for further exploration in some aspects. This study classified the land use types in the network into eight categories based on PLE dominant functions, and the different classification ways had a great impact on the complexity and practicality of the network structure; thus, we should try to classify them in a more precise way, and determine the optimal classification scheme for different characteristics of study areas in the future. Moreover, on this basis, it would be interesting to apply the method used in this study to compare several cities or regions with different characteristics as a way to test the validity of the model. In addition, this study mainly analyzed the changes in urban carbon metabolism from the perspective of land use, and lacked the in-depth investigation of the driving mechanism of the carbon metabolism network changes. The influence of more diverse factors, such as the scale and quality of urbanization [69,70], and the recessive functions of land, on the urban carbon metabolic network should be further considered in future studies.

Conclusions
This study constructed a carbon metabolic network model in urban PLE space, accounted for the carbon emissions and sequestration of each node in the network, analyzed the changes in horizontal carbon flow among nodes in the network, in combination with the land use transfer matrix, and evaluated the comprehensive effects of the PLE space changes on the network by using the ENA method.
The results showed that the total carbon sequestration in Suzhou from 2000 to 2018 showed a fluctuating change of rising and then declining, while the total carbon emissions grew rapidly. CU and IN were the key nodes for the generation of negative carbon flow in the network, and the spatial distribution of negative carbon flow showed a discrete-concentrated-discrete change trend. The generation of positive carbon flow mainly came from the mutual transformation between production and living space, and Gusu was the main region for its generation. The exploitation relationship was the dominant ecological relationship in the network, which mainly occurred between production and ecological space, and gradually shifted from Gusu and Kunshan to Wujiang, etc., in spatial distribution. The mutualism relationship was abundant in 2005-2010 and gradually decreased in 2010-2018. The frequency of competition relationship appeared gradually increased from 17.86% in 2000-2005 to 32.14% in 2010-2018. The mutualism index MI showed a trend of first increasing and then decreasing, while the synergy index SI continued to decrease, and the negative effect brought by horizontal carbon flow in the network was relatively significant. The ecological network hierarchy at the driving end evolved from an irregular shape dominated by primary consumers (especially CU) in 2000-2005, to a pyramidal-like shape dominated by producers (F, G, and W); at the pull end, although it presented an elliptical structure dominated by secondary consumers (IN), its pull weight kept shrinking from 61.91% in 2000-2005 to 48.21% in 2010-2018, and the weight of producers (F, G, and W) rose from 1.72% to 24.33%.
The contribution of this study was to provide a new idea on how to identify the nodes in an urban carbon metabolic network. The classification of network nodes based on PLE space can better measure the impact of socio-economic development on the urban carbon cycle. The urban carbon metabolic network model constructed in this study was able to quantify and spatialize the horizontal carbon flow due to land use transition, providing an example and methodology guidance for further investigation of the land-carbon relationship. Through the ENA method, we evaluated the integrated effect of PLE space changes on the urban carbon metabolic network and proposed more targeted low-carbon land management policies. The results of this study can provide a scientific reference for optimizing the spatial layout of Suzhou City to achieve low-carbon urban land use on the one hand, and have potential applications as guidance for spatial analysis of carbon transfer at other scales and regions, as well as for low-carbon adjustment of land use patterns, on the other hand.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/land11091445/s1, Table S1: Carbon emissions coefficients of each item of carbon sources; Table S2: Carbon sequestration coefficients of each secondary land type; Table S3: Direct carbon flow matrix F from 2000 to 2005/(t C/a); Table S4: Direct carbon flow matrix F from 2005 to 2010/(t C/a); Table S5: Direct carbon flow matrix F from 2010 to 2018/(t C/a); Table S6: Ecological relationships between PLE space types in the network from 2000 to 2005; Table S7: Ecological relationships between PLE space types in the network from 2005 to 2010; Table S8: Ecological relationships between PLE space types in the network from 2010 to 2018; Figure S1: Diagrammatic sketch of SDE.