Decoupling the Relationships between Carbon Footprint and Economic Growth within an Urban Agglomeration—A Case Study of the Yangtze River Delta in China

: Carbon footprint is emerging as an effective tool for carbon emission management, especially that from fossil energy consumption. In addition, decoupling analysis is important to keep a high pace of economic growth while reducing carbon emission and its carbon footprint. Taking the Yangtze River Delta (YRD) urban agglomeration in China as a case, this paper examined the changes in carbon footprint and carbon footprint pressure by incorporating land resource limits. On this basis, we further analyzed the decoupling relationships between carbon footprint, carbon footprint pressure and economic growth. The GeoDetector was also employed to detect the spatial heterogeneity of the carbon footprint pressure. The results showed that despite the decrease of carbon emissions from 2011 to 2019 in the YRD, carbon footprint pressure still revealed an increased trend in this period. As to the decoupling relationships between carbon footprint, carbon footprint pressure and economic growth, they were improved in most of the cities in the YRD, changing from expansive coupling to weak decoupling to strong decoupling. However, the descending trend of decoupling elasticity coefﬁcient for carbon footprint pressure is smaller than that of the carbon footprint. This result could be explained by the fact that not only carbon emission but also carbon sequestration (by productive lands including forests and grasslands) pose large impacts on carbon footprint pressure. The ﬁndings indicate the necessity not only to reduce carbon emission, but also to protect productive lands to realize low carbon economy.


Introduction
The global land surface temperature increased rapidly, at a rate of 1.59 • C between 1850 and 1900 and between 2011 and 2020 [1]. The warming is likely mainly due to a continuous increase in fossil energy consumption and thus its carbon emission [2,3]. The carbon emission from fossil energy consumption is of wide concern by scientific communities in recent years. For instance, Luderer et al. (2018) explored the residual carbon dioxide (CO 2 ) emissions from fossil fuels to hold global warming well below 2 • C while pursuing efforts to limit it below 1.5 • C to reach the goals of the Paris Agreement [4]. Christophe and Paul (2015) analyzed the emission limit from fossil fuel to limit global warming to 2 • C in different regions [5]. Taking the U.S. as an example, Deutch (2017) argued that it is misleading to avoid the risks of climate change by only reducing energy consumption and carbon intensity [6].
In recent years, carbon footprint became one of the widely recognized methods to evaluate the impact of carbon emission on environment pressure. For instance, Simion et al. (2013) proposed an ecological footprint indicator by integrating land occupation, CO 2 emissions from fossil energy and nuclear energy use to perform environmental sustainability assessment in European countries [7]. Zhao  carbon footprint to assess the impact of human activities on urban environment pressure in Nanjing City, China during 2000-2009 [8]. Wiśniewski and Kistowski (2017) assessed the role and importance of carbon footprint as a tool in local planning of a low carbon economy at local levels in Poland [9]. Ma et al. (2018) analyzed the carbon footprint in passenger transport in China over the period 2006-2015 [10]. Carbon footprint (CF) frequently reflects the measure of carbon dioxide emissions directly or indirectly caused by an activity or accumulated over the life stages of a product [11,12]. There are many ways to measure carbon emission: CO 2 physical emissions, CO 2 equivalents (CO 2 eq), or by translating them into biologically productive areas indirectly [13][14][15]. In comparison with direct carbon emission, indirect carbon footprint measurement tracks energy or resource throughput of economy development and translates it into biologically productive areas necessary to produce these flows [13,16,17]. In this paper, carbon footprint refers to the carbon emission measurement using productive lands indirectly. Concretely, the carbon footprint represents the productive land area to absorb the carbon emission from fossil energy consumption. Thus, as for different cities, natural resource endowments such as productive lands have different effects on carbon footprint [17].
To promote sustainability development, many studies explored the changes in carbon emission and their decoupling relationships with economic growth. For instance, Mikayilov et al. (2018) analyzed the decoupling relationships between CO 2 emissions and Gross Domestic Product (GDP) for a group of European economies [18]. Karakaya et al. (2019) analyzed the CO 2 emission trends and the decoupling performance between CO2 emissions and economic growth in Turkey [19]. Engo (2021) explored the decoupling indicators for carbon emissions in Egypt, Morocco, Algeria and Tunisia [20]. However, few works addressed carbon footprint, carbon footprint pressure and their decoupling relationships with economic growth by incorporating natural resource limits (e.g., productive lands) [15,21].
The Yangtze River Delta (YRD) is one of the most developed and densely populated regions in China. It is also one of the largest urban agglomerations worldwide. Taking the YRD as an example, this paper examined spatiotemporal changes in the carbon footprint and their carbon footprint pressure from fossil energy consumption. The decoupling effects between carbon footprint pressure and economic growth were analyzed at both urban and regional scales. Geodetector has many advantages for measuring spatial stratified heterogeneity and exploring their determinants [22,23]. In this study, the drivers of the change of carbon footprint pressure were analyzed using the Geodetector. The objectives of this study were to understand the decoupling relationships among economic growth, carbon footprint and carbon footprint pressure by accounting for land resource limits to realize low carbon development.

Study Area
In this paper, the study area is located in the Yangtze River Delta (YRD), China ( Figure 1). The YRD lies between 118-123 • E and 28-34 • N, with an area of 104.985 km 2 . Furthermore, it encompasses the entire Shanghai City, the southern part of Jiangsu Province (i.e., Suzhou, Wuxi, Changzhou, Nanjing, Zhenjiang, Nantong, Taizhou, Yangzhou) and the northern part of Zhejiang Province (i.e., Hangzhou, Ningbo, Shaoxing, Jiaxing, Huzhou, Taizhou and Zhoushan). The major vegetation type includes subtropical evergreen broadleaf forest, cropland, urban land, grasslands, and so on. Since economic reform in China in 1978, the YRD has witnessed fast industrialization and unprecedented urbanization. The gross domestic product (GDP) in the YRD reached US $1.74 trillion in 2016 (approximately 16.2% of the total GDP in China). The overall pattern of energy consumption in the YRD was dominated by raw coal and electricity, occupying about 40-45% and 22-33% of the total energy consumption, respectively [24].

Data Sources and Preprocessing
The data in this study includes the dataset on fossil energy consumption by cities, a land use/cover dataset, a terrestrial net primary productivity dataset (NPP) and statistical datasets ( Table 1). The fossil energy consumption dataset, which spans 2001-2019, was collected from statistical yearbooks for the 16 cities separately (e.g., Shanghai statistical yearbook, Suzhou statistical yearbook, Nanjing statistical yearbook, Hangzhou statistical yearbook, and so on). In addition, other statistical datasets, including the proportion of urban population, permanent population, urban population, GDP and built-up area, were collected by cities from several statistical yearbooks, including the Shanghai statistical yearbook, the Jiangsu statistical yearbook and the Zhejiang statistical yearbook (Table 1). Regarding the land use/cover data and the NPP data, the Moderate Resolution Imaging Spectroradiometer (MODIS) products (i.e., MCD 12Q1 and MOD17A3) were downloaded from NASA's Land Processes Distributed Active Archive Center (LP DAAC). The original data were obtained from their original sinusoidal projection to geographic grid cells using the MODIS Reprojection Tool (MRT).

Data Sources and Preprocessing
The data in this study includes the dataset on fossil energy consumption by cities, a land use/cover dataset, a terrestrial net primary productivity dataset (NPP) and statistical datasets ( Table 1). The fossil energy consumption dataset, which spans 2001-2019, was collected from statistical yearbooks for the 16 cities separately (e.g., Shanghai statistical yearbook, Suzhou statistical yearbook, Nanjing statistical yearbook, Hangzhou statistical yearbook, and so on). In addition, other statistical datasets, including the proportion of urban population, permanent population, urban population, GDP and built-up area, were collected by cities from several statistical yearbooks, including the Shanghai statistical yearbook, the Jiangsu statistical yearbook and the Zhejiang statistical yearbook (Table 1). Regarding the land use/cover data and the NPP data, the Moderate Resolution Imaging Spectroradiometer (MODIS) products (i.e., MCD 12Q1 and MOD17A3) were downloaded from NASA's Land Processes Distributed Active Archive Center (LP DAAC). The original data were obtained from their original sinusoidal projection to geographic grid cells using the MODIS Reprojection Tool (MRT).

IPCC Carbon Inventory Method to Estimate Carbon Emission from Fossil Energy Consumption
By using the data on fossil energy consumption from statistical yearbooks, this paper estimated the carbon emission from 2001 to 2019 in the YRD. The fossil energy consumption in the YRD was mainly raw coal, gasoline/gasoline products, coke, and so on [24]. For assuring the temporal continuity and consistency of the energy data for the 16 cities, raw coal, gasoline and diesel were selected for analyzing the carbon emission from fossil energy consumption. There are many ways to estimate carbon emissions, including life cycle assessment (LCA), the input output method (IO), the Kaya carbon emission identity method and the IPCC carbon inventory method [28][29][30][31]. In this paper, carbon emission from fossil energy consumption was estimated using the carbon inventory method recommended by IPCC. For the i-th energy, its carbon emission (E i , ton) is calculated as: where M i is the total amount of i-th type of energy consumption (10 4 ton of standard coal equivalent, tce); and F i is the carbon emission coefficient of the i-th energy (t C tce-1). Table 2 shows the standard coal coefficient and carbon emission coefficient of raw coal, gasoline and diesel.

Estimate the Carbon Footprint of Fossil Energy Consumption and Carbon Footprint Pressure
The terrestrial vegetation in the YRD is mainly forest, croplands, urban lands, grasslands, and so on. In particular, the carbon that is absorbed into crops is frequently released back into the atmosphere through food consumption. This implies that croplands do not contribute to a net carbon sequestration [32][33][34]. On the contrary, forests and grasslands show a mainly carbon sequestration effect while affecting carbon footprint changes. In this paper, a carbon footprint model was employed by incorporating the changes in natural resource endowment such as productive lands to address the spatial changes of land use/cover. The carbon footprint of fossil energy consumption was calculated as follows: where CF is the total carbon footprint of fossil energy consumption (hm 2 ) of the i-th region; C i is fossil energy consumption (1 × 10 4 tce) of the i-th energy type: raw coal, gasoline or diesel; F i is the carbon emission coefficient of the i-th energy type (tC tce −1 ); EP f and EP g are the net ecosystem production (NEP) of the forest and grasslands, respectively. In this study, EP f and EP g were determined as the global average NEP of the corresponding vegetation type (i.e., 3.809592 t C per hm 2 for forest and 0.948229 t C per hm 2 for grasslands) [35,36]; P erf and P erg are the carbon absorption ratios of forest and grasslands in the YRD respectively. The ratios of forest and grasslands were calculated to estimate the productive land as follows: where P erf and P erg are the forest and grasslands carbon absorption ratios, respectively, A f and A g are the total area of the forest and grasslands, respectively. The carbon footprint pressure index (CFP) was calculated based on estimated carbon footprint (CF) and existing productive lands as follows: where EF f and EF g are the equilibrium factors for forest and grasslands, respectively. The global equilibrium factor is used to facilitate the addition of the actual area of the forest and the grasslands. Concretely, the equilibrium factors of the forest (EF f ) and the grasslands (EF g ) are set to 1.34 and 0.49, respectively.

Tapio Decoupling Model to Analyze the Coupling or Decoupling Effects between Carbon Footprint, Carbon Footprint Pressure and Economic Growth
To decouple environmental pressure from economic growth, the Organisation for Economic Cooperation and Development (OECD) defined the term "decoupling" as breaking the link between "environmental bads" and "economic goods" [37]. That is, the decoupling effects of carbon emission with economic growth were assessed using the Tapio decoupling model [38,39]. Concretely, the elasticity coefficient (ε) is defined as the ratio of the change rate of carbon emission (∆CO 2 ) to the change rate of GDP (∆GDP) as follows: In this paper, environmental pressure (i.e., CO 2 emissions) (∆CO 2 ) was reflected by the carbon footprint and carbon footprint pressure. In addition, economic growth was represented by using per Capita GDP (∆GDP) as an indicator. In order not to over interpret slight changes as significant, Tapio (2005) put forward the Tapio decoupling model, which gives an improvement that ±20% variation of the elasticity values around 1.0 are regarded as coupling [40]. On this basis, decoupling states in the Tapio model can be further divided into three subcategories: negative decoupling, decoupling and coupling (Table 3). Concretely, decoupling occurred when the growth rate of environmental pressure is less than that of economic driver over a specific period. In particular, decoupling is strong when economic growth increases and environmental stress decreases (elasticity <0). On the contrary, decoupling is weak when both economic growth and environmental stress increase, but the growth rate of environmental stress is less than that of economic growth (0 < elasticity < 0.8). In addition, it is recessive when both economic growth and environmental stress decrease, but the decrease rate of environmental stress is bigger than that of economic growth (elasticity >1.2). Expansive negative decoupling occurs when both economic growth and environmental stress increase, but the increase rate of environmental stress is bigger than that of economic growth (elasticity >1.2). In strong negative decoupling, economic growth decreases while environmental stress increases (elasticity <0). When both economic growth and environmental stress decrease (0 < elasticity < 0.8) and the decreased rate of environmental stress is smaller than that of economic growth, there is weak negative decoupling [38].

Driving Forces of Carbon Footprint Pressure Based on the GeoDetector Method
Many methodologies were frequently employed to explore the drivers of carbon emission, including Kaya identity, structural decomposition analysis (SDA) and logarithmic mean Divisia index (LMDI) decompose model [41][42][43][44]. The GeoDetector is a statistical method to detect spatial stratified heterogeneity and elucidate the driving factors behind it. The method has many advantages: (1) it works without the assumption of linearity of the association; and (2) it has a straight physical meaning [22,23]. Four detectors, including risk detector, factor detector, ecological detector and interaction detector, are defined in the GeoDetector [22]. The risk detector compares the difference of average values between sub-regions. The factor detector compares the accumulated dispersion variance of each sub-region with the dispersion variance of the entire study region. The smaller the ratio, the stronger the contribution of the stratum. The ecological detector compares the variance calculated from each sub-region divided according to one determinant with that divided according to another determinant. The interaction detector compares the sum of the contribution of two individual attributes vs. the contribution of the two attributes when taken together.
The spatial stratified heterogeneity is adopted by the q-statistic in the GeoDetector. Assuming that a study area is composed of N units and is stratified into h = 1, 2, . . . , L stratum. It is composed of N h units in the stratum h. The q-statistic is defined as: where σ 2 h and σ 2 are the stratum variance of effect Y or determinant D for the layer h and the whole region, respectively. The value range of q is [0, 1]. When applying different driving forces of carbon footprint pressure, the explanatory power of the determinant is stronger when the q value is larger. In this paper, five drivers were examined, including the rate of gross domestic product (X 1 ), the rate of built-up area (X 2 ), the urbanization rate (X 3 ), the rate of energy consumption per unit of GDP (X 4 ) and the terrestrial NPP (X 5 ), were examined to analyze the driver forces of carbon footprint pressure from the perspective of economic development, social development, technological progress and terrestrial productivity, respectively.
In addition, interaction detection is also used to evaluate whether the interaction of different drivers can enhance or weaken their explanatory power to carbon footprint pressure. Concretely, the GeoDetector firstly calculates the q value of two judgemental factors X 1 and X 2 , and marks them as q(X 1 ) and q(X 2 ). The combined q-statistic of the two factors are then calculated and marked as q(X 1 ∩ X 2 ). The relationship between the two factors is shown in Table 4. The details of the GeoDetector can be found in the studies by Wang et al. (2010) and Wang et al. (2016) [22,23]. Table 4. The interactive relationships in the GeoDetector.

Description
Interaction cretely, both total carbon footprint and carbon footprint pressure in the YRD revealed increasing trends before 2011 and then peaked in 2013. That is, the carbon footprint reached 5.25 × 10 7 hm 2 in 2013, which is 2.46 times that in the year 2001. Thus, despite the decreased carbon emission in the YRD from 2011 to 2019, both carbon footprint and carbon footprint pressure were still enhanced simultaneously. The phenomena could be associated with the decreased productive lands over the period of 2011-2019 in the YRD (Figure 2a).

Changes of Productive Lands, Carbon Emission and Carbon Footprint Pressure at City Scale
To further elucidate the change of carbon footprint pressure, spatial heterogeneities of productive lands, carbon emission from fossil energy consumption and its carbon footprint pressure were analyzed by cities, respectively. As shown in Figure 3a, total productive lands show large differences in different cities in the YRD. For example, there is a large amount of productive lands (mainly forests) in the southern YRD (e.g., northern Zhejiang Province). This is especially the case in Hangzhou, Ningbo, Shaoxing and Taizhou cities. In contrary, it has a small amount of productive lands in the northern YRD (e.g., southern Jiangsu Province). The phenomenon can be associated with the wide distribution of forest lands in the northern Zhejiang Province and croplands in the southern Jiangsu Province (Figure 1). Regarding the average carbon emission from fossil energy consumption during 2001-2019, it is obviously large in several cities including Shanghai, Suzhou and Ningbo (Figure 3b). For average carbon footprint pressure, it is relatively small in the northern Zhejiang Province in the southern YRD (Figure 3c). The result could be associated with the wide distribution of forest lands there. In comparison, carbon footprint pressure is more severe in the northern YRD. This is especially the case in Shanghai City. We found that both the averages and the slopes of the carbon footprint pressure in Shanghai City are particularly bigger due to its excessively high carbon emissions in past decades. In addition, although the average carbon footprint pressure in Suzhou and Wuxi is not very large as in Shanghai City, the slope of the corresponding carbon footprint pressure reveals a growing trend over the two cities in comparison to most other cities over the period of 2001-2019 (Figure 3d). This phenomenon could be mainly caused by the sharp decrease of productive lands in the two cities. That is, productive lands witnessed a decrease, at rates of 4153 and 4040 hm 2 a −1 from 2001 to 2019 in Wuxi and Suzhou, respectively. The cumulative decreases in the past 19 years accounted for 33.64% and 26.86% of the total productive lands in the two cities in 2001. This result indicates the necessity to protect productive lands (e.g., forests and grasslands) to realize low carbon economy. sure reveals a growing trend over the two cities in comparison to most other cities over the period of 2001-2019 (Figure 3d). This phenomenon could be mainly caused by the sharp decrease of productive lands in the two cities. That is, productive lands witnessed a decrease, at rates of 4153 and 4040 hm 2 a −1 from 2001 to 2019 in Wuxi and Suzhou, respectively. The cumulative decreases in the past 19 years accounted for 33.64% and 26.86% of the total productive lands in the two cities in 2001. This result indicates the necessity to protect productive lands (e.g., forests and grasslands) to realize low carbon economy.

Decoupling analysis on the Relationships between Carbon Footprint, Carbon Footprint Pressure and Economic Growth
The decoupling relationships between carbon footprint and economic growth were analyzed by comparing their relative changes as well as the elasticity change of the carbon footprint and the carbon footprint pressure in the YRD from 2001 to 2018. As shown in Figure 4, the change rate of GDP in the YRD was all positive from 2001 to 2018. This indicates the continuous economic growth in the past 19 years. In addition, both the carbon footprint and the carbon footprint pressure revealed growing trends over the period of 2001-2019. The elasticity coefficients were therefore decreasing from 2001 to 2019. Consequently, the decoupling relationships between the carbon footprint and economic growth mainly changed from expansive coupling to weak decoupling to strong decoupling. Similar changes were also found for the carbon footprint pressure. Despite it, the elasticity coefficient of the carbon footprint pressure is bigger than that of the carbon footprint before 2009. However, the relationship was reversed after 2009. That is, the descending trend of the carbon footprint pressure is smaller than that of the carbon footprint. This difference could be associated with the changes in productive lands in the YRD.

Decoupling analysis on the Relationships between Carbon Footprint, Carbon Footprint Pressure and Economic Growth
The decoupling relationships between carbon footprint and economic growth were analyzed by comparing their relative changes as well as the elasticity change of the carbon footprint and the carbon footprint pressure in the YRD from 2001 to 2018. As shown in Figure 4, the change rate of GDP in the YRD was all positive from 2001 to 2018. This indicates the continuous economic growth in the past 19 years. In addition, both the carbon footprint and the carbon footprint pressure revealed growing trends over the period of 2001-2019. The elasticity coefficients were therefore decreasing from 2001 to 2019. Consequently, the decoupling relationships between the carbon footprint and economic growth mainly changed from expansive coupling to weak decoupling to strong decoupling. Similar changes were also found for the carbon footprint pressure. Despite it, the elasticity coefficient of the carbon footprint pressure is bigger than that of the carbon footprint before 2009. However, the relationship was reversed after 2009. That is, the descending trend of the carbon footprint pressure is smaller than that of the carbon footprint. This difference could be associated with the changes in productive lands in the YRD. 1, 10, x FOR PEER REVIEW 10 of 15

City-Based Decoupling Analysis on Carbon Footprint Pressure and Economic Growth
We further analyzed the decoupling relationships between carbon footprint pressure and economic growth in the YRD by cities. It is observed that the decoupling relationships were improved in most of the cities in the YRD in the past 19 years. For instance, Figure  5a shows the decoupling characteristics according to the average elasticity coefficient from 2001 to 2009 for different cities. We found that 12 of 16 cities were in a weak decoupling stage over this period. Meanwhile, 4 cities, including Taizhou in Jiangsu Province, Suzhou, Jiaxing and Taizhou in Zhejiang Province, were in an expansive coupling stage. However, the decoupling relationships were then improved in most of the cities over the period of 2010-2018 ( Figure 5b) as: (1) from expansive coupling to weak decoupling stage (e.g., in Taizhou in Jiangsu Province, Suzhou, Jiaxing and Taizhou in Zhejiang Province); (2) from weak coupling to strong decoupling stage (e.g., in Shanghai, Wuxi, Changzhou, Yangzhou, Hangzhou, Shaoxing and Ningbo). Despite these improvements, the decoupling relationships could also be degraded in the YRD. For instance, the decoupling relationship in Zhoushan experienced a conversion from the weak decoupling stage to the expansive coupling stage.

City-Based Decoupling Analysis on Carbon Footprint Pressure and Economic Growth
We further analyzed the decoupling relationships between carbon footprint pressure and economic growth in the YRD by cities. It is observed that the decoupling relationships were improved in most of the cities in the YRD in the past 19 years. For instance, Figure 5a shows the decoupling characteristics according to the average elasticity coefficient from 2001 to 2009 for different cities. We found that 12 of 16 cities were in a weak decoupling stage over this period. Meanwhile, 4 cities, including Taizhou in Jiangsu Province, Suzhou, Jiaxing and Taizhou in Zhejiang Province, were in an expansive coupling stage. However, the decoupling relationships were then improved in most of the cities over the period of 2010-2018 ( Figure 5b) as: (1) from expansive coupling to weak decoupling stage (e.g., in Taizhou in Jiangsu Province, Suzhou, Jiaxing and Taizhou in Zhejiang Province); (2) from weak coupling to strong decoupling stage (e.g., in Shanghai, Wuxi, Changzhou, Yangzhou, Hangzhou, Shaoxing and Ningbo). Despite these improvements, the decoupling relationships could also be degraded in the YRD. For instance, the decoupling relationship in Zhoushan experienced a conversion from the weak decoupling stage to the expansive coupling stage.

Driving Factors of Carbon Footprint Pressure Changes
In this paper, we also examined the impacts of natural, environmental and socioeconomic factors on the changes of carbon footprint pressure in the YRD. According to the results of the factor detector in the GeoDectector, these factors were ranked according to their effects on the carbon footprint pressure in the following order: the rate of built-up area (0.290) > terrestrial NPP (0.241) > urbanization rate (0.135) > energy consumption per GDP (0.134) > GDP rate (0.133) ( Table 5).

Driving Factors of Carbon Footprint Pressure Changes
In this paper, we also examined the impacts of natural, environmental and socioeconomic factors on the changes of carbon footprint pressure in the YRD. According to the results of the factor detector in the GeoDectector, these factors were ranked according to their effects on the carbon footprint pressure in the following order: the rate of built-up area (0.290) > terrestrial NPP (0.241) > urbanization rate (0.135) > energy consumption per GDP (0.134) > GDP rate (0.133) ( Table 5).  Table 5 also shows the response of the carbon footprint pressure to the interaction between two of the driving factors according to the interaction detector in the GeoDetector. According to the interactive detection analysis, all the interactive q values appeared to be higher than any q value of a sole factor. It is obvious that the interaction based on the urbanization rate is bigger than any other factors. In particular, the combined q value of the rate of built-up area and urbanization rate was 0.975, that of energy consumption per GDP and urbanization rate was 0.971, that of GDP rate and urbanization rate was 0.933, and that of terrestrial NPP and urbanization rate was 0.477. The combinations of the above-mentioned factors can better explain spatial variability of the changes in the carbon footprint pressure than most of the others. The results indicate that the urbanization rate   Table 5 also shows the response of the carbon footprint pressure to the interaction between two of the driving factors according to the interaction detector in the GeoDetector. According to the interactive detection analysis, all the interactive q values appeared to be higher than any q value of a sole factor. It is obvious that the interaction based on the urbanization rate is bigger than any other factors. In particular, the combined q value of the rate of built-up area and urbanization rate was 0.975, that of energy consumption per GDP and urbanization rate was 0.971, that of GDP rate and urbanization rate was 0.933, and that of terrestrial NPP and urbanization rate was 0.477. The combinations of the above-mentioned factors can better explain spatial variability of the changes in the carbon footprint pressure than most of the others. The results indicate that the urbanization rate played critical roles on the carbon footprint pressure in the YRD. In addition, the interactions between terrestrial NPP and other drivers are relatively strong. For instance, the q value of the rate of built-up area and terrestrial NPP was 0.546, that of energy consumption per GDP and the urbanization rate was 0.483, and that of the GDP rate and the urbanization rate was 0.482. This indicates that, besides the urbanization rate, carbon absorption that related to natural resource endowment (e.g., productive lands) can also be important to ease the carbon footprint pressure.

Discussions
In general, absolute decoupling means that economic growth coincides with absolute reduction in emission or resource use. On the contrary, relative decoupling denotes that resource use or emission increase less so than does the GDP. Actually, relative decoupling is frequent for most of the cases [45]. Regarding the carbon footprint and carbon footprint pressure, the decoupling relationships were also obvious for the cities in the YRD. We found that the elasticity values revealed decreased trends in all the cities in the YRD. However, the carbon footprint pressure revealed increased positive trends in most of the cities in the YRD from 2001 to 2019, except for Hangzhou and Nantong. Thus, it mainly revealed a relative decoupling for most of the cities in the YRD. This is consistent with the results of Haberl et al. (2020) [45]. Regarding Hangzhou City and Nantong City, carbon footprint pressure showed a continuously decreased trend after 2007 and 2010, respectively. Despite a little fluctuation, absolute decoupling existed in the past decade. This could be associated with the increased productive lands in the two cities. Thus, not only by reducing carbon emission, it is also critical to increase productive lands and thus the uptake of CO 2 by the terrestrial biosphere to ease the environmental pressure [6].
The decoupling state is also found to be associated with the stages of economic growth in different cities. Pilatowska and Wlodarczyk (2018) analyzed the decoupling relationships between economic growth and CO 2 emissions in the European Union (EU) countries [46]. The EKC hypothesis, which assumes that environmental degradation increases with per capita income during the early stages of economic growth, and then declines with per capita income after passing beyond an income turning point [47], was tested by dividing the EU countries into three groups depending on the level of knowledgebased economy. They found that the EKC hypothesis is valid for the most high-level and some middle level knowledge-based economies [46]. Moreover, these countries are characterized by either a decreasing trend in energy intensity or the diversified energy consumption mix with a significant share of nuclear and renewable energy and also a smaller share of solid fuels. This phenomenon is consistent with the results of this study on carbon footprint pressure. For example, four cities in the early stages of rapid economic growth or low-level development, including Taizhou in Jiangsu Province, and Suzhou, Jiaxing and Taizhou in Zhejiang Province, were in an expansive coupling stage from 2001 to 2009. In addition, other cities were in a weak decoupling stage over this period. However, the decoupling relationships between carbon footprint pressure and economic growth were then dramatically improved, especially in the cities with rapid economic growth and increased productive lands (e.g., Shanghai, Wuxi, Hangzhou, Ningbo, and so on). This indicates the necessity not only to reduce carbon emission, but also to protect productive lands to realize low carbon economy.
In addition, Li et al. (2021) classified the cities in the Yangtze River Economic Belt according to the differences in the decoupling stage; they identified three categories: pre-decoupling, post-decoupling, and unstable decoupling [48]. They inferred that predecoupling cities are mainly dominated by high-tech industries or service industries, such as Shanghai City. Our findings partly agree with the results from Li et al. (2021). For instance, most of the cities in the YRD were in a weak decoupling stage between carbon footprint pressure and economic growth over the period 2001-2009. After 2009, Shanghai, Wuxi, Yangzhou, Hangzhou, Shaoxing and Ningbo experienced a strong decoupling stage. The phenomena could be associated with decreased carbon emission, as well as the increase of productive lands (especially in the northern Zhejiang Province) in these cities. By comparison, post-decoupling cities are characterized by resources and geographical advantages, and some cities have accelerated the decoupling by eliminating high-emission industries, extending industrial chains, and nurturing new pillar industries (e.g., the coal chemical industry) [48]. In the study area, similar results were found as well. For example, several cities, including Taizhou in Jiangsu Province, Suzhou, Jiaxing and Taizhou in Zhejiang Province, changed from expansive decoupling during 2001-2009 to weak decoupling over the period 2010-2018. In addition, unstable decoupling cities were in a critical period of industrialization with high-emissions and low-value-added industries [48]. This was especially the case for Zhoushan City in the YRD. We found that the carbon emission in this city did not change from 2001 to 2010. However, its total carbon emission in 2019 reached 4.2 times that in 2010.
Decoupling relationships between economic growth and carbon emission were frequently attributed to per capita GDP, urbanization rate, industrial structure, energy intensity, etc. [49]. In this study, we found that the built-up area rate and terrestrial NPP explained large amounts of the changes in the carbon footprint pressure. In addition, the joint explanatory power mostly reached more than 93% based on the interaction analysis between the urbanization rate and other drivers. The phenomena indicated that not only carbon emission, but also carbon sequestration, posed large impacts on the carbon footprint, as well as carbon footprint pressure. Further concerns should be addressed about the drivers of carbon sequestration such as changes in productive lands.

Conclusions
In past studies, few works addressed the decoupling relationship between the carbon footprint and economic growth by incorporating natural resource limits. This paper found that the decoupling relationships in the YRD were improved when incorporating the changes in productive lands in the past decades. However, the descending trend of decoupling the elasticity coefficient for carbon footprint pressure was smaller than that of the carbon footprint. In addition, not only carbon emission, but also carbon sequestration, posed large impacts on carbon footprint pressure. The findings appealed for a differential strategy to coordinate economic development and environment protection according to each city's characteristic to realize a low carbon economy. Considering the electricity transmission from one place to another, this study used only statistical data on fossil energy consumption. However, it was subject to data limitations on fossil energy consumption at the city scale. That is, for assuring the temporal continuity and consistency of the energy data at different cities, we selected only raw coal, gasoline and diesel to reflect the changes in carbon emission. More cases should be conducted using complete and comparable fossil energy data on multiple energy types. Furthermore, the Tapio model only provides the decoupling state from an elastic perspective in a certain period. However, the relationships between economic growth and carbon footprint are variable and reflects different patterns of economic development. An in-depth analysis (e.g., wavelet transform), is essential to be conducted in the time-frequency domain for richer results [50]. This paper mainly focused on the relationships between economic development and environmental factors (i.e., carbon emission and carbon footprint); it is essential to pay more attention to the role and importance of the indirect carbon footprint (e.g., biologically productive area as an indicator) as a tool in spatial planning of a low carbon economy in future studies [9].  Data Availability Statement: The datasets in this study are available within this article. The data which were derived from the original datasets but not aforementioned are available upon requests.