Scenario Modeling of Urbanization Development and Water Scarcity Based on System Dynamics: A Case Study of Beijing–Tianjin–Hebei Urban Agglomeration, China

Due to the accelerated process of urbanization in China, urban agglomerations have become the core areas for human settlement and economic development. High population and economic density has brought great pressure on water supply. Water scarcity is increasingly becoming one of the most important issues for the sustainable and healthy development of China’s urban agglomerations. In this paper, a system dynamics model was constructed to simulate the current conditions and future scenarios of urbanization development and water scarcity in the Beijing–Tianjin–Hebei (BTH) urban agglomeration in 2000–2030, by examining the interaction and feedback between the six major subsystems: water supply, water demand, water pollution, population urbanization, economic urbanization, and land urbanization. It is found that the South-to-North Water Diversion Project and the improved Reclaimed Water Reuse System may greatly increase the water supply. However, the speed of population urbanization and economic growth, the spatial structure of urban agglomeration and the water consumption pattern may determine the water demand. Although all scenarios may risk water scarcity in the future at some point, we could detect a comprehensive and relatively rational scenario to balance water scarcity, regional equity, and efficiency. It might help to synthetically understand the coordinated development mode between urbanization and water resources in Beijing–Tianjin–Hebei (BTH) urban agglomeration, and provide a useful analytical and decision support tool for scientists and policy-makers to achieve the sustainable urbanization development and water resource management.


Introduction
Water is a fundamental and indispensable resource for human survival and ecosystem evolution [1][2][3]. It also plays a supporting role in the sustainable development of the socio-economic system [4][5][6]. Due to rapid urbanization and industrialization, urban areas have undergone steady population growth and economic development. These areas face serious water scarcity and pollution both in developed and developing countries, which are the major threats to public health [7][8][9]. Specifically, urbanization could lead to water consumption increase in urban areas. It usually means more wastewater effluent and less water in rivers and lakes, both resulting in water quality

Study Area and Data Sources
BTH urban agglomeration, one of the world's largest agglomerations, is located in the Haihe River Basin and the North China Plain. It has 13 cities at prefecture level and above, including Beijing, Tianjin, Shijiazhuang, Tangshan, Qinhuangdao, Handan, Xingtai, Baoding, Zhangjiakou, Chengde, Cangzhou, Langfang, and Hengshui ( Figure 1). It covers an area of 0.2168 million km 2 , which accounts for 2.26% of China. It belongs to a typical continental monsoon climate type, which is characterized by hot and rainy summers, and cold and dry winters. As a semi-humid and semi-arid region, the mean annual precipitation varies from 400 to 800 mm, and the mean annual potential evaporation ranges from 1100 to 2000 mm [45,46]. According to statistics, in 2014, it has a population of 112 million, which accounts for 8.18% of China. The urbanization rate, which is the ratio of urban population to total population, is 61.0%. The gross domestic product (GDP) is 6649 billion yuan at current price, which accounts for 10.45% of China. It is the political, economic and cultural center of China, as well as one of the most competitive support platforms for China's international economic system [47][48][49]. However, the gross amount of water resources is 20.369 billion m 3 , only 0.75% of China. Its per capita water resources are 182 m 3 , only 9.14% of the national average and much lower than 500 m 3 , meaning an extremely severe water scarcity [12]. The total water consumption is 25 billion m 3 , accounting for 4.10% of China. Its per capita water consumption is 224 m 3 , which is 50% of the national average. Water consumption per ten thousand yuan GDP is 38 m 3 , which is 40% of the national average, meaning high water use efficiency [45]. In a word, BTH urban agglomeration is not only one of the areas owning the highest population density and the most vigorous economics, but also one of the most prominent regions of contradiction between human and water in China and even throughout the world [50,51]. The basic data in this paper include socio-economic data, water resources data, land use data and eco-environmental data for the 13 cities at prefecture level and above in 2000-2014. The The basic data in this paper include socio-economic data, water resources data, land use data and eco-environmental data for the 13 cities at prefecture level and above in 2000-2014. The socio-economic data, land use data, and eco-environmental data are all obtained from the past years statistical yearbooks of Beijing, Tianjin, and Hebei, respectively. The water resources data are all obtained in the past years water resources communique of Beijing, Tianjin, and Hebei respectively. To make economic data comparable in time series, all the economic data are calculated at comparable prices in 2000. For the total and average indexes of the whole region of BTH urban agglomeration, we firstly sum up the total indexes of the 13 cities at prefecture level and above, and then calculate the average indexes.

General Framework and Concept Model
A system dynamics model starts with the development of a dynamic hypothesis, generally referred to a causal loop diagram, which is a useful qualitative analytical tool for representing the relationships among system variables that produce dynamic feedback structure [32,41]. We use Vensim 5.11 software (Ventana Systems Inc., http://www.ventanasystems.com/) to draw the causal loop diagram of BTH urban agglomeration to examine the feedback processes between water resource and urbanization development ( Figure 2). socio-economic data, land use data, and eco-environmental data are all obtained from the past years statistical yearbooks of Beijing, Tianjin, and Hebei, respectively. The water resources data are all obtained in the past years water resources communique of Beijing, Tianjin, and Hebei respectively. To make economic data comparable in time series, all the economic data are calculated at comparable prices in 2000. For the total and average indexes of the whole region of BTH urban agglomeration, we firstly sum up the total indexes of the 13 cities at prefecture level and above, and then calculate the average indexes.

General Framework and Concept Model
A system dynamics model starts with the development of a dynamic hypothesis, generally referred to a causal loop diagram, which is a useful qualitative analytical tool for representing the relationships among system variables that produce dynamic feedback structure [32,41]. We use Vensim 5.11 software (Ventana Systems Inc., http://www.ventanasystems.com/) to draw the causal loop diagram of BTH urban agglomeration to examine the feedback processes between water resource and urbanization development ( Figure 2  As shown in Figure 2, the positive (+) or negative (−) impact of a variable on another variable is indicated by the solid arrow. The variables are directly connected to each other by links, generating positive and negative feedback loops. The positive feedback loops indicate a reinforcing change within the water resources and urbanization development system. The negative feedback loops indicate a depressed and balancing process. In this concept model, the main loops of cause and effect are as follows.
(1) A negative feedback loop between population and water demand: Total population →+ Domestic water demand →+ Total water demand →+ Water supply and demand gap →-Total population.
(2) A negative feedback loop between population and water environment: Total population→+ Urban domestic water demand →+ Wastewater effluent →+ Water environment pollution→-Total population. As shown in Figure 2, the positive (+) or negative (−) impact of a variable on another variable is indicated by the solid arrow. The variables are directly connected to each other by links, generating positive and negative feedback loops. The positive feedback loops indicate a reinforcing change within the water resources and urbanization development system. The negative feedback loops indicate a depressed and balancing process. In this concept model, the main loops of cause and effect are as follows.
(1) A negative feedback loop between population and water demand: Total population →+ Domestic water demand →+ Total water demand →+ Water supply and demand gap →-Total population.
(2) A negative feedback loop between population and water environment: Total population→+ Urban domestic water demand →+ Wastewater effluent →+ Water environment pollution→-Total population.
(3) A negative feedback loop between economy and water demand: GDP →+ Water demand for production →+ Total water demand →+ Water supply and demand gap →-GDP.
(4) A negative feedback loop between economy and water environment: GDP →+ Water demand for production →+ Total water demand →+ Wastewater effluent →+ Water environmental pollution →-GDP.
(5) A positive feedback loop between economy and water demand: GDP →+ Water conservancy investment →+ Total water supply →-Water supply and demand gap →+ GDP.
The above negative feedback loops indicate that population growth and economic development may increase water demand, but water scarcity and water pollution may restrain the population growth and economic development. The above positive feedback loops indicate that economic development may increase water supply by increasing investment in water conservancy, and lessen the constraints caused by water scarcity and water pollution. Besides the above-mentioned main feedback loops, there are other complex interactions between water resources and urbanization development system through its structure and efficiency changes ( Figure 2).

Dynamic Simulation Model Settings and Description
Based on the causal loop diagram which is conceptual and qualitative, we divide the water resources and urbanization development system into six major subsystems: water supply, water demand, water pollution, population urbanization, economic urbanization, and land urbanization. Subsequently, we use Vensim 5.11 software to draw a stock and flow diagram (Figure 3), which can quantificationally simulate system behavior for each city of BTH urban agglomeration over a 30-year period (2000-2030). We consider not only the interactions among the stocks (level variables) and flows (rate variables) of one city's water resources and urbanization development system, but also the impacts of the water resources and population flows among its neighboring cities. Through the scenario simulations of the 13 cities at prefecture level and above, we sum the total indexes for the whole region of BTH urban agglomeration, and reveal the effects of the spatial structure and the urban-rural structure changes on water scarcity. (3) A negative feedback loop between economy and water demand: GDP →+ Water demand for production →+ Total water demand →+ Water supply and demand gap →-GDP.
(4) A negative feedback loop between economy and water environment: GDP →+ Water demand for production →+ Total water demand →+ Wastewater effluent →+ Water environmental pollution →-GDP.
(5) A positive feedback loop between economy and water demand: GDP → + Water conservancy investment →+ Total water supply →-Water supply and demand gap→+ GDP.
The above negative feedback loops indicate that population growth and economic development may increase water demand, but water scarcity and water pollution may restrain the population growth and economic development. The above positive feedback loops indicate that economic development may increase water supply by increasing investment in water conservancy, and lessen the constraints caused by water scarcity and water pollution. Besides the above-mentioned main feedback loops, there are other complex interactions between water resources and urbanization development system through its structure and efficiency changes ( Figure 2).

Dynamic Simulation Model Settings and Description
Based on the causal loop diagram which is conceptual and qualitative, we divide the water resources and urbanization development system into six major subsystems: water supply, water demand, water pollution, population urbanization, economic urbanization, and land urbanization. Subsequently, we use Vensim 5.11 software to draw a stock and flow diagram (Figure 3), which can quantificationally simulate system behavior for each city of BTH urban agglomeration over a 30-year period (2000-2030). We consider not only the interactions among the stocks (level variables) and flows (rate variables) of one city's water resources and urbanization development system, but also the impacts of the water resources and population flows among its neighboring cities. Through the scenario simulations of the 13 cities at prefecture level and above, we sum the total indexes for the whole region of BTH urban agglomeration, and reveal the effects of the spatial structure and the urban-rural structure changes on water scarcity.

Water Supply Subsystem
The core variable of water supply subsystem is total water supply (TWS), which is composed of local water resources (LWR), transit water resources (TWR), utilization rate of water resources (RT), transfer water resources (TFWR), and unconventional water resources (UWR). The description function is denoted by Equation (1): where TFWR t is the water supply by the South-to-North Diversion Project at time t. UWR t contains the desalinated seawater (DSW) and reused wastewater (RWW) at time t.

Water Demand Subsystem
The core variable of water demand subsystem is total water demand (TWD), which is composed of agricultural water demand (AWD), industrial water demand (IWD), domestic water demand (DWD), and ecological water demand (EWD). The description function is denoted by Equation (2): Agricultural water demand (AWD) is composed of water demand for farmland (FWD), woodland (WWD), orchard (OWD), and grassland (GWD). Each kind of water demand is the product of land use area (A i ), the corresponding coefficient of irrigation (α i ), and the irrigation water-demand quota (β i ). Thus, the water demand subsystem is connected to the land urbanization subsystem.
Industrial water demand (IWD) is the product of industrial added value (IAV t ) and industrial water-demand quota, which is defined as water demand per unit of industrial added value (IWQ t ). Thus, the water demand subsystem is connected to the economic urbanization subsystem.
Domestic water demand (DWD) is composed of urban domestic water demand (UDWD) and rural domestic water demand (RDWD), which is the product of urban (UP) and rural population (RP) and the corresponding water-demand quota (UPQ and RPQ), respectively. Thus, the water demand subsystem is connected to the population urbanization subsystem.
Ecological water is mainly used for irrigation in urban green space and water discharge to rivers or lakes for sustaining the functionality of freshwater ecosystems. To achieve sustainable development, we calculate ecological water demand (EWD) as Equation (6): where t is the current time, t 0 is the initial time, and r e is the growth of EWD at time t.

Water Pollution Subsystem
The core variable of water pollution subsystem is wastewater effluent (WWE). It is composed of industrial wastewater effluent (IWWE) and domestic wastewater effluent (DWWE), which is the product of industrial (IWD) and domestic water demand (DWD) and the corresponding effluent coefficient (γ and δ), respectively. Thus, the water pollution subsystem is connected to the water demand subsystem.
As reused wastewater (RWW) is the product of wastewater effluent (WWE) and wastewater reuse rate (WWRR), the water pollution subsystem is also connected to the water supply subsystem.

Population Urbanization Subsystem
The core variables of population urbanization subsystem are total population (P), urban population (UP) and rural population (RP), which are direct or indirect factors that affect the water demand and supply as mentioned above. The description functions are mainly as follows: where P t is population in the current time; P t0 is population in the initial time; Birth(t) and Death(t) are the births and deaths at time t, respectively; InP(t) and OutP(t) are the immigration and emigration at time t, respectively; and µ t is population urbanization ratio.

Economic Urbanization Subsystem
The core variable of economic urbanization subsystem is gross domestic product (GDP), which not only determines the water demand, but also affects the water supply through water conservancy investment (WCI). It is composed of added value of the primary industry (AVPI), added value of the secondary industry (AVSI), and added value of the tertiary industry (AVTI).
The added value of the primary industry (AVPI) is the product of agricultural output value (AOV) and the corresponding coefficient (η). The AOV is composed of output value of planting (OVP), forestry (OVF), livestock (OVL), fishery (OVFY), and sideline (OVS), which is calculated by growth rate method and connected to agricultural water demand (AWD) and different kinds of land use area (A i ), respectively.
The added value of the secondary industry (AVSI) is the product of industrial added value (IAV) and the corresponding coefficient (θ). The IAV is calculated by growth rate method and connected to industrial water demand (IWD) as mentioned above.
The added value of the tertiary industry (AVTI) is calculated by growth rate method as Equation (15): where t is the current time, t 0 is the initial time, and RT is the growth of AVTI at time t.

Land Urbanization Subsystem
The core variables of land urbanization subsystem are farmland area (A farm ), woodland area (A wood ), orchard (A orch ), grassland (A grass ), and construction land area (A con ). The description functions are mainly as follows: where A it is the ith kind of land use area in the current time; A it0 is the ith kind of land use area in the initial time; and AInc i (t) and ADcr i (t) are the increase and decrease of the ith kind of land use area at time t, respectively. As mentioned above, they are connected to the water demand subsystem and the economic urbanization subsystem.

Parameters and Scenario Settings
The objective of the SD model is to provide a decision support tool for policy-makers to explore plausible policy scenarios and grasp the sustainable water resources management. Due to the complexity of the water resources and urbanization development system, it is difficult to determine how it works. To overcome this difficulty, for general parameters, we assume that their changes are consistent with their respective historical development trends during 2000-2014. However, for key parameters, we design some scenarios from three major aspects to simplify the simulations.

Scenario Design of Water Supply
Water supply is determined by complex natural factors and human activities [2,3]. Parameters such as transfer water resources (TFWR), the utilization rates of water resources (RT), the industrial and domestic wastewater effluent coefficient (γ and δ), the desalinated seawater (DSW), and wastewater reuse rate (WWRR), are modeled by numerical experiments according to the historical data during 2000-2014, or by table functions according to the Government's planning. For parameters such as local water resources (LWR) and transit water resources (TWR), we use their average annual values during 2000-2014 as a baseline. However, local water resources and transit water resources in BTH urban agglomeration have greatly decreased since 1956 [45,46]. Moreover, the actual water withdrawal of the South-to-North Water Diversion Project is usually lower than the planning due to engineering and financial situation. Therefore, we assume three scenarios for water supply: (1) High water supply scheme, in which local water resources (LWR) and transit water resources (TWR) are 100% of the average annual values during 2000-2014, and transfer water resources (TFWR) are 100% of the planning values.
(2) Medium water supply scheme, in which local water resources (LWR) and transit water resources (TWR) are 75% of the average annual values during 2000-2014, and transfer water resources (TFWR) are 75% of the planning values.
(3) Low-water supply scheme, in which local water resources (LWR) and transit water resources (TWR) are 50% of the average annual values during 2000-2014, and transfer water resources (TFWR) are 50% of the planning values.

Scenario Design of Water Consumption Mode
Water-demand quotas serve as the control variables of water consumption mode. They are affected by complex factors such as local policies, regulations, water prices, living standards, technologies, public awareness of water conservation, and water resources management [37,38,43]. To simulate them succinctly and rationally, we assume two kinds of water consumption mode: (1) Water-consuming mode, in which all the water-demand quotas are forecasted by the trend extrapolation of the historical data during 2000-2014, and the results are then used as a baseline.
(2) Water-saving mode, in which all the water-demand quotas are less than the baseline according to the possibility of water saving in each city. For instance, the irrigation water-demand quotas for farmland, woodland, orchard, and grassland are 5-10 m 3 /ha less than the baseline. The industrial water-demand quotas are 1-2 m 3 / ten thousand yuan less than the baseline. The urban and rural domestic water-demand quotas are 5-10 L per capita per day less than the baseline. On the whole, through the water-saving measures scenario and no water-saving measures scenario, we may detect the impacts of water consumption mode on the changes of water demand in BTH urban agglomeration.

Scenario Design of Urbanization Development
The socio-economic development system, including the population urbanization subsystem, economic urbanization subsystem, and land urbanization subsystem, is affected by various policies and uncertain factors. We may set numerous scenarios for it. To be concise and grasp the focal point, we assume three kinds of urbanization development modes: (1) Core development mode, in which BTH urban agglomeration continues the previous development mode. The growth rates of population and economy of all cities are forecasted by the trend extrapolation of the historical data during 2000-2014, and the results are then used as a baseline. Beijing is continuing to be the core city to gather more population and industries due to its agglomeration effects. It has relatively high growth rates of population and economy.
(2) Subcore development mode, in which the population and urban function of Beijing are relieved. Its population urbanization ratio is one percentage point lower than its baseline, and its growth rate of each industry is 0.5 percentage point lower than its baseline. Tianjin, Shijiazhuang, Baoding, Tangshan, Langfang, and Cangzhou are considered as subcores and key development areas. Their population urbanization ratios are one percentage point higher than the baseline, and their growth rates of each industry are 0.5 percentage point higher than the baseline, respectively.
(3) Multinode development mode, in which other cities such as Qinhuangdao, Handan, Xingtai, Hengshui, Chengde, and Zhangjiakou are considered as key development areas to pursue a balanced development in BTH urban agglomeration. Their population urbanization ratios are one percentage point higher than the baseline, and their growth rates of each industry are 0.5 percentage points higher than the baseline, respectively.

Validation Results of the SD Model
The validity of the SD model is the precondition for forecasting and analyzing the future scenarios of urbanization development and water scarcity in BTH urban agglomeration. It could be indicated by the errors between the simulated results and the existing historical data. If the errors between the simulated and true values are in the interval −10% to 10%, the results could be acceptable. However, if the errors are larger, the model must be modified. Therefore, we select some important variables to test the SD model, including total population (P), urban population (UP), gross domestic product (GDP), added value of the secondary industry (AVSI), added value of the tertiary industry (AVTI), farmland area (A farm ), construction land area (A con ), total water demand (TWD), agricultural water demand

Water Supply under Different Scenarios
As the South-to-North Water Diversion Project is taking effect step by step, water supply planning for Beijing, Tianjin, and Hebei will increase 12 × 10 8 m 3 , 15 × 10 8 m 3 , and 30 × 10 8 m 3 in 2030, respectively. Water supply planning for each city in the south-central of Hebei province will increasẽ 3 × 10 8 m 3 . Thus total water supply in BTH urban agglomeration under different scenarios will be increased by 28.5 × 10 8 m 3 to 57 × 10 8 m 3 . On the other hand, the wastewater reuse rate in BTH urban agglomeration will increase from 21%, in 2015, to~36%, in 2030. Specifically, Beijing will increase from 60% in 2015 to~70% in 2030. Tianjin will increase from 30% in 2015 to~50% in 2030. Most cities in Hebei province will reach 15% to 30% in 2030. Thus total reused wastewater in BTH urban agglomeration will be~50 × 10 8 m 3 to 90 × 10 8 m 3 . The South-to-North Water Diversion Project and the improved Reclaimed Water Reuse System will greatly increase water supply. On the whole, in the high water supply scheme, total water supply in BTH urban agglomeration will increase from 265.  Table 1.

Urbanization Development under Different Scenarios
In the core development mode, total population in BTH urban agglomeration will increase from 118.9 million in 2020 to 132.6 million in 2030. The urbanization rate will increase from 67.5% in 2020 to 75.0% in 2030. The gross domestic product (GDP) will increase from 7707 billion yuan in 2020 to 15,371 billion yuan in 2030. In the subcore development mode, total population in BTH urban agglomeration will increase from 118.9 million in 2020 to 132.8 million in 2030. The urbanization rate will increase from 68.6% in 2020 to 75.7% in 2030. The gross domestic product (GDP) will increase from 7769 billion yuan in 2020 to 15,986 billion yuan in 2030. In the multinode development mode, total population in BTH urban agglomeration will increase from 118.9 million in 2020 to 133.6 million in 2030. The urbanization rate will increase from 67.8% in 2020 to 74.2% in 2030. The gross domestic product (GDP) will increase from 7752 billion yuan in 2020 to 15,255 billion yuan in 2030. Among the three kinds of scenarios, the subcore development mode has the largest urbanization rate and economic growth, whereas the multinode development mode has the lowest urbanization rate and economic growth. On the whole, the subcore development mode strikes a balance between development efficiency and regional equity while the other two are both one-sided. Great difference can be seen from the 13 cities in BTH urban agglomeration (Table 2).

Water Scarcity under Different Scenarios
As shown in Figures 6 and 7, water scarcities in water-consuming and water-saving modes under different urbanization development scenarios in BTH urban agglomeration are also different. Almost all scenarios have risks of water shortage for the 13 cities as at least one pillar grows downward in Figures 6 and 7, which means that the water supply and demand gap is negative. However, the water shortage ratios and water use efficiencies are various under different scenarios. From them, we could choose a comprehensive and relatively rational scenario to balance water scarcity, regional equity, and efficiency.
under different urbanization development scenarios in BTH urban agglomeration are also different. Almost all scenarios have risks of water shortage for the 13 cities as at least one pillar grows downward in Figures 6 and 7, which means that the water supply and demand gap is negative. However, the water shortage ratios and water use efficiencies are various under different scenarios. From them, we could choose a comprehensive and relatively rational scenario to balance water scarcity, regional equity, and efficiency.

Scenarios of Water-consuming Mode
If the water-demand quotas develop in accordance with historical trends, total water demand in 2030 in BTH urban agglomeration will be 259.6 × 10 8 m 3 , 269.8 × 10 8 m 3 , and 269.4 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. It is the lowest in the core development mode, and they are about the same in the subcore and multinode development mode. Water use efficiency will be 592.2 yuan/m 3 , 592.4 yuan/m 3 , and 566.2 yuan/m 3 correspondingly. It is the lowest in the multinode development mode, and they are about the same in the core and subcore development mode. It is seen that the core and subcore development mode may save more water and ease water scarcity due to higher water use efficiency. However, they may increase regional disparities and injustice according to the distribution of GDP (Table 2).

Scenarios of Water-Consuming Mode
If the water-demand quotas develop in accordance with historical trends, total water demand in 2030 in BTH urban agglomeration will be 259.6 × 10 8 m 3 , 269.8 × 10 8 m 3 , and 269.4 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. It is the lowest in the core development mode, and they are about the same in the subcore and multinode development mode. Water use efficiency will be 592.2 yuan/m 3 , 592.4 yuan/m 3 , and 566.2 yuan/m 3 correspondingly. It is the lowest in the multinode development mode, and they are about the same in the core and subcore development mode. It is seen that the core and subcore development mode may save more water and ease water scarcity due to higher water use efficiency. However, they may increase regional disparities and injustice according to the distribution of GDP (Table 2).
Moreover, in the low-water supply scheme, the water supply and demand gap in 2030 in BTH urban agglomeration will be −14.5 × 10 8 m 3 , −24.8 × 10 8 m 3 , and −24.4 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. The pillars in most cities grow downward in Figure 6c. It indicates that if BTH urban agglomeration does not make full use of the South-to-North Water Diversion Project and the improved Reclaimed Water Reuse System; it may still face water scarcity in the future 15 years no matter what kind of development mode it takes.
Finally, even in the high and medium water supply scheme, though the water supply and demand gaps of the whole urban agglomeration are above zero, those of some cities are still negative and the absolute values are relatively large (Figure 6a,b). It indicates that water deficiency in local regions will exist for a long time in BTH urban agglomeration unless extreme measures are taken to limit the development of these regions. Therefore, to reduce the risk of water shortage, a more severe water-saving mode had better be considered.
It indicates that if BTH urban agglomeration does not make full use of the South-to-North Water Diversion Project and the improved Reclaimed Water Reuse System; it may still face water scarcity in the future 15 years no matter what kind of development mode it takes.
Finally, even in the high and medium water supply scheme, though the water supply and demand gaps of the whole urban agglomeration are above zero, those of some cities are still negative and the absolute values are relatively large (Figure 6a,b). It indicates that water deficiency in local regions will exist for a long time in BTH urban agglomeration unless extreme measures are taken to limit the development of these regions. Therefore, to reduce the risk of water shortage, a more severe water-saving mode had better be considered.

Scenarios of Water-saving Mode
If the water-demand quotas are reduced according to the possibility of water saving in each city, total water demand in 2030 in BTH urban agglomeration will be 244.7 × 10 8 m 3 , 251.6 × 10 8 m 3 , and 253.0 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. Compared with the water-consuming mode, it may save 14.8 × 10 8 m 3 , 18.3 × 10 8 m 3 , and 16.4 × 10 8 m 3 of water resources. Correspondingly, water use efficiency will be 628.0 yuan/m 3 , 635.5 yuan/m 3 , and 602.9 yuan/m 3 . It is seen that water-saving potential and water use efficiency are both the highest in the subcore development mode. In other words, the subcore development mode could give due consideration to efficiency and fairness, as well as make full use of its advantages. However, it is undesirable to focus only on one city or equally on all cities.

Scenarios of Water-Saving Mode
If the water-demand quotas are reduced according to the possibility of water saving in each city, total water demand in 2030 in BTH urban agglomeration will be 244.7 × 10 8 m 3 , 251.6 × 10 8 m 3 , and 253.0 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. Compared with the water-consuming mode, it may save 14.8 × 10 8 m 3 , 18.3 × 10 8 m 3 , and 16.4 × 10 8 m 3 of water resources. Correspondingly, water use efficiency will be 628.0 yuan/m 3 , 635.5 yuan/m 3 , and 602.9 yuan/m 3 . It is seen that water-saving potential and water use efficiency are both the highest in the subcore development mode. In other words, the subcore development mode could give due consideration to efficiency and fairness, as well as make full use of its advantages. However, it is undesirable to focus only on one city or equally on all cities.
Specifically, in the low water supply scheme, the water supply and demand gap in 2030 in the whole urban agglomeration will be 0.3 × 10 8 m 3 , −6.5 × 10 8 m 3 , and −8.0 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. The risks of water scarcity in the future 15 years may gradually decrease to a low level. Although the water supply and demand gaps in some of the 13 cities are still negative (Figure 7c), the water shortage ratios are mostly below 10%. It indicates that water-saving is also one of the important keys to reduce water scarcity for BTH urban agglomeration.
As long as BTH urban agglomeration adopted strict water-saving measures, the risks of water scarcity could be controlled at a safety level.
Besides, in the high water supply scheme, the water supply and demand gap in 2030 in the whole urban agglomeration will be 72.9 × 10 8 m 3 , 66.1 × 10 8 m 3 , and 64.6 × 10 8 m 3 in the core, subcore, and multinode development modes, respectively. In the medium water supply scheme, they will be 42.1 × 10 8 m 3 , 35.3 × 10 8 m 3 , and 33.8 × 10 8 m 3 , respectively. Water supplies of most cities are larger than water demands as the pillars grow upward in Figure 7a,b. It indicates that, if BTH urban agglomeration adopts strict water-saving measures and makes full use of various water resources, there will be surplus water resources to restore the degraded ecosystem.

Conclusions
This paper constructed a SD model to simulate the current conditions and future scenarios of urbanization development and water scarcity in BTH urban agglomeration in 2000-2030, considering the nexus of water supply, water demand, water pollution, population urbanization, economic urbanization, and land urbanization, which are all closely related to public health. Based on the analysis and comparison of various scenarios, a relatively rational scheme to balance water scarcity and urbanization development was explored. The following conclusions were obtained.
(1) In the future 15 years, total water supply in BTH urban agglomeration has great uncertainty due to climate change, decrease of upstream water inflow, increase of water transfer, and utilization of unconventional water resources. It ranges from 245.0 × 10 8 m 3 to 317.7 × 10 8 m 3 in 2030 under different scenarios. The South-to-North Water Diversion Project and the improved Reclaimed Water Reuse System may greatly increase the water supply. Therefore, it is necessary to make full use of them to guarantee a steady water supply.
(2) The speed of population urbanization and economic growth, the spatial structure of urban agglomeration and the water consumption pattern may determine the water demand. All scenarios may have risks of water scarcity for the 13 cities, especially in the low-water supply scheme. However, if the local governments adjust the urbanization development and water consumption mode, and adopt vigorous measures to support a high and medium water supply, the water shortage ratios of the 13 cities may gradually decrease to a safety level. Some surplus water resources may be used to restore the degraded ecosystem.
(3) Among different scenarios, water-saving potential and water use efficiency in the water-saving and subcore development modes are the highest. It may be chosen as a comprehensive and relatively rational scenario to balance water scarcity, regional equity, and efficiency. In this scenario, total population in BTH urban agglomeration will increase from 118. 9 million, in 2020, to 132.8 million, in 2030. The urbanization rate will increase from 68.6% in 2020 to 75.7% in 2030. GDP will increase from 7769 billion yuan in 2020 to 15,986 billion yuan in 2030. The whole region will be not short of water in the high and medium water supply scheme. Even in the low-water supply scheme, the water shortage ratio will decrease from 15.0% in 2020 to 2.6% in 2030. However, water-saving potential and water use efficiency in the water-consuming and multinode development modes are the lowest on the whole, suggesting that high water consumption and completely balanced development strategy should be avoided.