*1.1. Motivation*

Since the industrial civilization, under the combined influence of human activities and natural factors, the global warming trend has become increasingly significant. How to deal with the challenges posed by climate change to sustainable development has become a major scientific issue facing mankind [1]. Building a low-carbon future development mode has gradually become a global consensus. The United Nations Framework Convention on Climate Change (UNFCCC), as the world's first international convention to control carbon dioxide (CO2) emissions, provides a basic framework for international cooperation on climate change [2]. The 21st United Nations Climate Change Conference (UNFCCC COP21) held in 2015 formally adopted the Paris Agreement, which sets the global average temperature increase within 2 ◦C as an explicit goal [3]. However, with the current trend of CO<sup>2</sup> emissions, the temperature control targets of the Paris Agreement will be difficult to achieve. Deep CO<sup>2</sup> reductions in the coming years are key to achieving that goal [4,5]. Currently, 133 countries have made carbon neutral commitments. China has adopted "carbon neutrality" as a long-term national strategy to address climate change. Compared with developed countries that have achieved carbon peak, developing countries are currently facing the dual pressures of low-carbon transformation and economic development [6,7].

Water resources are the material basis for human survival and the key support for social development and ecological protection. The field of water resources is an important area for implementing the goal of "carbon neutrality" and supporting sustainable development. The "2030 Carbon Peak Action Plan" issued by the Chinese government in 2021 regards hydropower generation, water ecological protection, and efficient utilization of water resources as important ways to promote carbon neutrality [8]. In addition, improving the carbon emission accounting mechanism in different fields and carrying out research on CO<sup>2</sup> emission accounting methods are also important contents of the action plan. Therefore, it is of great significance to explore the CO<sup>2</sup> emission equivalent analysis (CEEA) method for different water resource behaviors (WRBs), and to find a reference "ruler" for the accounting of CO<sup>2</sup> emission equivalent (CEE) in the field of water resources. This "ruler" also has certain positive significance for the global control of CO<sup>2</sup> emissions.

#### *1.2. Literature Review*

The identification of source-sink relationships and the assessment of emission intensity of CO<sup>2</sup> are the basis for scientific research in the field of climate change. Accounting for carbon emission and sink effects has been a popular research topic in this field. In terms of carbon emission assessment, most studies have focused on the carbon emission intensity of human activities and carbon footprint accounting in different fields. Carbon footprint can be simply defined as the total amount of greenhouse gases (GHG), mainly carbon dioxide, released directly or indirectly from human activities [9]. Carbon footprint accounting can be divided into macro and micro levels. The methods involved are mainly input-output analysis (IOA) [10] and life cycle assessment (LCA) [11]. In recent years, many scholars have carried out multidimensional accounting of carbon footprints at the mesoscale and macroscale, such as countries [12], cities [13], and industries [14]. Carbon footprint accounting research on the microscale such as enterprise [15], product [16], and technology [17] is also ongoing. For example, Chai et al. [17] used a life-cycle approach to compare the carbon footprints of three mainstream wastewater treatment technologies in China. Accounting for CO<sup>2</sup> emissions caused by land use change [18] is also a popular research topic. In addition, some studies have explored the carbon emission effects of lake wetland [19], reservoir [20], farmland [21], and other ecosystems. For example, Keller et al. conducted a study on carbon emissions from reservoir fallout zones and concluded that reservoirs are a source rather than a sink of carbon in the global carbon cycle [20].

In terms of carbon sink effect assessment, relevant studies have mostly focused on carbon sink effects in terrestrial ecosystems such as forest, grassland, and wetland. The research methods include ground investigation, eddy covariance carbon flux observation [22], ecosystem process model simulation [23], etc. It is worth mentioning that in 2019, the IPCC added the "top-down" atmospheric inversion methodological system to the basic methodological framework for future global GHG accounting [24]. Atmospheric inversion methods have received more attention in recent years in the study of ecosystem carbon sinks. For example, Fernández-Martínez et al. [25] analyzed the trend of global carbon sinks based on atmospheric inversion and vegetation models and explored the relationship between CO<sup>2</sup> emissions and temperature. The research on the carbon sink effects of terrestrial ecosystems has formed a sound theoretical and methodological system. However, compared to terrestrial ecosystems, research on carbon sinks in marine ecosystems is still in the developmental stage.

In particular, studies on energy consumption and CO<sup>2</sup> emission accounting in the field of water resources have been carried out by relevant institutions. In 2005, California released a report on California's water–energy nexus [26], which systematically studied the energy consumption of water supply, water transmission, water utilization, and water

treatment in California. Further, the River Network, a U.S. research organization, released The Carbon Footprint of Water [27] in 2009 comprehensively assessed the various waterrelated carbon footprints in the United States. In addition, the issue of energy consumption and carbon emissions in the field of water resources has been actively discussed by scholars from different countries. LCA and IOA are still the most mainstream research methods under this research topic. The research results basically cover all aspects of social water cycle and urban water system. Although more research has been done in urban water systems, the scale of research has involved countries [28], regions [29], cities [30], schools [31], etc. For example, Wakeel et al. [28] analyzed the energy consumption of different countries in various segments of the social water cycle and compared different methods for measuring energy consumption in the water sector. Using energy consumption as a bridge to quantitatively assess the relationship between water and carbon emissions, Rothausen and Conway [29] systematically explored the GHG emissions in the water sector in different countries and regions. At the urban scale, Valek et al. [32] quantified the CO<sup>2</sup> emissions associated with the water system in Mexico City based on survey data. Based on the LCA method, Friedrich et al. [30] assessed the carbon footprint of different parts of the urban water system (storage, treatment, distribution, collection, and wastewater treatment) in Durban, South Africa. Similarly, Sambito and Freni [33] used the LCA method to quantify the carbon footprint of a metropolitan water system in Italy. In addition, Li et al. [31] quantified the water–energy–carbon relationship on a campus in northern China and explored the spatial distribution pattern of carbon sources/sinks at a small scale.

In addition to the overall study of the carbon emissions from the social water cycle and urban water system, some scholars have carried out targeted discussions on different links of the water system (water production and supply, desalination, water utilization, wastewater treatment, etc.). In terms of water production and supply, relevant studies mainly focus on carbon footprint accounting of water supply system and water distribution system. For example, Fang and Newell [34] used the LCA method to assess the carbon footprint of Southern California's water supply system, arguing that the carbon footprint of local reclaimed water is much lower than that of long-distance water supply. Boulos and Bros [35] proposed a WNEE (Water network energy efficiency) method for measuring the carbon footprint of energy consumption in a water distribution system, which was applied in a European city. Moreover, Heihsel and Lenzen [36] constructed a multi-regional inputoutput model (MIOA) for measuring GHG emissions from desalination in Australia, which provides a solution for the calculation of the carbon footprint of desalination at a macroscale. In terms of water end-use, the studies mainly cover energy consumption and carbon emission measurement for domestic and agricultural water use. For example, Siddiqi and Fletcher [37] summarized the range of energy intensity of domestic water and agricultural water in the end-use process. Escriva-Bou et al. [38] simulated GHG emissions associated with domestic water use in California using probability distribution models and emission factors. Wang et al. [39] evaluated the carbon footprint of agricultural groundwater use in 31 provinces of China based on statistical survey data. In terms of wastewater collection and treatment, the carbon emission effects and measurement methods of wastewater treatment plants and municipal wastewater sectors in different countries such as China [40], the United States [41], and Italy [42] have been deeply discussed. Further, the carbon emission effects of different wastewater treatment technologies and options have been studied in comparison [43]. It is worth mentioning that research on carbon emissions accounting for water saving behavior has also been carried out, covering different scales such as city [44] and campus [45]. In addition, Wang et al. [46] explored the water footprint and carbon footprint in hydropower stations in China and made recommendations for carbon emission reduction of hydropower stations. Some of the studies addressing the carbon emission effects in the water resources sector are summarized in Table 1.


**Table 1.** Selected representative literature of carbon emission effect studies in the field of water resources.

In general, relevant research results provide important reference value for the quantitative identification of water-carbon relationship and carbon neutrality in the field of water resources. However, some of the studies are too targeted, difficult to obtain data, and the experimental methods are not easily reproducible to meet the demand for systematic research on CO<sup>2</sup> emission accounting in the field of water resources. In addition, carbon dioxide emissions related to water resource behaviors involve many links and are not limited to the scope discussed above. How to make a comprehensive and feasible "ruler" to provide convenience and reference for the estimation of water-related CO<sup>2</sup> emissions is still a problem to be further explored. To facilitate the discussion of CO<sup>2</sup> emission or absorption effects in the field of water resources, this paper is devoted to the study of water resource behaviors, that is, a series of activities related to the development, allocation, utilization, and protection of water resources. Different links of the water cycle or water resources system can be understood as different WRBs. Researching the methodologies for quantifying the CO<sup>2</sup> emission effects of different WRBs is a further refinement and extension of the carbon source/sink effects accounting in the field of water resources.

### *1.3. Contribution and Objectives*

Based on the literature review, this study proposes a carbon dioxide emission equivalent analysis (CEEA) method for several common water resource behaviors (WRBs) from four dimensions: water resources development, water resources allocation, water resources utilization, and water resources protection. The function table of CO<sup>2</sup> emission equivalent analysis (FT-CEEA) of WRBs is constructed for the first time, which provides a method set for researchers in different regions and industries to evaluate the CO<sup>2</sup> emission equivalent (CEE) of WRBs. Compared to existing studies, the contributions of this study are: (1) The CEEA method is proposed to realize the quantitative calculation of CEE for different WRBs; (2) the FT-CEEA is developed to provide a convenient and feasible "ruler" for the measurement of CEE in the field of water resources; (3) based on the FT-CEEA, the spatial distribution characteristics of CO<sup>2</sup> emission or absorption effects of WRBs in 31 provinces in China are clarified.

This paper is organized as follows: Section 2 is the introduction of the CEEA method and FT-CEEA; Section 3 is the study area and data description, as well as results analysis and discussion; Section 4 is the main conclusion and research prospect.

#### **2. Methodology**

#### *2.1. CEEA Method Framework of Water Resource Behaviors*

Water resource behavior (WRB) is a collective term for a range of activities related to the development, allocation, utilization, and protection of water resources. The carbon dioxide emission equivalent (CEE) of water resource behaviors refers to the CO<sup>2</sup> emission or absorption effects directly or indirectly caused by water resource behaviors. In this study, the method to quantify the CEE generated by WRBs is called the carbon dioxide emission equivalent analysis method (CEEA) of WRBs. Most WRBs do not emit CO<sup>2</sup> themselves and are not explicitly linked to CO2. However, WRBs are often accompanied by energy consumption, which in turn leads to CO<sup>2</sup> emissions. Therefore, compared to "carbon dioxide emission", "carbon dioxide emission equivalent" is more accurate to represent the CO<sup>2</sup> emission or absorption effects of WRBs.

This study proposes the CEEA method and develops the FT-CEEA (the function table of carbon dioxide emission equivalent analysis), aiming to find a reference "ruler" to provide methodological reference and technical support for the accounting of CEE related to WRBs. The general idea of the CEEA method is to develop diversified CEE functions for WRBs in different dimensions by direct reference, refinement, and innovation, and finally integrate them into a unified calculation platform to form a relatively complete "ruler", namely FT-CEEA. The general idea diagram of the CEEA method is shown in Figure 1.

**Figure 1.** The general idea of the CEEA method. **Figure 1.** The general idea of the CEEA method.

The WRBs involve a wide range of fields and factors, and the CEE accounting of WRBs should have a clear system boundary to avoid the infinite extension of indirect calculations. The principle of this study for system boundary formulation is to focus on CEE directly caused by WRBs, with appropriate consideration of indirect CEE that are closely related to such WRBs. Based on the definition of WRBs, the system boundary of CEE accounting is determined, as shown in Figure 2. The categories of WRBs can be roughly divided into water resource development behaviors (WRDBs), water resource allocation behaviors (WRABs), water resource utilization behaviors (WRUBs), and water resource protection behaviors (WRPBs). Each category contains a variety of typical WRBs, each WRB has a corresponding CEEA method. The WRBs involve a wide range of fields and factors, and the CEE accounting of WRBs should have a clear system boundary to avoid the infinite extension of indirect calculations. The principle of this study for system boundary formulation is to focus on CEE directly caused by WRBs, with appropriate consideration of indirect CEE that are closely related to such WRBs. Based on the definition of WRBs, the system boundary of CEE accounting is determined, as shown in Figure 2. The categories of WRBs can be roughly divided into water resource development behaviors (WRDBs), water resource allocation behaviors (WRABs), water resource utilization behaviors (WRUBs), and water resource protection behaviors (WRPBs). Each category contains a variety of typical WRBs, each WRB has a corresponding CEEA method.

*Water* **2023**, *15*, 431 7 of 25

a method set for researchers in different regions and industries to evaluate the CO<sup>2</sup> emission equivalent (CEE) of WRBs. Compared to existing studies, the contributions of this study are: (1) The CEEA method is proposed to realize the quantitative calculation of CEE for different WRBs; (2) the FT-CEEA is developed to provide a convenient and feasible "ruler" for the measurement of CEE in the field of water resources; (3) based on the FT-CEEA, the spatial distribution characteristics of CO<sup>2</sup> emission or absorption effects of

This paper is organized as follows: Section 2 is the introduction of the CEEA method and FT-CEEA; Section 3 is the study area and data description, as well as results analysis

Water resource behavior (WRB) is a collective term for a range of activities related to the development, allocation, utilization, and protection of water resources. The carbon dioxide emission equivalent (CEE) of water resource behaviors refers to the CO<sup>2</sup> emission or absorption effects directly or indirectly caused by water resource behaviors. In this study, the method to quantify the CEE generated by WRBs is called the carbon dioxide emission equivalent analysis method (CEEA) of WRBs. Most WRBs do not emit CO<sup>2</sup> themselves and are not explicitly linked to CO2. However, WRBs are often accompanied by energy consumption, which in turn leads to CO<sup>2</sup> emissions. Therefore, compared to "carbon dioxide emission", "carbon dioxide emission equivalent" is more accurate to repre-

This study proposes the CEEA method and develops the FT-CEEA (the function table of carbon dioxide emission equivalent analysis), aiming to find a reference "ruler" to provide methodological reference and technical support for the accounting of CEE related to WRBs. The general idea of the CEEA method is to develop diversified CEE functions for WRBs in different dimensions by direct reference, refinement, and innovation, and finally integrate them into a unified calculation platform to form a relatively complete "ruler", namely FT-CEEA. The general idea diagram of the CEEA method is shown in Figure 1.

and discussion; Section 4 is the main conclusion and research prospect.

WRBs in 31 provinces in China are clarified.

*2.1. CEEA Method Framework of Water Resource Behaviors*

sent the CO<sup>2</sup> emission or absorption effects of WRBs.

**2. Methodology**

*2.2. CEEA Method to Water Resource Development Behaviors*

Water resource development behaviors (WRDBs) refer to a series of activities related to water resources development. In this study, WRDBs are preliminarily defined as sur-

(1) Surface water lifting (WRDB1): Surface water resources are extracted from natural rivers or lakes to higher elevations using water extraction projects to achieve centralized treatment and unified distribution of "raw freshwater". The electric energy consumed by the water lifting project is converted into the mechanical energy needed for water resources lifting, so the CO<sup>2</sup> emissions of this behavior are mainly concentrated in the energy consumption link of the water lifting project. The CEEA formula of WRDB1 based on

1

*h*

6

(1)

(2)

(WRDB3), raw water treatment (WRDB4), and seawater desalination (WRDB5).

*E Q EI EF*

*E*

1 1 1

= 

1

=

*I*

3.6 10

 *g*

 

 

**Figure 2.** Schematic diagram of CEE accounting for WRBs. **Figure 2.** Schematic diagram of CEE accounting for WRBs.

emission factor [52] is as follows:

#### *2.2. CEEA Method to Water Resource Development Behaviors*

Water resource development behaviors (WRDBs) refer to a series of activities related to water resources development. In this study, WRDBs are preliminarily defined as surface water lifting (WRDB1), groundwater extraction (WRDB2), reservoir storage (WRDB3), raw water treatment (WRDB4), and seawater desalination (WRDB5).

(1) Surface water lifting (WRDB1): Surface water resources are extracted from natural rivers or lakes to higher elevations using water extraction projects to achieve centralized treatment and unified distribution of "raw freshwater". The electric energy consumed by the water lifting project is converted into the mechanical energy needed for water resources lifting, so the CO<sup>2</sup> emissions of this behavior are mainly concentrated in the energy consumption link of the water lifting project. The CEEA formula of WRDB1 based on emission factor [52] is as follows:

$$E\_1 = Q\_1 \times EI\_1 \times EF \tag{1}$$

$$EI\_1 = \frac{\rho \times g \times h\_1}{3.6 \times 10^6 \times \eta} \tag{2}$$

$$EF = \frac{\sum\_{i} \left( FC\_{i,y} \times NCV\_{i,y} \times EF\_{\text{CO}\_2, i, y} \right)}{EG\_y} \tag{3}$$

where *E*<sup>1</sup> is the carbon dioxide emission equivalent of surface water lifting behavior, kg; *Q*<sup>1</sup> is the amount of surface water lifting, m<sup>3</sup> ; *EI*<sup>1</sup> is the energy intensity of WRDB1 (the amount of electricity required to lift per unit of surface water) kWh/m<sup>3</sup> ; *ρ* is the density of surface water (typically 1000), kg/m<sup>3</sup> ; *g* is the acceleration of gravity (typically 9.8), m/s<sup>2</sup> ; *h*<sup>1</sup> is the surface water lifting head, m; *η* is the efficiency of the water lifting project; *EI*<sup>1</sup> is the power system CO<sup>2</sup> emission factor (the amount of CO<sup>2</sup> emitted per unit of electricity consumed), kg/kWh; *EG* is the total net power generation during the calculation period of the power system, kWh; *FC* is the consumption of fuel by the generator set during the calculation period, in mass or volume units; *NCV* is the average low-level heat content of the fuel during the calculation period, in GJ/mass or volume units; *EF*CO<sup>2</sup> is the CO<sup>2</sup> emission factor of the fuel (amount of CO<sup>2</sup> emitted per unit of energy) during the calculation period, kgCO2/GJ; *i* is the type of fossil fuels consumed to generate electricity; and *y* is the year. Power-related departments in different countries will regularly release the *EF* of the power system. For example, the Ministry of Ecology and Environment of China has issued *EF* reference values for different provinces in China. In China, *EI*<sup>1</sup> is mainly related to the water head, which can be 0.2 kWh/m<sup>3</sup> [53,54] on average. The global level can refer to the value range given by relevant research: 0.0002–1.74 kWh/m<sup>3</sup> [55].

(2) Groundwater extraction (WRDB2): Similar to the principle of WRDB1, groundwater extraction behavior also needs to convert the electrical energy of pumping equipment into the mechanical energy required for groundwater rise. CO<sup>2</sup> emissions are mainly concentrated in the energy consumption of pumping equipment:

$$E\_2 = Q\_2 \times EI\_2 \times EF \tag{4}$$

$$EI\_2 = \frac{9.8 \times \rho \times h\_2}{3.6 \times 10^6 \times \eta} \tag{5}$$

where *E*<sup>2</sup> is the CEE of WRDB2, kg; *Q*<sup>2</sup> is the amount of groundwater extraction, m<sup>3</sup> ; *EI*<sup>2</sup> is the energy intensity of WRDB2 (the amount of electricity required to extract per unit of groundwater) kWh/m<sup>3</sup> ; *h*<sup>2</sup> is the groundwater depth, m. Other variables have the same meaning as above. Unlike surface water, the energy intensity of WRDB2 varies considerably with different groundwater burial depths. The value of *EI*<sup>2</sup> can be obtained according to the actual situation of the study area, and *EI*<sup>2</sup> in different regions of China can also refer to Table 4 [39]. In addition, the *EI*<sup>2</sup> of different countries is available in studies: 0.18–0.49 kWh/m<sup>3</sup> (USA) [56], 0.48–0.53 kWh/m<sup>3</sup> (Australia) [57], and 0.37–1.44 kWh/m<sup>3</sup> (Global) [55].

(3) Reservoir storage (WRDB3): The energy consumption of WRDB3 mainly comes from the daily operation and management of water storage infrastructure [54], such as gate control, lighting, and monitoring equipment operation. This process also produces CO<sup>2</sup> emissions. The calculation formula is as follows:

$$E\_3 = Q\_3 \times EI\_3 \times EF \tag{6}$$

where *E*<sup>3</sup> is the CEE of WRDB3, kg; *Q*<sup>3</sup> is the actual volume of water stored in the reservoir, m<sup>3</sup> ; *EI*<sup>3</sup> is the energy intensity of WRDB3 (the amount of electricity required to store per unit of water in the reservoir) kWh/m<sup>3</sup> . *EI*<sup>3</sup> varies due to differences in reservoir conditions in different regions. Field visits to reservoirs can be conducted to obtain the value of *EI*3. *EI*<sup>3</sup> can also refer to existing research. Studies have shown that the energy intensity range of WRDB3 in China is [0.07,0.2] kWh/m<sup>3</sup> [54], and 0.14 kWh/m<sup>3</sup> can be used to study the average state of China [58].

(4) Raw water treatment (WRDB4): After taking raw water from the water source, it needs to be treated by the waterworks, including coagulation, sedimentation, filtration, and disinfection [41]. Each process relies mainly on electricity to maintain the normal operation of the processing equipment, so WRDB4 also produces CO<sup>2</sup> emissions [51]:

$$E\_4 = Q\_4 \times EI\_4 \times EF \tag{7}$$

where *E*<sup>4</sup> is the CEE of WRDB4, kg; *Q*<sup>4</sup> is the volume of raw water treatment, m<sup>3</sup> ; *EI*<sup>4</sup> is the energy intensity of WRDB4 (the amount of electricity required to treat per unit of raw water) kWh/m<sup>3</sup> . *EI*<sup>4</sup> can be determined by the statistical calculation of energy consumption data of each link of WRDB4. According to the yearbook of Chinese urban water supply, the national average is 0.31 kWh/m<sup>3</sup> [54,59], which can be used for reference. Existing studies have also given the values of different countries for reference: 0.371–0.392 kWh/m<sup>3</sup> (USA) [56]; 0.1–0.6 kWh/m<sup>3</sup> (Australia) [60]; 0.38–1.44 kWh/m<sup>3</sup> (Canada) [61]; 0.11–1.5 kWh/m<sup>3</sup> (Spain) [62]; 0.15–0.44 kWh/m<sup>3</sup> (New Zealand) [63].

(5) Seawater Desalination (WRDB5): Nine coastal provinces in China have large-scale seawater desalination capacity. Although the industrialization process of desalination in China is slow, seawater desalination is an important behavior in the process of sustainable development of water resources in the future [51].

$$E\_5 = Q\_5 \times EI\_5 \times EF \tag{8}$$

where *E*<sup>5</sup> is the CEE of WRDB5, kg; *Q*<sup>5</sup> is the volume of seawater desalination, m<sup>3</sup> ; *EI*<sup>5</sup> is the energy intensity of WRDB5 (the amount of electricity required to treat per unit of seawater) kWh/m<sup>3</sup> ; *EI*<sup>5</sup> should be obtained based on the survey data of desalination plants, and can also refer to existing studies: 5.9 kWh/m<sup>3</sup> (China) [64–66]; 4 kWh/m<sup>3</sup> (Australia) [57]; 2.4–8.5 kWh/m<sup>3</sup> (Global) [55,67].

#### *2.3. CEEA Method to Water Resource Allocation Behaviors*

Water resource allocation behaviors (WRABs) refer to a series of activities related to water resource transportation and distribution. Representative WRABs include urban-rural tap water allocation (WRAB1) and inter-regional water transfer (WRAB2).

(1) Tap water allocation (WRAB1): The treated water from the waterworks is distributed to individual water users through the urban and rural water distribution system. The energy consumption of WRAB1 is mainly the head loss in the water transmission and distribution process, and the CEE is focused on the power consumption in the pressurization process [68]:

$$E\_6 = Q\_6 \times EI\_6 \times EF \tag{9}$$

$$EI\_6 = \frac{9.8 \times \rho \times (h\_f + h\_j)}{3.6 \times 10^6 \times \eta} \tag{10}$$

$$h\_f = \lambda \frac{l}{4R} \frac{v^2}{2g} \tag{11}$$

$$h\_{\dot{\jmath}} = \zeta \frac{v^2}{2g} \tag{12}$$

where *E*<sup>6</sup> is the CEE of WRAB1, kg; *Q*<sup>6</sup> is the amount of urban and rural tap water allocation, m<sup>3</sup> ; *EI*<sup>6</sup> is the energy intensity of WRAB1 (the amount of electricity required to distribute per unit of tap water) kWh/m<sup>3</sup> ; *h<sup>f</sup>* is head loss along the path, *λ* is the drag coefficient along the path, *l* is the length of tap water allocation, R is the hydraulic radius, m; *v* is the average velocity of tap water transmission and distribution, m3/s; *h<sup>j</sup>* is local head loss, m; *ζ* is local drag coefficient; *η* is the efficiency of the pressurized pump station. Head loss can be calculated by the Darcy formula. *EI*<sup>6</sup> can be obtained according to the investigation and statistics of unit water distribution power consumption data of water supply company. The energy intensity of tap water companies in different regions of China is quite different [69]. Combined with the China Urban Water Supply Yearbook and related research [68–70], the recommended value is 0.2 kWh/m<sup>3</sup> for reference. Reference values of *EI*<sup>6</sup> in other countries: 0.2–0.32 kWh/m<sup>3</sup> (California, USA) [26]; 0.12–0.22 kWh/m<sup>3</sup> (Spain) [71]; 0.1 kWh/m<sup>3</sup> (South Africa) [72].

(2) Inter-regional water transfer (WRAB2): Most of the inter-regional water transfer projects require pumping stations for pressurized delivery to overcome the energy loss from head loss. The CEE calculation principle of WRAB2 is similar to that of WRAB1. The difference is that the urban and rural tap water allocation system is mostly pressure pipe flow, while the inter-regional water transfer is mostly open channel constant flow:

$$E\_7 = Q\_7 \times EI\_7 \times EF \tag{13}$$

$$EI\_7 = \frac{\rho \times \text{g} \times (h\_f + h\_{\bar{f}})}{3.6 \times 10^6 \times \eta} \tag{14}$$

where *E*<sup>7</sup> is the CEE of WRAB2, kg; *Q*<sup>7</sup> is the amount of water transferred across regions, m<sup>3</sup> ; *EI*<sup>7</sup> is the energy intensity of WRAB2 (the amount of electricity required to transfer per unit water resources across regions) kWh/m<sup>3</sup> ; *h<sup>f</sup>* and *h<sup>j</sup>* are the head loss along the open channel and local head loss, m; the specific calculation can be referred to the relevant formula of open channel hydraulics [73]. In the absence of the necessary investigation conditions, *EI*<sup>7</sup> can refer to the value of 0.815 kWh/m<sup>3</sup> (China) taken in existing studies [54,74].

#### *2.4. CEEA Method to Water Resource Utilization Behaviors*

Water resource utilization behaviors (WRUBs) refer to a series of activities related to water use. WRUBs include domestic water utilization (WRUB1), industrial water utilization (WRUB2), agricultural water utilization (WRUB3), ecological water utilization (WRUB4), and hydroelectric power generation (WRUB5). Many carbon emission studies based on LCA methods do not consider the end-use process, because the emission effects caused by end-use are not part of the life cycle [29]. However, some indirect emission effects closely related to WRBs are generated or caused by these behaviors, and end-use often results in a high proportion of CEE [27]. Therefore, based on the definition of WRBs, this study also includes CEE in the end-use process of water resources in the calculation range.

(1) Domestic water utilization (WRUB1): WRUB1 does not include public domestic water because the end-use purpose of public domestic water is so broad that it is difficult to achieve a relatively accurate quantification. The main source of CO<sup>2</sup> emissions from WRUB1 is the energy consumption of the heating process [27]. Combined with the actual domestic water consumption in China, CO<sup>2</sup> emissions in the energy-consuming process of cooking and bath heating can be taken as the CEE of WRUB1, and its CEEA method is as follows:

$$E\_8 = Q\_8 \times EI\_8 \times EF \tag{15}$$

$$\rm{EI}\_8 = \rho \times R\_{household} \times (R\_{heat1} + R\_{heat2}) \times C\_w \times \Delta T \times 1/\eta \tag{16}$$

where *E*<sup>8</sup> is the CEE of WRUB1, kg; *Q*<sup>8</sup> is the total amount of domestic water consumption, m<sup>3</sup> ; *ρ* is the density of surface water (typically 1000), kg/m<sup>3</sup> ; *C<sup>w</sup>* is the heat capacity of the water (generally 1.162 <sup>×</sup> <sup>10</sup>−<sup>3</sup> kWh/(kg· ◦C) [74]); ∆*T* is the temperature difference before and after heating, ◦C; *η* is the efficiency of the heating equipment (generally 95% [74]). *Rhousehold* is the proportion of residential household domestic water consumption in total domestic water consumption; *Rheat* is the proportion of water used for heating in residential household domestic water consumption, where *Rheat*<sup>1</sup> is the proportion of cooking and drinking water, and *Rheat*<sup>2</sup> is the proportion of bathing water. Depending on different research needs, *Rhousehold* can be obtained according to the actual investigation, or according to the proportion in the water resources bulletin. In addition, studies have examined the energy intensity (*EI*8) of household water use in different regions for reference: 7.43 kWh/m<sup>3</sup> (China) [75], 24.6 kWh/m<sup>3</sup> (Ontario, Canada) [61].

(2) Industrial water utilization (WRUB2): China has a wide range of industrial sectors, and the water use processes in different sectors have different CO<sup>2</sup> emission characteristics. The energy consumption of WRUB2 is mainly concentrated in the link of water cooling and water heating [59], which is also the main source of CO<sup>2</sup> emission. There are two ideas for calculating the CEE of WRUB2:

$$E\_{\mathfrak{P}} = Q\_{\mathfrak{P}} \times EI\_{\mathfrak{P}} \times EF \tag{17}$$

$$E\_9 = \mathbb{C}\_{industry} \times R\_{water} \times EF \tag{18}$$

where *E*<sup>9</sup> is the CEE of WRUB2, kg; *Q*<sup>9</sup> is the total amount of industrial water consumption, m<sup>3</sup> ; *EI*<sup>9</sup> is the energy intensity of WRUB2 (energy consumption per unit of industrial water) kWh/m<sup>3</sup> . *EI*<sup>9</sup> can be determined from field surveys, and relevant studies have concluded that the energy intensity of industrial water use in a typical Chinese city is 5.033 kWh/m<sup>3</sup> [76]. Another idea is to calculate CEE by determining the power consumption of WRUB2 through the power consumption structure of the industrial sector [59]. A study suggests that water-related electricity consumption in the industrial sector in typical Chinese cities accounts for about 10% [59]. *Cindustry* is total industrial electricity consumption, kWh; *Rwater* is the ratio of water cooling and water heating power consumption to total power consumption in the industrial sector, %.

(3) Agricultural water utilization (WRUB3): Unlike domestic and industrial water, CO<sup>2</sup> emissions from agricultural water utilization are mainly concentrated in the irrigation process. There are five main sources of carbon emissions from farmland ecosystems: chemical fertilizers, pesticides, agricultural films, agricultural machinery, and agricultural irrigation [77]. In this study, CO<sup>2</sup> emissions from agricultural irrigation are used as the CEE of WRUB3. In addition, the carbon sink effect occurs on farmland due to photosynthesis during crop growth [78]. Therefore, the CO<sup>2</sup> absorption effect of WRUB3 should be considered [79]. The three elements of crop growth are: sunlight, water, and fertilizer, and the carbon sink effect in farmland is the result of the joint action of these three elements. Obviously, it is not appropriate to consider the entire amount of CO<sup>2</sup> absorbed by the farmland as the CO<sup>2</sup> absorption effect of WRUB3. Therefore, the CO<sup>2</sup> absorption effect of WRUB3 is separated from the overall CO<sup>2</sup> absorption effect of farmland by setting weights. Assuming that the three elements of sunlight, water, and fertilizer are equally important for the crop growth process [80], the contribution of these three elements to the carbon sink effect can be distributed by equal weight method. Of course, the weight distribution scheme can be discussed and adjusted according to the actual situation of crop planting. The CEE calculation method of WRUB3 is as follows:

$$E\_{10} = E\_{10\text{emission}} - E\_{10\text{absorption}}\tag{19}$$

$$E\_{10emision} = A \times \delta\_e \times \frac{44}{12} \tag{20}$$

$$E\_{10absorption} = \omega \times A \times \delta\_a \times \frac{44}{12} \tag{21}$$

where *E*<sup>10</sup> is the CEE of WRUB3, kg; *E*10*emission* is the total CO<sup>2</sup> emissions of WRUB3, kg; *E*10*absorption* is the amount of CO<sup>2</sup> absorbed by agricultural water utilization; *A* is the actual agricultural irrigation area, ha; *δ<sup>e</sup>* and *δ<sup>a</sup>* are CO<sup>2</sup> emission and absorption coefficient per unit irrigated area, t/ha; *ω* is the weight, which is initially set to 1/3.

(4) Ecological water utilization (WRUB4): Water resources are the foundation and core of ecosystem functions. The function of ecosystems such as woodlands, grasslands, wetlands, and watersheds cannot be performed without the maintenance of ecological water [81]. WRUB4 refers to artificial ecological water, that is, urban environmental water and rivers, lakes, and wetland replenishment water supplied by human measures [82]. Different from domestic and production water utilization behaviors, the CEE of WRUB4 cannot be directly quantified by energy as a medium. Therefore, in this study, the CO<sup>2</sup> absorbed by four land types closely related to ecological water use, namely, urban garden, urban green space (excluding garden area), water area within the jurisdiction, and wetland within the jurisdiction, is roughly taken as the CEE of WRUB4. Of course, the actual process of CO<sup>2</sup> absorption from WRUB4 is far more complicated than described.

$$E\_{11} = -\sum\_{i}^{n} A\_i \times \delta\_i \times \frac{44}{12} \tag{22}$$

where *E*<sup>11</sup> is the CEE of WRUB4, kg; *A* is the area of ecological water land type, ha; *δ* is the CO<sup>2</sup> absorption coefficient of ecological water land type (the amount of CO<sup>2</sup> absorbed per unit area of ecological water land), t/ha. *i* is the type of land. *δ* can be obtained by field measurements in the study area, or by referring to existing studies [78].

(5) Hydroelectric power generation (WRUB5): CO<sup>2</sup> emissions from hydropower generation are much lower than those from thermal power [83]. Based on the UN CDM (United Nations' Clean Development Mechanism), GHG emissions from hydropower generation can be disregarded in the calculation of hydropower CDM projects [84]. Therefore, the relative carbon reduction effect of hydropower compared to thermal power is used in this study to quantify the CEE of WRUB5.

$$E\_{12} = -G \times \mathcal{C}PG \times EF\_{\mathcal{C}} \tag{23}$$

where *E*<sup>12</sup> is the CEE of WRUB5, kg; *G* is the total amount of hydroelectric power, kWh; *CPG* is the standard coal consumption of power generation unit, tce/kWh; *EF<sup>c</sup>* is the CO<sup>2</sup> emission coefficient of standard coal, kg/tce. *CPG* can be obtained from the investigation of the thermal power industry in the study area. Studies have shown that the average coal consumption of thermal power generating units in China is 3.7 <sup>×</sup> <sup>10</sup>−<sup>4</sup> tce/kWh [85]. *EF<sup>c</sup>* can refer to IPCC guidelines for national greenhouse gas inventories [52] or existing studies [85].

#### *2.5. CEEA Method to Water Resource Protection Behaviors*

Water resource protection behaviors (WRPBs) refer to a series of activities related to water resources protection, including water saving (WRPB1), wastewater collection (WRPB2), wastewater treatment (WRPB3), and reclaimed water reuse (WRPB4).

(1) Water saving (WRPB1): Water saving behavior directly avoids part of the energy consumed in the development and allocation of water resources, so it can be regarded as a carbon reduction behavior [32,86]. Its CEEA method is as follows:

$$E\_{13} = -Q\_{13} \times \left(EP\_{exploitation} + EP\_{distribution}\right) \tag{24}$$

$$EP\_{exploitation} = (E\_1 + E\_2) / (Q\_1 + Q\_2) \tag{25}$$

$$EP\_{\text{distribution}} = \text{E}\_7/\text{Q}\_7\tag{26}$$

where *E*<sup>13</sup> is the CEE of WRPB1, kg; *Q*<sup>13</sup> is the total amount of water saved, m<sup>3</sup> ; *EPexploitation* is the comprehensive CO<sup>2</sup> emission coefficient of water resource exploitation (CO<sup>2</sup> emissions per unit of water resource exploitation), kg/m<sup>3</sup> ; *EPdistribution* is the comprehensive

CO<sup>2</sup> emission coefficient of water resource allocation (CO<sup>2</sup> emissions per unit of water resource allocation), kg/m<sup>3</sup> . Other variables have the same meaning as above.

(2) Wastewater collection (WRPB2): Wastewater from different sources usually relies on gravity to converge to the wastewater network, and then is pressurized by the wastewater network pump to the wastewater treatment plant. Similar to WRAB1, the CEE of WRPB2 is mainly generated by energy consumption to overcome head loss [49]:

$$E\_{14} = Q\_{14} \times EI\_{14} \times EF \tag{27}$$

$$EI\_{14} = \frac{9.8 \times \rho \times (h\_f + h\_{\bar{j}})}{3.6 \times 10^6 \times \eta} \tag{28}$$

where *E*<sup>14</sup> is the CEE of WRPB2, kg; *Q*<sup>14</sup> is the total amount of wastewater collected, m<sup>3</sup> ; *EI*<sup>14</sup> is the energy intensity of WRPB2 (electricity consumption by collecting unit of wastewater), kWh/m<sup>3</sup> . *EI*<sup>14</sup> should be obtained based on the investigation and statistics of the wastewater collection system in the study area, and can also refer to the values in related studies: 0.013 kWh/m<sup>3</sup> (China) [86].

(3) Wastewater treatment (WRPB3): The treatment methods of wastewater treatment plants in different countries are different, but generally include three stages: primary treatment, secondary treatment, and tertiary treatment. Each stage has different processes, and the energy consumption intensity of each process is different. The main CO<sup>2</sup> emissions are concentrated in the secondary and tertiary treatment stages [87]. On the other hand, untreated wastewater contains more pollutants such as *COD* and *BOD*5, which can produce large amounts of carbon emissions. WRPB3 has a positive CO<sup>2</sup> reduction effect by reducing the concentration of such pollutants [88]. In addition, wastewater treatment plants generally use sludge in wastewater for power generation [89], and its carbon reduction effect should also be considered. In this study, the CO<sup>2</sup> absorption effect of WRPB3 is considered based on the concentration difference of major carbon emission pollutants before and after wastewater treatment and the sludge power generation:

$$E\_{15} = E\_{15
emission} - E\_{15
absorption} \tag{29}$$

$$E\_{15emission} = Q\_{15} \times EI\_{15} \times EF - Q\_{15} \times R\_s \times P\_s \times EF \tag{30}$$

$$EI\_{15} = \sum\_{i=1}^{3} \sum\_{j} EI\_{ij} \tag{31}$$

$$E\_{15\text{absorption}} = Q\_{15} \times \Delta R\_{\text{COD}} \times EF\_{\text{COD}} + Q\_{15} \times \Delta R\_{\text{BOD5}} \times EF\_{\text{BOD5}}\tag{32}$$

where *E*<sup>15</sup> is the CEE of wastewater treatment behavior, kg; *Q*<sup>15</sup> is the total amount of wastewater treatment, m<sup>3</sup> ; *EI*<sup>15</sup> is the energy intensity of WRPB3 (electricity consumption by treating unit of wastewater), kWh/m<sup>3</sup> . *EIij* is the energy consumption intensity of the process *j* in stage *i*, kWh/m<sup>3</sup> . The energy intensity or emission factor of unit wastewater treatment can be obtained by investigating the energy consumption and treatment capacity of the wastewater treatment plant [28,29]. *EI*<sup>15</sup> from relevant studies are available for reference: 0.24 kWh/m<sup>3</sup> (China) [74]; 0.8–1.5 kWh/m<sup>3</sup> (Australia) [60]; 0.177–0.78 kWh/m<sup>3</sup> (USA) [56]; 0.41–0.61 kWh/m<sup>3</sup> (Spain) [71]; 0.44 kWh/m<sup>3</sup> (South Africa) [72]; 0.38–1.122 kWh/m<sup>3</sup> (Global) [55]. *Rs* is the sludge concentration in wastewater, generally 0.3–0.5% [90]; *Ps* is the power generation of unit sludge, and the coefficient in related research is 14.27 kWh/m<sup>3</sup> for reference [89]. ∆*RCOD* and ∆*RBOD*<sup>5</sup> are the concentration differences of *COD* and *BOD*5 before and after wastewater treatment, respectively. When the measurement conditions are available, the measurement results shall prevail. When conducting large-scale research, ∆*RCOD* and ∆*RBOD*<sup>5</sup> can also be determined according to relevant emission standards. According to China's comprehensive wastewater discharge standard, the concentration difference between *COD* and *BOD*5 before wastewater treatment (Level 3 standard) and after wastewater treatment (Level 1 standard) is 0.94 kg/m<sup>3</sup> and 0.58 kg/m<sup>3</sup> . *EFCOD* and *EFBOD*<sup>5</sup> are the amount of CO<sup>2</sup> reduced by removing unit *COD* and *BOD*5, and the units are kgCO2/kg*COD* and kgCO2/kg*BOD*5, respectively. According to the relevant emission factors released by IPCC [52], *EFCOD* and *EFBOD*<sup>5</sup> are 0.69 and 1.65, respectively.

(4) Reclaimed water reuse (WRPB4): Reclaimed water reuse reduces the extraction of surface water and groundwater, and can therefore be considered as a WRB to reduce CO<sup>2</sup> emissions. The calculation formula of CEE is as follows:

$$E\_{16} = -Q\_{16} \times EP\_{exploitation} \tag{33}$$

$$EP\_{exploitation} = (E\_1 + E\_2) / (Q\_1 + Q\_2) \tag{34}$$

where *E*<sup>16</sup> is the CEE of reclaimed water reuse behavior, kg; *Q*<sup>16</sup> is the amount of reclaimed water reuse, m<sup>3</sup> ; *EP* is the comprehensive CO<sup>2</sup> emission coefficient of water resources exploitation (CO<sup>2</sup> emissions per unit of water resource exploitation), kg/m<sup>3</sup> .

#### *2.6. Function Table of CEEA for Water Resource Behaviors*

The above methods and ideas are summarized and all the CEEA formulas are combined to form a table, which is the function table of CEEA (FT-CEEA) for WRBs (Table 2). In addition, in view of the large regional differences in the grid CO<sup>2</sup> emission factor and the energy intensity of groundwater extraction, the referenceable values (Tables 3 and 4) for different regions of China are given [54], which can be selected according to the actual situation of the study area. The instructions for using FT-CEEA are as follows.



**Table 2.** *Cont.*

**Table 3.** Average CO<sup>2</sup> emission factor of power grids in different regions of China (kgCO2/kWh).


**Table 4.** Energy intensity of unit groundwater extraction in different regions of China (kWh/m<sup>3</sup> ).


(1) FT-CEEA is a collection of formulas for estimating and cross-sectionally comparing the CEE of various WRBs. The CEEA formulas for different WRBs in FT-CEEA can be used selectively depending on the study purpose and study scale. The quantity, type, and calculation method of WRBs in FT-CEEA are not static and can be updated and improved according to the changing situation and new research progress.

(2) The results of each formula are not necessarily an absolute measurement of the emission or absorption effects of CO2, but the idea of each formula is relatively reasonable. FT-CEEA is equivalent to setting up a "ruler" as a relative comparison of CEE generated by WRBs calculated by different researchers. FT-CEEA has no scale limitation and can be applied to different scales with limited accuracy requirements. However, the specific parameters need to be adjusted according to the actual situation of the research object.

(3) Most of the formulas in FT-CEEA need to be supported by relevant parameters, but in most cases, it is difficult to carry out field investigations and measurements of the parameters. Given this situation, some valuable reference values are provided in this table. Of course, some changes can be made in the selection of parameter reference values according to different research needs and actual conditions.

#### **3. Case Study**

#### *3.1. Overview of the Study Area*

China has a vast territory, and there are significant spatial differences in industrial structure, water use mode, and carbon emission intensity in different regions. In terms of CO<sup>2</sup> emissions in 2019, Shanxi (the province with the highest emission intensity) is 37 times higher than Qinghai (the province with the lowest emission intensity) under different development orientation [91]. In the past 20 years, under the background of rapid economic and social development, some provinces in China are facing many challenges such as the insufficient capacity for sustainable utilization of water resources and prominent conflict between carbon emission reduction and economic development [92]. Since the 1990s, China has been in a new period of rapid growth in carbon emissions, lagging behind developed countries in time. Although China's total carbon emissions ranked first in the world in recent years, China's per capita carbon emissions are still far lower than developed countries. Many traditional industries in China still maintain a production mode with high consumption and high emission. Promoting the low-carbon transformation of traditional industries has become an urgent bottleneck to achieving China's carbon neutrality goal [7].

In this study, 31 provincial administrative regions in mainland of China are divided into 8 regions [93]. The regional division, elevation distribution, water supply structure, and CO<sup>2</sup> emission intensity of the study area are shown in Figure 3.

#### *3.2. Data source and Description*

In addition to the important parameters in FT-CEEA, the data used in the case study are mainly the data of indicators involved in different WRBs of 31 provinces in China in 2020. The data involved in WRABs include tap water allocation and inter-regional water transfer. The data involved in WRUBs include domestic water consumption, industrial water consumption, actual agricultural irrigation area, land area of four kinds of artificial ecological water utilization, and hydroelectric power generation. The data involved in WRPBs include water saving, wastewater treatment, and reclaimed water reuse.

The sources of the above data include China Water Resources Bulletin 2020, China Seawater Utilization Bulletin 2020, Water Resources Bulletin of 31 provinces in 2020, China Statistical Yearbook 2021, China Water Statistical Yearbook 2021, China Energy Statistical Yearbook 2021, China Environmental Statistical Yearbook 2021, and China Urban Construction Statistical Yearbook 2021.

the 1990s, China has been in a new period of rapid growth in carbon emissions, lagging behind developed countries in time. Although China's total carbon emissions ranked first in the world in recent years, China's per capita carbon emissions are still far lower than developed countries. Many traditional industries in China still maintain a production mode with high consumption and high emission. Promoting the low-carbon transformation of traditional industries has become an urgent bottleneck to achieving China's car-

In this study, 31 provincial administrative regions in mainland of China are divided into 8 regions [93]. The regional division, elevation distribution, water supply structure,

and CO<sup>2</sup> emission intensity of the study area are shown in Figure 3.

**Figure 3.** The study area. **Figure 3.** The study area.

bon neutrality goal [7].

#### *3.2. Data source and Description 3.3. Results and Discussion*

In addition to the important parameters in FT-CEEA, the data used in the case study 3.3.1. Carbon Dioxide Emission Equivalent Analysis of WRDBs

are mainly the data of indicators involved in different WRBs of 31 provinces in China in 2020. The data involved in WRABs include tap water allocation and inter-regional water transfer. The data involved in WRUBs include domestic water consumption, industrial Based on FT-CEEA and the above data, the CEE of WRBs in 31 provinces and 8 regions of China in 2020 was calculated. The calculation results of the eight regions are obtained by summing the included provinces.

water consumption, actual agricultural irrigation area, land area of four kinds of artificial ecological water utilization, and hydroelectric power generation. The data involved in WRPBs include water saving, wastewater treatment, and reclaimed water reuse. The sources of the above data include China Water Resources Bulletin 2020, China Seawater Utilization Bulletin 2020, Water Resources Bulletin of 31 provinces in 2020, China Statistical Yearbook 2021, China Water Statistical Yearbook 2021, China Energy Statistical Yearbook 2021, China Environmental Statistical Yearbook 2021, and China Urban Construction Statistical Yearbook 2021. The CEEA results of WRDBs are presented in Table 5. In 2020, the surface water lifting behavior (WRDB1) in eight regions of China generated 63.52 million tons of CEE, accounting for 29.8% of the total CEE produced by WRDBs. Among them, the WRDB1 in middle Yangtze River and east coast provinces produced higher CEE of 12.63 million tons and 11.9 million tons, respectively. The three provinces of Shanghai, Jiangsu, and Zhejiang in the east coast region are dominated by surface water utilization. The surface water supply of Jiangsu Province in 2020 is 55.6 billion cubic meters, resulting in the CEE generated by WRDB1 ranking first among 31 provinces (8.18 million tons). The region with the smallest CEE of WRDB1 is the north coast region (4.64 million tons). On the other hand, WRDB2 in the north coast region produced the most CEE (8.2 million tons). In contrast, groundwater extraction in the east coast region produced only 0.12 million tons of CEE in 2020. The spatial distribution characteristics of CEEA results of WRDB1 and WRDB2 are closely related to the water supply structure in different regions. Compared with the southern provinces of China, the northern provinces have a higher degree of groundwater exploitation and a larger proportion of groundwater utilization, which is also a manifestation of the uneven spatial distribution of water resources in China [94]. In addition, Xinjiang is the province with the most CEE generated by WRDB2 in 31 provinces (5.69 million tons). The reason is that Xinjiang has a large amount of groundwater supply. In 2020, the groundwater supply in Xinjiang is 12.43 billion cubic meters, second only to Heilongjiang (12.94 billion cubic meters). Another important factor is that Xinjiang's higher altitude means it takes much more energy to extract per unit of groundwater than the eastern provinces [39].


**Table 5.** CEE of WRDBs and WRABs in eight regions of China in 2020 (10,000 tons).

WRDB3 is the behavior that produces the most CEE in WRDBs, generating 75.9 million tons of CEE in 2020, accounting for 35.6% of the total CEE produced by WRDBs. Among them, the CEE produced in middle Yangtze River and southwest regions was significantly higher than that in other regions, and the CEE produced by WRDB3 in the east coast region was less (4.69 million tons). Raw water treatment behavior (WRDB4) produced 39.55 million tons of CEE in 2020. Due to the high proportion of domestic and industrial water, the east coast, the middle Yangtze River and the southern coastal provinces have become the main contributors to the CEE generated by WRDB4. In 2020, the CEE generated by seawater desalination behavior (WRDB5) was 2.92 million tons, accounting for 1.4% of the total CEE generated by WRDBs. China's desalination plants are mainly concentrated in 9 coastal provinces [64], which are Shandong, Hebei, Zhejiang, Tianjin, Liaoning, Guangdong, Fujian, Hainan, and Jiangsu in descending order according to CEE. The proportion of CEE generated by WRDB5 in the north and east coast provinces exceeded 87%.

#### 3.3.2. Carbon Dioxide Emission Equivalent Analysis of WRABs

The CEEA results of water resource allocation behaviors (WRABs) are shown in Table 5. WRAB1 produced 26.36 million tons of CEE in 2020, accounting for 60.9% of the total CEE produced by WRABs. The CEE of WRAB1 is similar to WRDB4 in spatial distribution. The difference in water resources utilization structure in different regions of China can explain the distribution characteristics to some extent. Compared with the eastern provinces of China, the northwest provinces have a higher proportion of agricultural water and a lower proportion of industrial and domestic water [95]. Tap water supply is mainly concentrated in industrial and domestic water. Therefore, the water use structure dominated by agricultural water has led to the CO<sup>2</sup> emission effect of WRAB1 in the northwest region being much lower than that in the eastern region. Cross-regional water transfer behavior produced 16.95 million tons of CEE in 2020. Due to the existence of large-scale water diversion projects such as the South-to-North Water Diversion Project and the Luanhe River Diversion Project, the CEE generated by WRAB2 in the north coast and the middle Yellow River provinces accounted for up to 89%. This spatial distribution feature is similar to the research results of Xiang and Jia [54].

#### 3.3.3. Carbon Dioxide Emission Equivalent Analysis of WRUBs

The CEEA results of WRUBs are shown in Table 6. Among the five kinds of WRUBs, the CEE value of domestic water utilization and industrial water utilization is positive, resulting in the CO<sup>2</sup> emission effect. The CEE value of agricultural water utilization, ecological water utilization, and hydroelectric power generation is negative, resulting in the CO<sup>2</sup> absorption effect. Among them, the CEE calculation of WRUB2 is based on the first calculation scheme (energy intensity scheme).


**Table 6.** CEE of WRUBs in eight regions of China in 2020 (10,000 tons).

In 2020, the CEE of WRUB1 (326.88 million tons) and WRUB2 (353.54 million tons) in 31 provinces of China are not very different in total, but there are large differences between regions. The CEE generated by WRUB2 in the east coast and middle Yangtze River provinces is higher than that generated by WRUB1, especially in the east coast provinces. The outstanding proportion of industrial and domestic water in Jiangsu Province leads to the highest CEE generated by WRUB2. The difference in water use structure is the main reason for the difference in CEE of industrial and domestic water in different regions [96]. The absorption effect of WRUB3 (288.84 million tons) is greater than the emission effect (57.02 million tons), so the CEE of agricultural water utilization behavior is negative in total (−231.83 million tons). The WRUB3 of the northwest provinces has produced a considerable CO<sup>2</sup> emission effect (6.48 million tons), which is consistent with the local water resource utilization structure [97]. The middle Yangtze River provinces have more agricultural irrigation area, and the CO<sup>2</sup> absorption effect produced by WRUB3 is also the highest among the eight regions (51.9 million tons).

Ecological water utilization behavior (WRUB4) produced −144.06 million tons of CEE in 2020. The CEE of WRUB4 in southwest and northwest provinces was nearly half of the total CEE produced by WRUB4. The main reason is that the wetland and water area of Sichuan, Tibet, Qinghai, Xinjiang, and other provinces is much higher than other regions. The strong guarantee of ecological water use in the above-mentioned provinces has played an important role in maintaining the carbon sink function of wetland and water ecosystem [98]. In 2020, the hydroelectric power generation behavior (WRUB5) in eight regions of China produced a total of −335.95 million tons of CEE with significant spatial differences. Southwest provinces have the most abundant hydropower resources [99], while the proportion of hydropower in the energy structure of the north coast provinces is very small. The distribution of hydropower resources in China is the main reason for the CEE spatial difference of WRUB5.

#### 3.3.4. Carbon Dioxide Emission Equivalent Analysis of WRPBs

The CEEA results of WRPBs are shown in Table 7. Among the four WRPBs, only the CEE value of wastewater collection behavior (WRPB2) is positive, resulting in CO<sup>2</sup> emission effect. The CEE values of the other three WRPBs are negative, resulting in the CO<sup>2</sup> absorption effect. Water saving behavior (WRPB1) can undoubtedly provide a positive impact on reducing CO<sup>2</sup> emissions [100]. If only the energy saving effect of WRPB1 on water resources development and allocation is considered, the CEE of WRPB1 in 2020 is −2.05 million tons. In general, Shanghai, Guangdong, Zhejiang, Jiangsu, and Beijing are at the forefront of the construction of water-saving society [101], and there is still a large room for improvement in the capacity of water-saving and emission reduction in northwest provinces. The CEE of WRPB2 is the smallest among all WRBs in FT-CEEA (0.5 million tons). The CO<sup>2</sup> absorption effect produced by wastewater treatment behavior (WRPB3) is significantly greater than the emission effect. The spatial distribution characteristics of CEE of WRPB3 are directly related to the wastewater treatment capacity of different regions.

The east coast and south coast provinces have a large amount of wastewater discharge and a strong wastewater treatment capacity [102], which correspondingly brings a higher carbon dioxide emission and absorption effect. If only the energy saving effect of reclaimed water reuse behavior on water resources development is considered, the CEE generated by WRPB4 in 2020 is −2.68 million tons. The CEEA results of WRPB4 are closely related to regional water resource endowment and water supply structure. Compared with the southern provinces, Beijing, Hebei, Shandong, Henan, and other northern provinces are relatively short of water, so the reuse of reclaimed water has become an effective means to alleviate the contradiction between local water supply and demand [103]. As a result, the amount of reclaimed water supplied by these provinces is much higher than that of other provinces, and correspondingly, more CO<sup>2</sup> absorption effect is generated.


**Table 7.** CEE of WRPBs in eight regions of China in 2020 (10,000 tons).

#### **4. Conclusions**

In this study, the carbon dioxide emission equivalent analysis (CEEA) method of water resource behaviors (WRBs) was developed, and a function table of carbon dioxide emission equivalent (FT-CEEA) was constructed. Based on the FT-CEEA, the CEE of different WRBs in 31 provinces of China in 2020 was analyzed. Some valuable conclusions are as follows:


Based on the above conclusions, some targeted measures and suggestions are discussed for the carbon neutrality goal in the field of water resources. Increasing the proportion of hydropower generation, improving the capacity of ecological water security, strengthening wastewater treatment and reclaimed water reuse, and promoting the construction of water-saving society can be considered as effective ways to promote carbon neutrality in this field.

However, there are still some limitations. The consideration of water resource behavior categories may not be comprehensive. In this study, the water resource behaviors were divided into four categories: development, allocation, utilization, and protection. However, water resource behaviors are not limited to the four categories, and the number of WRBs is far more than 16. Therefore, FT-CEEA is dynamic rather than static, and needs to be constantly updated. In addition, many CEE calculations of WRB are completed by using energy as an intermediate medium, which is the quantitative scheme adopted by most related studies. Although the energy consumption is the major factor in the generation of CEE by those WRBs, it cannot be excluded that there may be other potential factors contributing to carbon emissions. When these potential factors reach a certain scale, the resulting CEE also needs to be considered. Moreover, for some WRBs, the CEEA method may not be considered perfect. For example, the CO<sup>2</sup> absorbed by the four types of land closely related to ecological water utilization was roughly used as the CEE of WRUB3. In fact, the CO<sup>2</sup> absorbed by the lands is due to many factors, including ecological water utilization. How to separate the CEE of ecological water and CEE produced by other factors? Further exploration and refinement are still needed.

**Author Contributions:** Conceptualization, formal analysis, funding acquisition, project administration, Q.Z.; data curation, methodology, software, visualization, writing—original draft, Z.Z., C.Z. and X.Q.; investigation, J.M., Z.Z., C.Z., and X.Q.; supervision, Q.Z. and J.M.; resources, Q.Z. and Z.Z.; writing—review and editing, Q.Z., Z.Z., J.M., C.Z. and X.Q. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Science Foundation of China (no. 52279027), the National Key Research and Development Program of China (no. 2021YFC3200201), and the Major Science and Technology Projects for Public Welfare of Henan Province (no. 201300311500).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. These data can be found and downloaded on relevant websites.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

**Yuxin Liang <sup>1</sup> , Liping Zhang 1,2,3,\*, Mengsi Leng <sup>1</sup> , Yi Xiao 1,2,3,\* and Jun Xia 1,2,3**


**Abstract:** Green development is a low-carbon, sustainable model for the achievement of the harmonious development of the economy and nature. Nowadays, the problems of resource scarcity and environmental pollution in the process of economic development are pressing, and the promotion of green development is the general trend. As one of the three growth poles of China's Yangtze River economic belt, the Chengdu-Chongqing City Group is an important platform to lead toward green development in the western region of China. Based on the understanding of the connotation of green development, this study established a green development-level evaluation system, including 19 indicators in three dimensions: target level, criterion level, and indicator level, and used the entropy weight method to measure the green development level of the Chengdu-Chongqing City Group. In view of the dynamic nature of the green development process, this study constructed a system dynamics model of the green development level of the Chengdu-Chongqing City Group and simulated and compared it between 2022 and 2050 under five shared socio-economic pathway (SSP) scenarios so as to provide a reference basis for future development. The results show that the overall green development level of the Chengdu-Chongqing City Group is on an upward trend, with the highest green development level under the SSP1 path and the lowest under the SSP3 path, and the lagging distance tends to increase further. In the next 30 years, the Chengdu-Chongqing City Group should initially follow SSP2 as the basis for development and then gradually perform a transition to SSP1 by 2035 to achieve real sustainable development, after which it should continue to develop according to the SSP1 path until 2050.

**Keywords:** urban agglomeration; green development; system dynamics; shared socio-economic pathways; simulation prediction

### **1. Introduction**

Green development theory originated from the green movement in the West. British scholars first noticed the constraint relationship between socio-economic development and environmental capacity. In 1662, William Petty recognized that value originates from labor, whose ability to create wealth is not infinite, being constrained by conditions such as land and technology [1]. Although the theory of green development has been deeply studied by international scholars for decades, there is still no uniform definition.

The research on green development theory in China is relatively late compared to that in Western countries; it focuses on reducing the use of resources and environmental damage while ensuring economic growth.

The representative one is the "people-oriented" concept of green development, proposed by Hu Angang of the National Research Center of Tsinghua University, which emphasizes the win–win situation of economic development and ecological environmental

**Citation:** Liang, Y.; Zhang, L.; Leng, M.; Xiao, Y.; Xia, J. System Simulation and Prediction of the Green Development Level of the Chengdu-Chongqing City Group. *Water* **2022**, *14*, 3947. https:// doi.org/10.3390/w14233947

Academic Editor: Bahram Gharabaghi

Received: 9 November 2022 Accepted: 30 November 2022 Published: 4 December 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

protection, i.e., the rational exploitation of resources and the reduction in pollution under the premise of sustained and stable economic growth [2]. In synthesis, green development is a model of economic development that does not entail damage to the ecological environment or over-exploitation and consumption of resources. It differs from the traditional crude development mode as it is environmentally friendly and entails resource-saving; moreover, it increases green wealth, allows the realization of a green life, and promotes the harmonious coexistence of humans and nature.

In recent years, green development, as a transformative development concept, has gradually achieved global consensus. With the maturity of the green development theory, scholars have started to conduct research related to the green development level, mainly focusing on the relationship between the economy, resources, and the environment [3]. More specifically, the related research focused mainly on three aspects: the measurement and evaluation of the green development level [3,4]; the identification of its influence factors [5,6]; and its simulation prediction [7,8]. These studies focused on both single cities and the whole country.

At present, the green development level is measured and evaluated mainly through the comprehensive evaluation index system method, the TOPSIS model, and the data envelopment analysis method. At the end of the 20th century, the comprehensive green development index system was developed and was further improved [9–13] as a comprehensive evaluation index system to measure and evaluate the green development level in specific regional and urban scopes. In 2009, the Economist Intelligence Unit (London, UK) constructed eight categories of green city indicators to analyze the environmental situation in 30 European cities [14]. In 2018, Carli et al. constructed a multi-objective integrated indicator system to analyze the sustainable development of urban energy, water, and environmental systems [15]. In 2019, Tian et al. adopted the entropy method to measure the greening of the Yangtze River Delta city agglomeration in terms of green growth, green welfare, green wealth, and green governance [16]. In 2021, Hu et al. established an urban green development evaluation index system for 108 cities in the Yangtze River Economic Zone to analyze the role of technological innovation in promoting green development [17].

To fully reflect the dynamics of the green development process and the interaction of its subsystems, as well as to improve the guidance and support for realistic decision making, scholars introduced the system dynamics approach to study green development from a system perspective. Founded by Professor Forrester at MIT, system dynamics was initially called "industrial dynamics" [18]; after the 1960s, he transformed system dynamics into urban dynamics when he studied the rise and fall of cities. This concept was then gradually improved and developed by Mass and Alfeld [19,20].

In the field of green development, the research on system dynamics is mostly focused on future simulation, and system dynamics models have been constructed to analyze the change trend of the future green development level by setting different simulation scenarios so as to obtain the optimal development path or countermeasure suggestions. Rudneva et al. constructed a "green" economic system dynamics model and analyzed the scenarios of regional "green" economic development [21]. Yang Guangming et al. used Chongqing city as an example to simulate and predict the sustainable development of the water resource carrying capacity, based on a system dynamics model [22]. Zhou studied the optimal path of urban agglomeration development by simulating the system dynamics of green development in the Beijing-Tianjin-Hebei urban agglomeration under different paths, putting forward relevant policy suggestions [23]. Wang et al. simulated and predicted the future scenarios of the ecological security index of the Beijing-Tianjin-Hebei city agglomeration by constructing a system dynamics model [24]. Li et al. considered three scenarios of economic development, technological innovation, and government investment to simulate and predict the industrial green development of the Beijing-Tianjin-Hebei urban agglomeration [25]. Jing et al. developed a system dynamics model of green lowcarbon development for the Beijing-Tianjin-Hebei urban agglomeration considering four factors: economy, energy, environment, and urban aggregation, to explore the future system

effects under different scenarios [26]. O'Keeffe et al. proposed a natural capital system dynamics model framework to predict and assess the natural capital performance of urban development in London [27]. Gudlaugsson, B et al. apply a system dynamics approach to investigate the complexity of energy transition in the Tees Valley region of the UK under a sustainable transition pathway and make policy recommendations [28]. In synthesis, a large number of studies have shown the general applicability of system dynamics to perform simulation and prediction studies in the field of green development.

Through a literature review, it can be found that the simulation and prediction research on green development is relatively mature, although there are still some areas for improvement:


To address these shortcomings, this study selected the Chengdu-Chongqing City Group as the study area and divided the whole green development level system into four subsystems: social, economic, resource, and environmental subsystems. Next, a green development level system dynamics model was constructed, whose reliability and stability was verified through a historical test and a sensitivity test. On this basis, a simulation scheme of the system dynamics model of each city in the Chengdu-Chongqing City Group from 2019 to 2050 under five shared socio-economic pathways (SSPs) was set to analyze the change trend of the green development level, select the optimal pathway, and put forward corresponding suggestions and countermeasures.

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

#### *2.1. Description of the Study Area*

The Chengdu-Chongqing City Group, the Yangtze River midstream urban agglomeration, and the Yangtze River Delta urban agglomeration are the three major growth poles of China's Yangtze River economic belt, and the core areas for the promotion of urbanization and economic growth. The Chengdu-Chongqing City Group is located in the upper reaches of the Yangtze River; it spans between 28◦240 N~32◦270 N and 101◦560 E~110◦110 E, and it is a region with both development strength and potential in western China. As such, it is an important basis for the green development of the Yangtze River economic belt. The core of the Chengdu-Chongqing City Group is represented by the cities of Chengdu and Chongqing; it includes 27 districts (counties) in Chongqing and parts of Kaixian and Yunyang and 15 cities in Sichuan Province [29]. The location and extension of the Chengdu-Chongqing City Group is shown in Figure 1.

**Figure 1.** The Chengdu-Chongqing City Group scope diagram. **Figure 1.** The Chengdu-Chongqing City Group scope diagram.

In terms of natural resources, the Chengdu-Chongqing City Group has abundant precipitation, with an average annual precipitation of about 900~1300 mm and abundant underground thermal energy and drinking water, such that the groundwater and transit water can basically meet the residents' water needs for living and production. At the same time, there are 12 main streams, including the Minjiang River and the Tuo River, and dozens of tributaries; hence, there are sufficient natural water resources in the Chengdu-Chongqing City Group area. According to the incomplete statistics alone, the animal and plant resources reach 11 classes, 200 families, 764 genera, and more than 3000 species, including rare animals such as giant pandas, red pandas, and golden snub-nosed monkeys, which are protected by the state. Thus, it can be seen that the Chengdu-Chongqing City Group is rich in natural resources such as water resources, mineral resources, land In terms of natural resources, the Chengdu-Chongqing City Group has abundant precipitation, with an average annual precipitation of about 900~1300 mm and abundant underground thermal energy and drinking water, such that the groundwater and transit water can basically meet the residents' water needs for living and production. At the same time, there are 12 main streams, including the Minjiang River and the Tuo River, and dozens of tributaries; hence, there are sufficient natural water resources in the Chengdu-Chongqing City Group area. According to the incomplete statistics alone, the animal and plant resources reach 11 classes, 200 families, 764 genera, and more than 3000 species, including rare animals such as giant pandas, red pandas, and golden snub-nosed monkeys, which are protected by the state. Thus, it can be seen that the Chengdu-Chongqing City Group is rich in natural resources such as water resources, mineral resources, land resources, and biological resources.

resources, and biological resources. In terms of socio-economics, in 2018 the Chengdu-Chongqing City Group had a resident population of 95 million, accounting for 6.8% of the total population of the country. This area is rich in human resources and is a multi-ethnic area with more than 50 ethnic minorities, including Tibetans, Hui, Yi, and Tujia, in addition to the Han Chinese. Moreover, it is one of the regions with the best economic foundation and the strongest comprehensive strength in western China. Chengdu was awarded the first batch of national culture and tourism consumption demonstration cities in 2020; Chongqing has formed the world's largest electronic information industry cluster and China's largest domestic automotive industry cluster; the Chengdu and Chongqing region has strong strength in electronic information, equipment manufacturing, finance, and other industries; it has a better environment for innovation and entrepreneurship, and the open economic system is grad-In terms of socio-economics, in 2018 the Chengdu-Chongqing City Group had a resident population of 95 million, accounting for 6.8% of the total population of the country. This area is rich in human resources and is a multi-ethnic area with more than 50 ethnic minorities, including Tibetans, Hui, Yi, and Tujia, in addition to the Han Chinese. Moreover, it is one of the regions with the best economic foundation and the strongest comprehensive strength in western China. Chengdu was awarded the first batch of national culture and tourism consumption demonstration cities in 2020; Chongqing has formed the world's largest electronic information industry cluster and China's largest domestic automotive industry cluster; the Chengdu and Chongqing region has strong strength in electronic information, equipment manufacturing, finance, and other industries; it has a better environment for innovation and entrepreneurship, and the open economic system is gradually taking shape.

ually taking shape. In terms of spatial structure, the Chengdu-Chongqing City Group is located at the intersection of the horizontal axis along the Yangtze River corridor and the vertical axis of the Baokun corridor, with the Chengdu metropolitan area and the Chongqing metropolitan area as the main development axis, complemented by the urban belt along the river and the cities of Chengdu, Deyang, Mianyang, and Leshan, forming a spatial development pattern characterized by "one axis, two belts, two cores and three districts" [30]. As such, this region is key for the rapid development of central and western China and In terms of spatial structure, the Chengdu-Chongqing City Group is located at the intersection of the horizontal axis along the Yangtze River corridor and the vertical axis of the Baokun corridor, with the Chengdu metropolitan area and the Chongqing metropolitan area as the main development axis, complemented by the urban belt along the river and the cities of Chengdu, Deyang, Mianyang, and Leshan, forming a spatial development pattern characterized by "one axis, two belts, two cores and three districts" [30]. As such, this region is key for the rapid development of central and western China and the creation of new space for economic growth [31].

the creation of new space for economic growth [31].

#### *2.2. Data Sources*

The statistical data used in this article are mainly from the 1995–2018 "China Statistical Yearbook" (http://www.stats.gov.cn/tjsj./ndsj/) (accessed on 17 February 2022), released by the National Bureau of Statistics of China; the "Sichuan Statistical Yearbook" (http://tjj.sc.gov.cn/scstjj/c105855/nj.shtml) (accessed on 17 February 2022); the statistical yearbooks of 16 cities in the Chengdu-Chongqing City Group; the statistical bulletin of national economic and social development; the water resources bulletin; and other documents published in recent years. The interpolation method was employed to fill in the data gaps for single years.

#### *2.3. The Entropy Weight Method*

The entropy weight method is one of the most widely used objective weighting methods. This method determines the weights through a certain mathematical operation according to the relationships among the original data. It can avoid the non-objectivity and bias of the measurement results to a certain extent and has obvious advantages in calculating the level of green development [32]. The calculation steps are as follows:

Step 1: Data standardization processing

Because of the need to calculate the logarithm of the standardized index values in the calculation process of the entropy weighting method, in this study the index standardization method of non-zero transformation was adopted for data standardization [33]. The calculation formula is as follows:

Positive indicators:

$$\mathbf{x}'\_{ij} = \frac{(\mathbf{x}\_{ij} - \mathbf{x}\_{j\text{min}})}{(\mathbf{x}\_{j\text{max}} - \mathbf{x}\_{j\text{min}})} \times \mathbf{0.99} + \mathbf{0.01} \tag{1}$$

Negative indicator:

$$\mathbf{x}'\_{ij} = \frac{(\mathbf{x}\_{j\text{max}} - \mathbf{x}\_{ij})}{(\mathbf{x}\_{j\text{max}} - \mathbf{x}\_{j\text{min}})} \times \mathbf{0.99} + \mathbf{0.01} \tag{2}$$

where *x* 0 *ij* is the standardized value of the *j*-th indicator in the *i*-th year; *xij* is the original value of the *j*-th indicator in the *i*-th year; *xj*min is the minimum value of the *j*-th indicator; and *xj*max is the maximum value of the *j*-th indicator.

Step 2: Calculation of the specific gravity, as follows:

$$Y\_{ij} = \frac{\mathbf{x}'\_{ij}}{\sum\_{i=1}^{m} \mathbf{x}'\_{ij}} \tag{3}$$

where *Yij* is the proportion of the *j*-th indicator in the *i*-th year, and *m* is the number of years of the study.

Step 3: Calculation of the information entropy, as follows:

$$e\_{\dot{j}} = \frac{1}{-\ln m} \sum\_{i=1}^{m} \left( Y\_{\dot{i}\dot{j}} \times \ln Y\_{\dot{i}\dot{j}} \right) \tag{4}$$

where *e<sup>j</sup>* is the information entropy of the *j*-th index.

Step 4: Calculation of the information utility value, as follows:

$$d\_{\rangle} = 1 - e\_{\rangle} \tag{5}$$

where *w<sup>j</sup>* is the weight of the *j*-th indicator, and *n* is the number of indicators.

Step 5: Calculation of the level of green development

The subsystem of the green development level was calculated as follows:

$$w\_j = \frac{d\_j}{\sum\_{j=1}^n d\_j} \tag{6}$$

where *w<sup>j</sup>* is the weight of the *j*-th indicator, and *n* is the number of indicators.

Step 6: Calculation of the level of green development

The subsystem of the green development was calculated as follows:

$$GDL = \sum\_{j=1}^{n} \mathbf{x}\_{ij}' w\_j \tag{7}$$

where *GDL* is the overall green development level of the system, with a value range of (0,1). The larger the value, the higher the overall green development level of the system, and vice versa.

#### *2.4. System Dynamics Model Theory*

System dynamics (SD) is based on the feedback control theory. Through the qualitative analysis of the research problem, real data and information are used to quantitatively deduce the coupling relationship between variables, while computer simulation technology is employed to study the behavior of dynamic social systems [34], ensuring the integrated application of natural and social sciences [35]. The establishment of an SD model mainly includes the drawing of a causal loop diagram, a system flow chart, a model equation, and parameter alignment; the stability of the system is ensured through the historical test and the sensitivity test [36].

The commonly used software for modeling SD includes Vensim, Ithink, DYNAMO, and Stella [37]. Vensim is easy to operate, has a friendly interface, uses arrows to connect variables when building an SD model, and the relationship between the variables is written in the form of model equations, which allows users to easily modify the content of the model. Therefore, the Vensim PLE software was selected in this study to model the SD of the green development level of the Chengdu-Chongqing City Group in the Yangtze River economic belt.

#### Model Variables and Equations

The SD model mainly includes four types of variables: level variables, rate variables, auxiliary variables, and constant variables. The state variables and the rate variables are the key variables [38].

The four main types of model equations for SD are the state equation, the rate equation, the auxiliary equation, and the constant equation. A table function can be used when the nonlinear relationships between the variables are difficult to present through mathematical operations [39]. The following equations represent the basic mathematical expressions used by Vensim.

1. State equation:

$$LvS(t) = St\_0 + \int\_{t\_0}^{t} rateS(t)dt\tag{8}$$

where *LvS*(*t*): is the value of the state variable at time *t*; *St*<sup>0</sup> is the value of the state variable at initial time *t0*; and *rateS*(*t*) is the rate of change of the state variable. The state equation can also be written in a discrete form as follows:

$$\text{LEVEL.K} = \text{LEVEL.J} + \text{DT} \times (\text{INFLION.JK} - \text{OLITFLOAT} \text{.JK}) \tag{9}$$

where *LEVEL*.*K* is the value of the state variable at time *K*; *LEVEL*.*J* is the value of the state variable value at time *J*; *INFLOW* is the input rate; *OUTFLOW* is the output rate; and *DT* is the temporal interval from time *J* to time *K*.

2. Rate equation: where *rateS t*( ) : is the rate of change of the state variable; *LvS t*( ) : is the value of the state

2. Rate equation:

$$rateS(t) = g[LvS(t), aux(t), const] \tag{10}$$

*rateS t g LvS t aux t const* () () () = , , (10)

where *rateS*(*t*): is the rate of change of the state variable; *LvS*(*t*): is the value of the state variable at time *t*; and *aux*(*t*) is the value of the auxiliary variable at time *t*. 3. Table Functions: The table function is a function that describes the relationship between the independ-

where *LEVEL K*. is the value of the state variable at time *K*; *LEVEL J*. is the value of the state variable value at time *J*; *INFLOW* is the input rate; *OUTFLOW* is the output

#### 3. Table Functions: ent and dependent variables by means of a list. The table function can be edited in Vensim

The table function is a function that describes the relationship between the independent and dependent variables by means of a list. The table function can be edited in Vensim PLE software using the "With lookup" function; its mathematical description is as follows: PLE software using the "With lookup" function; its mathematical description is as follows: ( ) ( ) ( ) ( )( ) ( ) min max min max 1 1 2 2 , , , ,, ,, , *n n lookupname X X Y Y X Y X Y X Y* − (11)

$$\text{lookup} \left( \left[ (X\_{\text{min}}, X\_{\text{max}}) - (Y\_{\text{min}}, Y\_{\text{max}}) \right] / (X\_1, Y\_1) / (X\_2, Y\_2) , \dots (X\_{\text{fl}}, Y\_{\text{fl}}) \right) \tag{11}$$

where (*X*1,*Y*1),(*X*2,*Y*2), . . .(*Xn*,*Yn*) are the values of the independent and dependent variables that have been assigned. iables that have been assigned. *2.5. Shared Socio-economic Pathways (SSPs)* 

#### *2.5. Shared Socio-Economic Pathways (SSPs)* 2.5.1. Overview of Shared Socio-Economic Pathways (SSPs)

2.5.1. Overview of Shared Socio-Economic Pathways (SSPs) SSPs are a powerful tool to describe global socio-economic development scenarios.

*Water* **2022**, *14*, 3947 7 of 26

rate; and *DT* is the temporal interval from time *J* to time *K*.

SSPs are a powerful tool to describe global socio-economic development scenarios. This tool was launched in 2010 by the IPCC on the basis of Representative Concentration Pathways (RCP)s, taking into account factors such as population growth, economic development, environmental conditions, and technological progress [40–42]. They describe the possible future development of society without the impact of climate change or climate policy [43]. This tool was launched in 2010 by the IPCC on the basis of Representative Concentration Pathways (RCP)s, taking into account factors such as population growth, economic development, environmental conditions, and technological progress [40–42]. They describe the possible future development of society without the impact of climate change or climate policy [43]. Through a comprehensive analysis of existing global development frameworks, the

Through a comprehensive analysis of existing global development frameworks, the IPCC has identified five typical SSPs from the perspectives of both climate challenge adaptation and climate challenge mitigation, namely SSP1 (sustainable development pathway), SSP2 (intermediate pathway), SSP3 (regional competitive pathway), SSP4 (uneven pathway), and SSP5 (fossil fuel-based development pathway), as shown in Figure 2. The mitigation challenge in Figure 2 refers to slowing down the pace of climate change by multiple means, including carbon dioxide emissions reduction; the adaptation challenge refers to the adoption of measures to adapt to the various impacts of climate change and avoid catastrophic consequences. IPCC has identified five typical SSPs from the perspectives of both climate challenge adaptation and climate challenge mitigation, namely SSP1 (sustainable development pathway), SSP2 (intermediate pathway), SSP3 (regional competitive pathway), SSP4 (uneven pathway), and SSP5 (fossil fuel-based development pathway), as shown in Figure 2. The mitigation challenge in Figure 2 refers to slowing down the pace of climate change by multiple means, including carbon dioxide emissions reduction; the adaptation challenge refers to the adoption of measures to adapt to the various impacts of climate change and avoid catastrophic consequences.

**Figure 2. Figure 2.** Shared social economy path diagram. Shared social economy path diagram.

2.5.2. Simulation Scheme of the Future Green Development Level of the Chengdu-Chongqing City Group under the SSPs

The green development of each city in the Chengdu-Chongqing City Group entails a process of continuous change; accordingly, the indicators that characterize the level of urban green development are also continuously changing. In order to ensure the validity of the SD model, it was necessary to carry out a historical statistical description of the model parameters of each city in the Chengdu-Chongqing City Group [44]. In this study, the maximum, minimum, and average values of the six rate variables in each city from 2004 to 2018 were assessed to provide a reference interval for the setting of the model parameters for each city. The results for Chengdu and Chongqing are detailed in Appendix A.

The SSPs introduced by the IPCC contain almost all the possible future development trends. A research team led by Jiang Tong conducted a study on the future population and economy of China under the five SSPs, based on domestic and international population and economic forecasting studies [45–49]. At present, the research of Jiang Tong's team on China's population and economic future prediction has reached a temporal scale up to 2100 and has proved to be very mature [45,50]. Therefore, in this study, when setting the parameters of the model for the prediction of the future green development level of each city in the Chengdu-Chongqing City Group, the population growth rate and the GDP growth rate referred to the research results of Jiang Tong's team. The rate of the increase in the total energy consumption, arable land area, public green space area, and total water supply were assumed to be either high, medium-high, medium, or low according to the characteristics of the five SSPs considered. The details are shown in Table 1.

**Table 1.** Hypothetical scenarios of model parameters under the shared social economy path.


In synthesis, the following principles were followed to formulate the future simulation plan of each city in the Chengdu-Chongqing City Group based on the five SSPs. The model parameters of each city in the Chengdu-Chongqing City Group under the SSP2 path were mainly based on the development trend of the statistical data in the historical periods and cities. Moreover, the values of the model parameters under the SSP1, SSP3, SSP4, and SSP5 paths were either higher or lower than a certain level of the SSP2 path, according to the parameter assumptions and historical data fluctuations.

#### **3. Results and Discussion**

#### *3.1. Construction of the Green Development Level Indicator System and Indicator Weights*

Based on the understanding of the connotation of green development, this study constructed an indicator system to evaluate the level of green development of the urban agglomerations by referring to the evaluation index system used in previous studies on green development and the 56 green development indicators suggested by the National Development and Reform Commission of China.

The developed evaluation index system consists of three levels: the target level, which is the green development level of the urban agglomerations; the guideline level, which includes the socio-economic, resource, and environment levels; and the quality-of-life level. It is composed of 19 indicators at the indicator level. The results of the calculation of the indicator weights are shown in the Table 2.


**Table 2.** Evaluation indicator system of green development level of the Chengdu-Chongqing City Group.

w1: weight of the indicator layer relative to the target layer; w2: weight of the indicator layer relative to the guideline layer.

### *3.2. Green Development Level SD Model Construction*

### 3.2.1. Boundary Definition and Systematic Analysis

The SD model constructed in this study considered the administrative boundaries of 16 cities in the Chengdu-Chongqing City Group as the spatial boundaries. Considering China's national construction goals and the reliability of the historical data, the temporal limit of the system was set as covering the period 2004–2050; 2004–2018 was set as the historical data period, while 2019–2050 was set as the model simulation period, and the time step was set as one year.

The systematic analysis of the green development level of the Chengdu-Chongqing City Group allowed the identification of the key factors affecting the green development level. The complex green development level system was divided into four subsystems: society, economy, resources, and environment. The indicator types of each subsystem are shown in Figure 3.

After the analysis of the four subsystems, an analysis framework of the subsystem of the green development level of the urban agglomeration was constructed, and the mutual influence relationship among the subsystems was clarified, as shown in Figure 4.

*Water* **2022**, *14*, 3947 10 of 26

**Figure 3.** Systematic analysis of SD model of the green development level in the Chengdu-Chongqing City Group. **Figure 3.** Systematic analysis of SD model of the green development level in the Chengdu-Chongqing City Group. the green development level of the urban agglomeration was constructed, and the mutual influence relationship among the subsystems was clarified, as shown in Figure 4.

**Figure 4.** Analytical framework of SD model of the green development level in the Chengdu-Chongqing City Group. **Figure 4.** Analytical framework of SD model of the green development level in the Chengdu-Chongqing City Group.

#### 3.2.2. Causal Feedback Relation

**Figure 4.** Analytical framework of SD model of the green development level in the Chengdu-Chongqing City Group. 3.2.2. Causal Feedback Relation A system can be divided into an open system and a feedback system [51]. The SD model of the green development level of the Chengdu-Chongqing City Group belongs to the feedback system; accordingly, the feedback loop of this system and the causal rela-3.2.2. Causal Feedback Relation A system can be divided into an open system and a feedback system [51]. The SD model of the green development level of the Chengdu-Chongqing City Group belongs to the feedback system; accordingly, the feedback loop of this system and the causal relationship between the variables could be clarified through the causality diagram. Combining the model index analysis and the green development level system analysis framework, the causal relationship between the subsystems and the variables was assessed, and the causal loop diagram of the SD model of the green development level of the Chengdu-A system can be divided into an open system and a feedback system [51]. The SD model of the green development level of the Chengdu-Chongqing City Group belongs to the feedback system; accordingly, the feedback loop of this system and the causal relationship between the variables could be clarified through the causality diagram. Combining the model index analysis and the green development level system analysis framework, the causal relationship between the subsystems and the variables was assessed, and the causal loop diagram of the SD model of the green development level of the Chengdu-Chongqing City Group in the Yangtze economic belt was drawn using the SD software Vensim PLE, as shown in Figure 5.

tionship between the variables could be clarified through the causality diagram. Combining the model index analysis and the green development level system analysis framework, the causal relationship between the subsystems and the variables was assessed, and the

Chongqing City Group in the Yangtze economic belt was drawn using the SD software

Chongqing City Group in the Yangtze economic belt was drawn using the SD software

Vensim PLE, as shown in Figure 5.

Vensim PLE, as shown in Figure 5.

*Water* **2022**, *14*, 3947 11 of 26

**Figure 5.** The causal loop of the SD model of green development level in the Chengdu-Chongqing City Group. **Figure 5.** The causal loop of the SD model of green development level in the Chengdu-Chongqing City Group. City Group.

#### 3.2.3. System Flow Chart 3.2.3. System Flow Chart 3.2.3. System Flow Chart The system flow chart is the core of an SD model. Based on the causal loop diagram

The system flow chart is the core of an SD model. Based on the causal loop diagram of the green development level of the Chengdu-Chongqing City Group, the system flow chart was further drawn. It includes 48 variables, including 6 state variables (i.e., total population, total GDP, total energy consumption, cultivated land area, public green area, and total water supply), 6 rate variables (i.e., population growth rate, GDP growth rate, increase rate of total energy consumption, increase rate of cultivated land area, increase rate of public green space, and increase rate of total water supply), and 36 auxiliary variables, as shown in Figure 6. The system flow chart is the core of an SD model. Based on the causal loop diagram of the green development level of the Chengdu-Chongqing City Group, the system flow chart was further drawn. It includes 48 variables, including 6 state variables (i.e., total population, total GDP, total energy consumption, cultivated land area, public green area, and total water supply), 6 rate variables (i.e., population growth rate, GDP growth rate, increase rate of total energy consumption, increase rate of cultivated land area, increase rate of public green space, and increase rate of total water supply), and 36 auxiliary variables, as shown in Figure 6. of the green development level of the Chengdu-Chongqing City Group, the system flow chart was further drawn. It includes 48 variables, including 6 state variables (i.e., total population, total GDP, total energy consumption, cultivated land area, public green area, and total water supply), 6 rate variables (i.e., population growth rate, GDP growth rate, increase rate of total energy consumption, increase rate of cultivated land area, increase rate of public green space, and increase rate of total water supply), and 36 auxiliary variables, as shown in Figure 6.

**Figure 6.** System flow chart of SD model of the green development level in the Chengdu-Chongqing City Group.

#### 3.2.4. Establishment of Model Parameters and Equations

The variables of the SD model were divided into four main types: state variables, rate variables, auxiliary variables, and constants. The variables used in this model and their properties are shown in Table 3.

**Table 3.** SD model variables of the green development level in the Chengdu-Chongqing City Group.


L: state variables; R: rate variables; A: auxiliary variables; C: constants.

The model equations were established in the following three ways: for variables with clear mathematical relationships, they were established directly; for variables with random distributions and no clear mathematical relationships, the dataset was fitted using fitted regression analysis; for variables without functional relationships, or with functional relationships that were or difficult to establish, the table functions were created with the help of the "With lookup" function in the Vensim software. The model equations are shown in Table 4.


**Table 4.** Model equations for explicit mathematical relationships.

The parameters in the SD model generally include the constant values and the initial values of the state variables in the model. The constant values in the SD model of the green development level of the Chengdu-Chongqing City Group refer to six rate variables. Considering that the various cities were at different stages of development during the period of 2004–2018, the rate variables were assigned fixed values that were not in line with the actual situation. Therefore, a time-continuous function was adopted in this study to represent these variables; the specific functional relationship is shown in Table 5. The model also includes six state variables, whose initial values were taken from the statistical data of each city in 2003, as detailed in Appendix A.


**Table 5.** Chengdu and Chongqing rate variable time-continuous functions.

#### *3.3. Model Validity Test*

#### 3.3.1. Historical Test

The historical test aims to verify whether the relative error between the model simulation results and the historical values is within an acceptable range [52]; the formula employed in this study is as follows:

$$D\_t = \frac{X\_t' - X\_t}{X\_t} \times 100\% \tag{12}$$

where *D<sup>t</sup>* is the relative error value at time *t*; *X* 0 *t* is the model simulation value at time *t*; and *Xt* is the historically true value at time *t*.

Six core indicators in the SD model of the green development level of the Chengdu-Chongqing City Group were selected to perform the historical test: total population, total GDP, total energy consumption, total water supply, arable land area, and public green space area. The passing rate of the historical test of each city in the Chengdu-Chongqing City Group is shown in Appendix A. The results of the historical test show that the SD model established in this study can fit the historical data well.

#### 3.3.2. Sensitivity Check

In a stable SD model, the change in its simulation values when the model parameters change is relatively low, i.e., the SD model has low sensitivity [8,53]. To verify whether the SD model constructed in this study could operate effectively, six rate variables were selected to test its sensitivity. The sensitivity of the rate variables of each city in the Chengdu-Chongqing City Group from 2004 to 2018 was derived based on the sensitivity formula with a 10% increase and a 10% decrease. The sensitivity check formula employed is as follows:

$$\mathcal{S}\_{Lij} = \left| \frac{\Delta L\_{it}}{L\_{it}} \times \frac{\mathbf{X}\_{jt}}{\Delta \mathbf{X}\_{jt}} \right| \tag{13}$$

$$\mathcal{S}\_{\mathbf{X}j} = \frac{1}{\mathcal{N}} \sum\_{i=1}^{N} \mathcal{S}\_{Lij} \tag{14}$$

where *Lit* is the value of the *i*-th state variable L at time *t*; ∆*Lit* is the change in the *i*-th state variable L at time t; *Xjt* is the value of the *j*-th parameter X at time t; ∆*Xjt* is the change in the *j*-th parameter X at time t; *SLij* is the sensitivity of the *i*-th state variable L to the *j*-th parameter X; and *SXj* is the sensitivity of the model to the *j*-th parameter X.

The results of the sensitivity test showed that the average sensitivity of the six rate variables to the SD model was lower than 0.01 in the two cases of the 10% increase and the 10% decrease. Hence, the structure of the SD model of the green development level of the Chengdu-Chongqing City Group constructed in this study was stable and could be used for the purposes of simulation and prediction.
