*Article* **A SEEC Model Based on the DPSIR Framework Approach for Watershed Ecological Security Risk Assessment: A Case Study in Northwest China**

**Bin Wang 1,2, Fang Yu 2,\*, Yanguo Teng 1,\*, Guozhi Cao 2, Dan Zhao <sup>2</sup> and Mingyan Zhao <sup>3</sup>**


**Abstract:** The DPSIR model is a conceptual model established by the European Environment Agency to solve environmental problems. It provides an overall framework for analysis of environmental problems from five aspects: driving force (D), pressure (P), state (S), impact (I), and response (R). Through use of the DPSIR model framework, this paper presents the SEEC model approach for evaluating watershed ecological security. The SEEC model considers four aspects: socioeconomic impact (S), ecological health (E), ecosystem services function (E), and control management (C). Through screening, 38 evaluation indicators of the SEEC model were determined. The evaluation results showed that the ecological security index of the study area was >80, indicating a generally safe level. The lowest score was mainly attributable to the low rate of treatment of rural domestic sewage. The water quality status was used to evaluate the applicability of the SEEC model, and the calculation results indicated that the higher the score of the ecological security evaluation results, the better the water quality status. The findings show that the SEEC model demonstrates satisfactory applicability to evaluation of watershed ecological security.

**Keywords:** watershed ecological security assessment; DPSIR model framework; environmental management

#### **1. Introduction**

The footprints of human activities have covered the world [1]. In the process of rapid development of both industry and agriculture, the ecological environment has suffered unprecedented damage [2,3]. Globally, the soil [4,5], water [6,7], air [8,9], and other environmental media in areas with frequent human activity are in a state of continuous deterioration [10,11]. Ecosystem degradation and environmental pollution are gradually threatening and destroying human socioeconomic progress, survival, and development [12,13]. In recent years, researchers have attempted to evaluate the consequences and degree of risk associated with current changes of the ecological environment but without reaching consensus [14,15].

China remains in the process of rapid economic growth and urbanization. However, various ecological and environmental problems continue to emerge, threatening to destroy China's sustainable development and affecting the living conditions of the population [16]. China has experienced ecological and environmental crises in the past and it now faces many new challenges regarding environmental protection. Therefore, President Xi proposed the idea of an ecological civilization, which means that China's model of development has changed from that of "grow first, clean up later" to one of sustainable

**Citation:** Wang, B.; Yu, F.; Teng, Y.; Cao, G.; Zhao, D.; Zhao, M. A SEEC Model Based on the DPSIR Framework Approach for Watershed Ecological Security Risk Assessment: A Case Study in Northwest China. *Water* **2022**, *14*, 106. https://doi.org/ 10.3390/w14010106

Academic Editor: Tammo Steenhuis

Received: 30 November 2021 Accepted: 31 December 2021 Published: 4 January 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/).

development [17,18]. At the policy level and in everyday life, the expectation is for a safer and cleaner living environment. In the past, researchers often used the concept of environmental risk to assess whether the environment of a region might pose a threat to human health. Specifically, such assessments can be used to evaluate whether the current ecological situation of an area in which humans survive continues to be safe and whether it can ensure the environmental needs of human life and development. Therefore, evaluation of ecological safety is vital.

The main objective of ecological security assessment is to determine the ecological status and ecological pressure faced by a region under normal human activities [19]. It was formally proposed by the International Institute for Applied Systems Analysis in 1989 [20]. Evaluation results can be expressed using the ecological security index (ESI). A high ESI value indicates that the ecological state of the evaluation receptor is able to not only ensure the needs of human survival and development but also resist the pressures brought by human development. Therefore, through scientific evaluation of the factors on which the ESI is based, policies can be formulated to improve the situation. This approach also makes the work of ecological protection more refined and targeted, which is of great importance considering China's current state of development and environmental protection [21,22].

In accordance with different evaluation objects, ecological security can be divided into water ecological security [23], land ecological security [24], coastal ecological security [25], and urban ecological security [26], and all these aspects of ecological security assessment have been widely studied and applied. Broadly, ecological security includes natural ecological security, economic ecological security, and social ecological security, which mainly refers to a state in which human life, health, and resources are not threatened in terms of the above aspects. In this study, we were interested in the ecological security of a watershed, which refers to the ecological state of the lakes, rivers, and other areas within a catchment, and its ability to resist the ecological pressure brought by human activities.

Currently, model evaluation methods are used in ecological security assessment, and the most commonly used evaluation models include the PSR model, DSR model, and DPSIR model [27]. In 1979, Rapport and Friend proposed a model framework for analyzing and describing the interaction between socioeconomic development and the ecological environment, which was further improved by the Organization for Economic Co-operation and Development and United Nations Environment Programme, forming the PSR model framework [28]. The basic connotation of the PSR model is that human activities exert pressure on the environment and its natural resources (P-pressure), which changes the state of both the environment and the quality of the natural resources (S-state), forcing human society to respond to these state changes through adoption of policies, decisions, or management measures that affect the environment, economy, and land (Rresponse) [29]. The PSR model is suitable for ecological security evaluation on a small spatial scale and with few influencing factors. However, because it simplifies the causal relationship between indicators, it ignores the complexity of the system, especially the driving force factors of ecological security [30]. To overcome this weakness, the United Nations Conference on Sustainable Development established the DSR framework in 1996, in which the driving force factors (D) refer to the regional socioeconomic objectives which represent the fundamental environmental pressure. The DSR model can better characterize the impact of the driving force factors on ecosystem evolution, but the definitions of the driving force factors and the response factors in the model were vague. Therefore, to improve the applicability of the DSR model, in 1999, the European Environment Agency officially adopted the DPSIR model (driving force–pressure–state–influence–response), in which influence refers to the impact of changes in environmental status on environmental receptors [31]. The model combines the characteristics of the PSR and DSR frameworks. It has the advantages of comprehensiveness, systematicness, and flexibility, and covers five assessment factors and constructs a causal network between them to reflect the impact of socioeconomic development and human activities on the system state and the human response to adverse impacts [32].

In the DPSIR model framework suitable for watershed ecological security assessment, factor D generally includes population, socioeconomic, and other indicators; factor P generally includes pollutant discharge and other indicators; factor S generally includes water quality status, sediment status, and other indicators; factor I generally includes water service function and other indicators; and factor R generally includes river protection policy, ecological restoration, and other indicators. Generally, the DPSIR model framework is a circular system, i.e., the driving force leads to pressure, then the pressure changes the state, and the change of state has a consequential impact, which promotes a response that leads to adjustment of the driving force [33]. In recent years, the DPSIR model has been widely accepted and used in the process of ecological security research because it can reveal causal relationships between ecology and human activities [33–40]. It provides a conceptual model for a research scheme for evaluation of human activities, resources, the ecology, and sustainable development [41,42] and is also applied to interdisciplinary research [36,43]. Through application of the DPSIR model framework, many studies have performed ecological security evaluation of lakes, rivers, land, and oceans, thereby providing a scientific basis for further expansion of the connotation and application of the DPSIR model framework [44,45].

Although the DPSIR model has been used widely in many fields, applicability of the evaluation method has been limited owing to inconsistent selection of indicators, poor analysis of the reasons for the selection of indicators, and unclear determination of the process of index weighting [27,46–48]. Additionally, in previous watershed ecological security assessments, factors D and P were usually evaluated using socioeconomic and other related indicators, and ecological indicators were ignored, which resulted in overestimation of the impact of policy and economic development and underestimation of the ability of the ecosystem to deal with the pressure (factor P). Furthermore, the existing evaluation method lacks verification of the evaluation results, thereby diminishing the reliability and guidance of the evaluation results [49,50]. To resolve the problems of poor applicability of the DPSIR model to watershed ecological security evaluation and lack of a verification method, this study adopted the following research methods. By identifying the key factors affecting the ecological security of a watershed, and through analyzing the DPSIR model framework, the SEEC model including the process of indicator selection and the determination of weights was established. A study area was selected for application of the new model, and a method of water quality evaluation and analysis was innovatively used to evaluate the applicability of the SEEC model. The study area, located in an arid area in Northwest China, comprised a watershed that represents an important water supply source for a large city. However, owing to the specific geographical location and the harsh natural environment of the watershed, research data were scant, and therefore a watershed ecological security assessment was not undertaken. The ESI of the study area was obtained, and the reasons for low ESI values were analyzed, on the basis of which, suggestions for improvement of the ESI of the watershed were proposed. The results could serve both as a reference for subsequent environmental planning and management and as a scientific basis for comprehensive pollution control and ecological environmental protection of the study area.

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

*2.1. Construction of the Evaluation Model*

#### 2.1.1. Identification of the Model

Currently, methods used for watershed ecological security assessment are not unified. By identifying the connotation of the model, this article established a method suitable for watershed ecological security assessment under the framework of the DPSIR model. However, the DPSIR model only provides an evaluation framework, and it does not offer methods for selecting and evaluating the applicability of indicators. Therefore, using the DPSIR model framework, this study identified the primary indicators of each factor and, in combination with consideration of the key issues of watershed ecological security assessment, the SEEC model was constructed. The essence of the SEEC model is that it is a

representation of the DPSIR model framework specifically suited to watershed ecological security assessment. Therefore, in building the SEEC model, the connotation of the DPSIR model should be identified first, and an index system suitable for watershed ecological security assessment should also be established.

The essence of the DPSIR model is to identify the main factors affecting ecological security under the influence of human activities. It needs to determine the ecological state under the action of these factors, identify the impact, select relevant indicators of the response, evaluate the state of ecological security in terms of the five aspects, and obtain the ESI. Therefore, the key to using the DPSIR model to evaluate watershed ecological security is to accurately identify representative indicators that can characterize watershed ecological security.

The driving force factors (D) in the DPSIR model mainly represent social development and economic growth, reflecting the trends of population change, socioeconomic activity, and industrial economic development. These factors represent potential causes of environmental change, and they are also the most primitive and important indicators of change of the water environment security system.

Pressure factors (P) refer to the pressure applied directly to the water ecosystem through the driving force (D). Similar to D, P is an external force that affects the development and change of water ecosystem security. In previous research on urban ecological security, D and P were regarded as two separate factors, because the driving force and pressure can directly affect urban ecological security [51]. However, it is difficult to observe and calculate the impact of driving forces on watershed ecological security, mainly because most areas of many watersheds are not located in urban built-up areas and there are few human activities around. In this case, it needs to redefine the meaning of P used for watershed ecological security assessment under the DPSIR framework. Therefore, P mainly reflects the pollution load in this article.

State factors (S) refer to the state of the water ecosystem under the influence of both D and P. Thus, S can be illustrated directly through the characteristics of environmental media such as water quality and sediment, which are indicators of ecological health.

Impact factors (I) refer to the impact of the state of the water environment ecosystem on the economy and the livelihood and health of the population, which is the inevitable result of the interaction of the first three factors (D, P, and S). The ecological state changes caused by the above factors are mainly reflected in changes of the watershed ecological service functions. Therefore, the impact factors mainly include watershed ecological service functions such as water resources supply and the cultural landscape.

Response factors (R) refer to the countermeasures taken by humans to improve or adapt to the ecological state, which reflect the process of human regulation and management. Therefore, R mainly includes supervision of the ecological environment, pollution control, and industrial adjustment.

Because the DPSIR model overestimates nonecological factors such as the economy and population change, it is considered that the nature of the ecosystem itself is a more direct factor of ecological security, and thus it should receive greater attention. Moreover, with improvement of the level of ecological management, government departments also consider ecological factors when determining economic and population objectives, resulting in gradual closing of the relationship between the driving force (D) and pressure (P) factors. Considering D and P as a comprehensive factor (i.e., S—socioeconomic impact) allows more detailed analysis of the impact of human social activities on ecology. At the same time, by choosing ecological health (E) as the characterizing factor of S, ecosystem services function (E) as the characterizing factor of I, and control management (C) as the characterizing factor of R, the SEEC model can be established. The framework of the SEEC model is shown in Figure 1.

**Figure 1.** Framework of the SEEC model.

2.1.2. Construction of the System of Indicators

After determining the four factors of the SEEC model, the representative indicators must be screened to obtain the evaluation results of the SEEC model.

(1) Socioeconomic impact

As mentioned in the introduction, the factors D and P in the original DPSIR model were generally related to the socioeconomic factors in the conventional sense [52,53]. However, when assessing the ecological security of the watershed, the particularity of this assessment objective needs to be considered, that is, for the ecological security of the watershed, the driving force and pressure are not very intuitive in the general sense. This paper needs to consider the driving force and pressure indicators that have a more significant impact on the ecological security of the watershed. Moreover, in the assessment of watershed ecological security, the applicability of each assessment factor should be considered in a balanced manner, and the original DPSIR model cannot be relied on completely, which will over amplify the impact of D and P. Therefore, socioeconomic impact factors mainly include the socioeconomic development and pollution load attributable to human activities in this article. Socioeconomic development mainly includes population, economic, and social indicators. Population indicators include the quantity, density, and natural growth rate of the population. Economic indicators are mainly used to determine the level of economic development and the intensity of economic activities. Therefore, economic indicators that can fully represent the structure and scale of the economy should be selected, e.g., GDP, per capita GDP, and total industrial and agricultural output value. Social indicators are relatively comprehensive and can be expressed by the urbanization rate.

Watershed pollution load is the primary way in which human activities affect water ecology. Generally, point source or nonpoint source discharge of pollutants can have serious impact on a watershed.

(2) Ecological health

Ecological health is reflected by the water quality and water ecology. In the process of watershed monitoring, water quality is commonly assessed using dissolved oxygen, total nitrogen, total phosphorus, the permanganate index, ammonia nitrogen, transparency, suspended solids, chlorophyll a, and heavy metal indicators.

In addition to the above indexes, water ecological indicators also include zooplankton, phytoplankton, and the benthic biomass. However, because it is very difficult to monitor the above indicators within many watersheds, the quality of sediment is often used to characterize the ecological status of water.

(3) Ecosystem services

The ecological services function of a watershed is mainly reflected in purifying the water quality, providing aquatic products, and providing cultural tourism services. Considering that aquaculture and fishing are no longer allowed in many watersheds, inclusion of this indicator was dependent on the specific situation of the watershed. The water purification function of a watershed is generally realized through natural shoreline filtration on both banks of the river. Cultural tourism services are also related to the geographical location of the watershed.

(4) Control management

Control management is mainly reflected in policies related to the economy, ecology, and environment, including ecological protection capital investment, industrial structure adjustment, and ecological supervision capacity building.

Through analysis of the connotation, scientificity, representativeness, and applicability of each indicator and following literature-based research, 34 evaluation indicators were determined. See Table A1 for the definition, reasons for selection, and calculation or acquisition method of each indicator.

#### 2.1.3. Data Processing

#### (1) Data standardization

All indicators must be standardized to facilitate ecological security assessment. The concept of a reference standard should be introduced, meaning the value of each evaluation index in the ideal state (conducive to ecological security) or at the average level in a large-scale region [54].

The process of standardization of the indicators is conducted as follows:

$$\text{Positive irradiance: } R\_{i\bar{j}} = X\_{i\bar{j}} / S\_{i\bar{j}\prime} \tag{1}$$

$$\text{Negative inducator: } R\_{\vec{i}\vec{j}} = S\_{\vec{i}\vec{j}} / X\_{\vec{i}\vec{j}}.\tag{2}$$

where positive indicator means that the larger the value of the indicator, the more favorable it is to ecological security; negative indicator means that the larger the value of the indicator, the more unfavorable it is to ecological security; *Xij* is the measured value of indicator *i* at sampling point *j*; *Sij* is the reference standard of indicator *i*; and *Rij* is the dimensionless value of the evaluation indicator, where 0 < *Rij* < 1; when *Rij* > 1, we take *R*ij = 1.

(2) Weight calculation

The main methods used to determine the weights are the subjective weight method and the objective weight method. The most common subjective weighting method is the Delphi method, also known as the expert scoring method [55]. Its advantages are its clear concept, simplicity, and ease of operation, which can grasp the main factors of ecological security assessment, but it needs a certain number of experts with experience to produce the scores. The objective weighting method determines the index weights using a judgment matrix composed of evaluation index values. The most commonly used objective weighting method is the entropy method, which uses the utility values of the index information in the calculation; the higher the utility value, the more important it is to the evaluation. The SEEC model involves four factors, each of which contains information with differing complexity. Thus, the objective weighting method will lead to some factors with less impact on ecological security obtaining higher weighting, making it difficult to truly reflect the importance of the factors. Therefore, the weights of the four factors were determined using the expert scoring method. The full score was set at 10 and the higher the score, the more important the factor. Using a judgment matrix after sorting the consultation results, the weighting coefficient of the index was calculated, and the entropy method was used to determine the weights of the 34 indicators.

The process of application of the entropy method is as follows [56]:

(a) Construct the judgment matrix *Z* for *n* samples and *m* evaluation indicators:

$$Z = \begin{bmatrix} X\_{11} & X\_{12} & \dots & X\_{1m} \\ X\_{21} & X\_{22} & \dots & X\_{2m} \\ \vdots & \dots & \dots & \dots & \dots \\ X\_{n1} & X\_{n2} & \dots & X\_{nm} \end{bmatrix} \tag{3}$$

(b) The dimensionless data are used to obtain a new judgment matrix, in which the expression of the element is:

$$R = (r\_{\vec{i}\vec{j}} \, n \times m) \tag{4}$$

(c) According to the definition of entropy, for n samples and m evaluation indicators, the entropy of the evaluation indicators can be determined as follows:

$$H\_{\bar{i}} = -\frac{1}{1\,\mathrm{n}(n)} \left[ \sum\_{\bar{i}=1}^{n} f\_{\bar{i}\bar{j}} \,\mathrm{ln} f\_{\bar{i}\bar{j}} \right] \tag{5}$$

$$f\_{ij} = \frac{r\_{ij}}{\sum\_{i=1}^{n} r\_{ij}} \tag{6}$$

where 0 ≤ *Hi* ≤ 1.

To make ln*fij* meaningful, it is assumed that *fij* = 0, *fij*ln*fij* = 0, *i* = 1,2, ... ,*m*, and *j* = 1,2, . . . ,*n*.

(d) Calculate the entropy weight (*Wi*) of the evaluation indicators:

$$\mathcal{W}\_{\bar{i}} = \frac{1 - H\_{\bar{i}}}{m - \Sigma H\_{\bar{i}}} \tag{7}$$

where *Wi* is the weighting coefficients of the evaluation indicators that meet the following requirement:

$$
\Sigma \mathcal{W}\_{\text{l}} = 1 \tag{8}
$$

(3) Expression of the evaluation results

The evaluation results are expressed using the ESI; the higher the ESI value, the safer the ecological status. On the basis of the general expression of river and lake ecosystem assessment results in China, the ESI is divided into five levels from 0–100, as shown in Table 1.

**Table 1.** Classification standard of the ESI.


#### *2.2. Study Area*

The study area is located in Northwest China. The watershed is 214 km long and 25–50 km wide from east to west. The drainage area is 4684 km<sup>2</sup> and the annual runoff is 237 million m3. The average elevation is 2545.63 m (Figure 2).

**Figure 2.** Overview of the study area.

The study area is affected mainly by the mid-latitude near-surface atmospheric circulation. The annual average temperature of the watershed is 3.5 °C. The temporal distribution of precipitation is uneven, falling mainly in June–August. The annual average precipitation of the watershed is approximately 420 mm, and its spatial distribution is very uneven. Solid precipitation accounts for approximately 54.2% of the annual total precipitation. According to the results of the "Water Resources Bulletin," the average annual precipitation in the study area was 1.28 billion m<sup>3</sup> during 2006–2013, with no obvious trend of increase or decrease. In 2015, the amount of surface water resources within the study area was 1.160 billion m3, and the amount of groundwater resources was 574 million m3. The total amount of the water supply was maintained at more than 1.1 billion m3, and more than half of the water supply was derived from surface water.

The soil distribution in the upper reaches of the study area has vertical zonation. In addition to ice and snow cover, at elevations above 3600 m, exposed rocks and stone mounds are widely distributed, although some areas have poorly developed soil. The soil in the elevation range of 3100–3400 m is mainly alpine meadow soil; 2000–3100 m is mainly subalpine meadow soil; 1700–2000 m is mountain chernozem; 1200–1700 m is mountain chestnut soil; and 800–1200 m is brown calcareous soil. The soil at 1700–2900 m has mosaic distribution characteristics. Taupe forest soil is mainly found on shady slopes at elevation of 1700–2900 m, and it is distributed in a compound area with chernozem and subalpine meadow soil. The climate and terrain in the watershed vary markedly, and the corresponding vegetation types are relatively complete with obvious regularity in terms of geographical distribution. The main vegetation types are coniferous forest, broadleaved forest, shrub, grassland, meadow, and desert grassland. We used ArcGIS 10.3 to interpret the land use of the study area in 2014, and the results showed that grassland was the main vegetation cover type, followed by woodland and cultivated land (Figure 3).

**Figure 3.** Land use types within the study area (2014).

In addition to the main stream, the study area includes one main glacier and three reservoirs (Figure 2). The main glacier, which is the source of the river, presently covers an area of 1.62 km2. Taking this glacier as the center, 109 large and small modern glaciers are developed within an area of 300 km2, comprising a total glacier area of 38.3 km2. Among them, 22 glaciers are distributed within the study area, covering a total area of 9.7 km2. The annual runoff of glacier meltwater into the trunk stream is 17.69 million m3. Therefore, glaciers are not only an important water source in the study area but they also represent a "solid reservoir" regulating runoff. Reservoir 1 (Figure 2) is located in the upper reaches of the study area. It is a water conservancy project that has flood control and irrigation as its primary and secondary purposes. Its operation is mainly divided into reservoir closure, dam flood control, and sluice water storage. The dam crest elevation is 2189.2 m, and the total storage capacity is 69.9 million m3. The highest dam height is 98 m. Reservoir 1 can exploit its flood-retention capability to cut the peak flood and regulate and store the flood during the flood season. Reservoir 2 (Figure 2) is located in a mountain depression of the middle reaches of the study area. It is a series regulating reservoir upstream

of the urban area. It undertakes the three tasks of flood control, drought relief, and urban water supply for the city. The catchment area above the section of Reservoir 1 is 2596 km2, including 1070 km2 in the mountainous area of the main stream, 950 km2 in the piedmont plain area, and 576 km2 in the mountainous area of tributary gullies. Reservoir 3 (Figure 2) is located 6 km downstream of Reservoir 2. It is a through-injection reservoir built using natural depressions. The current water surface area is 4.5 km<sup>2</sup> and the maximum storage capacity is 53 million m3. It is a large reservoir for comprehensive flood discharge irrigation, power generation, fish aquaculture, and urban water supply.

To obtain accurate data, a number of investigations were conducted in the study watershed during 2016–2017. These investigations included discussion with local research departments, data acquisition through visits to local government departments, collection of water and sediment samples, and investigation of pollution sources. Thus, data regarding the 34 evaluation indicators were obtained. Details were given in Table A2, which described the source of original data. Tables A3 and A4, respectively, introduced the weights and reference standards of each indicator.

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

*3.1. Results of the Ecological Security Assessment*

Through evaluation of the SEEC model, it was determined that the ESI value of the study area is 80.9–94.2. In accordance with the geographical characteristics, the watershed was divided into three sections:


The assessment results indicate that the status of the entire watershed was in the "safe" state, as defined in Table 1, indicating that the ecological security level was high. The ESI values of the entire study area are shown in Figure 4. It can be seen that the highest ESI values in the entire watershed are in the upstream area, and that the lowest ESI values are near Reservoir 2 in the downstream area. The ESI values of the midstream area are at an intermediate level. To identify the causes of the low ESI scores, we examined the scores for the four evaluation factors in the SEEC model using radar charts, as shown in Figure 5. It can be seen that the scores for the four evaluation factors are uneven. The lowest score is for C, indicating deficiencies in watershed regulation and management, and the second lowest value is for E (ecological health), indicating the need for attention to improve the state of ecological health. The scores for S and E (ecological services function) are high, indicating that the current socioeconomic factors have not impacted negatively on the ecological security of the watershed. Although the score for ecological health is low, it might not have affected the ecological services function of the watershed. Nevertheless, to improve the ESI score of the watershed, factors C and E (ecological health) need to be the focus of attention.

**Figure 4.** Ecological security index (ESI) values throughout the entire study area.

**Figure 5.** Radar diagram for the four factors in the SEEC model.

#### *3.2. Score for Each Evaluation Factor*

(a) Assessment Results of Socioeconomic Impact (S)

The score for S is approximately 96, indicating the positive role it plays in ensuring ecological security. The score distribution of the 10 indicators is shown in Figure 6a. Among the 10 indicators, those with relatively low scores are population growth rate and per capita GDP. The lowest values of both are in the upstream and midstream areas, which are regions with poor natural conditions, sparse population, and low per capita GDP. In the past, it was often believed that if the population growth rate and per capita GDP were low, the pressure on ecology would be small and ecological security would not be threatened. The assessment reveals that levels of population growth rate and per capita GDP that are too low are not conducive to ecological security and thus they should be maintained at reasonable levels.

**Figure 6.** Radar diagrams for the four factors in the SEEC model: (**a**) socioeconomic impact assessment results, (**b**) ecological health assessment results, (**c**) ecological services function assessment results, and (**d**) control management assessment results.

(b) Assessment Results of Ecological Health (E)

The score for E (ecological health) is about 78, indicating the negative role it has in ecological security. The score distribution of the 13 indicators is shown in Figure 6b. Among the 13 indicators, the value for the total nitrogen in the water is too high, followed by the comprehensive nutritional index, organic matter, and heavy metal risk index. In terms of the spatial distribution, the comprehensive nutritional index in the upstream area is relatively low, but it increases gradually from the midstream area to the downstream area, reaching its highest value near Reservoir 3, which mainly reflects the intensity of human activity and tourism development [57]. The heavy metal risk index in the upstream and midstream areas is low, while the highest value is near Reservoir 2. The evaluation results of the comprehensive nutritional index and heavy metal risk index are shown in Figure 7.

#### (c) Assessment Results of Ecosystem Services Function (E)

The calculation result of E (ecosystem services function) is approximately 86, indicating its negative role in ecological security. The score distribution for the four indicators is shown in Figure 6c. The main reason for the low score for E is the low coverage of forest and grass, especially in the upstream and midstream areas, which mainly reflects the relatively intense water and soil erosion in these areas in recent years. Moreover, human activities such as livestock grazing and free-range poultry breeding have aggravated grassland degradation [58]. Additionally, the natural shoreline rate of the river is also low, which is mainly attributable to the construction of embankments, diversion channels, and other human projects in the downstream area [59].

**Figure 7.** Scoring results for the comprehensive nutritional index (**left**) and potential heavy metals ecological risk index (**right**).

(d) Assessment Results of Control Management (C)

The score for C is approximately 64, indicating its negative role in ecological security. The score distribution for the seven indicators is shown in Figure 6d. The low rate of treatment of rural domestic sewage is the main problem in relation to the low C score. The population of the watershed is sparse and scattered, meaning that the high cost of construction of environmental infrastructure make it difficult to collect and treat domestic sewage. The score for the soil and water erosion control rate is the second biggest problem in relation to the low C score. Water and soil erosion in the upstream and midstream areas is serious, and the grassland is deteriorating steadily. The effect of implemented mitigation measures has been limited because the dry climate and steep terrain in the upstream and midstream areas are not conducive to the control of water and soil erosion.

#### *3.3. Applicability Evaluation of the SEEC Model*

To evaluate the applicability of the SEEC model, factors that could characterize the level of watershed ecological security were selected. If the assessment results of the SEEC model indicated that these factors were relevant, then the evaluation method could be considered to have satisfactory applicability. In assessment of watershed ecology, water quality status is often used as an important comprehensive assessment factor with which to characterize whether the current ecology is threatened, i.e., whether the ecological status is safe. Therefore, the water quality of the study watershed was selected as the assessment factor for evaluation of the applicability of the SEEC method. To keep the assessment independent, we gave priority to the correlation analysis between the selected water quality indicators for evaluation (hereafter, the water evaluation indicators) and the water quality indicators participating in the evaluation of the SEEC model (hereafter, the water model indicators). The water evaluation indicators comprised pH, conductivity, biochemical oxygen demand, petroleum, volatile phenol, mercury, lead, copper, zinc, fluoride, selenium, arsenic, cadmium, hexavalent chromium, cyanide, anionic surfactant, sulfide, fecal Escherichia coli, sulfate, chloride, nitrate, and suspended solids. The water model indicators comprised dissolved oxygen, transparency, ammonia nitrogen, total phosphorus, total nitrogen, the permanganate index, mineralization of water, chlorophyll a, and the comprehensive nutritional index. The water quality data of the evaluation indicators were obtained from the local environmental monitoring department.

First, the monitoring data at the same point and for the same period were selected, and IBM SPSS 20.0 software was used to calculate the Pearson correlation coefficient of the two datasets (i.e., water evaluation indicators and water model indicators). Evaluation indicators that showed obvious correlation with the water model indicators were eliminated. Additionally, because the SEEC model included the heavy metal risk index of sediment, the heavy metal indicators were also removed from the water evaluation indicators to avoid affecting the independence of the index. Through the screening process, the evaluation indicators were determined as follows: pH, petroleum, volatile phenol, fluoride, cyanide, anionic surfactant, sulfide, fecal Escherichia coli, sulfate, chloride, and suspended solids. There was no obvious correlation between the final water evaluation indicators and the water model indicators, indicating that the final water evaluation indicators were independent and could be used to evaluate the applicability of the SEEC model.

Second, using the water evaluation indicators, the Nemero index method [60] was used to evaluate water quality, and the evaluation results were expressed in terms of China's surface water environmental quality standard. The evaluation results were characterized as the higher the score, the worse the water quality. Obviously, when the water quality of the watershed is poor, its ecological security level is low. In other words, under natural conditions, the water quality score results should be negatively correlated with the watershed ecological security index. If the evaluation result of this paper also conforms to the above rules, it shows that the evaluation result of this paper is relatively accurate, and the evaluation method has good applicability.

Then, ArcGIS 10.3 was used for spatial interpolation to obtain the spatial distribution characteristics of water quality, as shown in Figure 8. Finally, the Pearson correlation coefficient between the data illustrated in Figure 8 and the SEEC evaluation results (Figure 4) was calculated using the ArcGIS 10.3 spatial analysis module. Through calculation, the Pearson correlation coefficient was determined to be approximately −0.4, which was consistent with the rules described above, indicating that the comprehensive evaluation results of water quality were better in areas with high ecological security scores. The evaluation results showed that the SEEC model has satisfactory performance when applied to evaluation of watershed ecological security.

**Figure 8.** Comprehensive evaluation results of water quality in the study area in 2014 (the lower the score, the better the water quality).

#### **4. Conclusions**

Using the framework of the DPSIR model, and considering the connotation of watershed ecological security assessment, this article established the SEEC model that incorporates four factors: socioeconomic impact (S), ecological health (E), ecosystem services function (E), and control management (C). The SEEC model contains 34 evaluation indicators. We selected a watershed in the arid region of Northwest China for application of the research method. Through evaluation, the ESI value of the study area was approximately 81–94. Through analysis of the evaluation results, it was elucidated that the factors leading to the low score were mainly the low rate of treatment of rural domestic sewage and the low rate of mitigation of soil and water erosion. The results of the evaluation of the applicability of the SEEC model showed that the SEEC model has satisfactory performance in evaluation of watershed ecological security.

This article is successful from these aspects: it puts forward the method approach for watershed ecological security assessment and gives the specific evaluation index, index weight, and other key information; through application, it identifies the key factors affecting the ecological security of the study area, which has good guiding significance for local environmental management and can also provide reference for similar studies. However, owing to limited research funds and other reasons, this article fails to verify the SEEC model, such as through comparison of analysis of the ecological security status in different years. Therefore, more in-depth research is needed regarding the applicability and verification of the SEEC model.

**Author Contributions:** Conceptualization, F.Y., B.W. and G.C.; methodology, B.W. and F.Y.; software, B.W.; validation, B.W.; formal analysis, B.W.; investigation, B.W.; data curation, Y.T.; writing—original draft preparation, B.W.; writing—review and editing, Y.T. and M.Z.; visualization, B.W. and D.Z.; supervision, Y.T.; project administration, G.C. and D.Z.; funding acquisition, M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China (2019YFC1804404), Beijing Advanced Innovation Program for Land Surface Science of China and Hebei Province Water Conservancy Science and Technology Project (2018-09).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** During the investigation, the local Environmental Protection Bureau, local environmental monitoring station, local Agriculture and Animal Husbandry Bureau, local Water Affairs Bureau, local Health Bureau, local water industry group, local Institute of Ecology and Geography, and the Chinese Academy of Sciences provided help regarding the data used in the study. In the process of establishing research methods and analyzing data, experts in environmental science, hydrology, water ecology, biology, and statistics from the local Institute of Ecology and Geography, Chinese Academy of Sciences, Ecological Center of the Chinese Academy of Sciences, Chinese Academy of Environmental Sciences, and other institutions provided assistance. We thank James Buxton, from Liwen Bianji (Edanz) (accessed on 29 November 2021), for editing the English text of a draft of this manuscript.

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

#### **Appendix A**

**Table A1.** Definition, reasons for selection, and calculation or acquisition method of each indicator.



**Table A1.** *Cont.*


**Table A1.** *Cont.*


**Table A1.** *Cont.*


**Table A1.** *Cont.*

**Table A2.** Data acquisition process.



#### **Table A2.** *Cont.*


**Table A3.** Weight of each factor and indicator.

**Table A4.** Reference standard and basis for determination of each indicator.



#### **Table A4.** *Cont.*

#### **References**


### *Article* **Research on Water Resources Allocation System Based on Rational Utilization of Brackish Water**

**Dasheng Zhang 1,2,3, Xinmin Xie 2, Ting Wang 1,2, Boxin Wang <sup>3</sup> and Shasha Pei 3,\***


**Abstract:** The rational utilization of unconventional water sources is of great significance to areas where conventional water resources are scarce, and water resource allocation is an important way to realize the rational distribution of multiple water sources. This paper constructs a water resources allocation system integrating model data parameter database, water resources supply and demand prediction module, groundwater numerical simulation module and water resources allocation module. Taking brackish water as the main research object and final goal of achieving the best comprehensive optimization of social, economic and ecological benefits. The brackish water is incorporated as an independent water source into the water resource allocation model, and the stratum structure model and groundwater numerical model are constructed to simulate the brackish water level in the planning target year. The water resources allocation system is applied to Guantao County, China. The results show that increasing the development and utilization of brackish water under the recommended scheme can significantly reduce the water supply pressure of local fresh water resources in agriculture and industry. Compared with the current year, the overall water shortage in the region will be reduced by 4.493 <sup>×</sup> 106 m3 in 2030, and meanwhile, the brackish water level will be decreased by 12.69 m in 2035, which plays a positive role in improving soil salinization.

**Keywords:** water allocation; brackish water; groundwater numerical simulation; general algebraic modeling system

#### **1. Introduction**

Water resources are one of the important natural resources that protect people's livelihood and promote social and economic development, especially for developing countries that have a large population and need to improve their water-saving level and practices. Since the 21st century, rapid population growth and rapid development of the economy has increased the urgent demand for the quantity and quality of water resources [1]. The increasing demand for water has gradually exceeded the supply capacity of local water resources, resulting in a large-scale water resources shortage [2,3]. Especially for agriculture, under the current situation of insufficient local surface water, it is often necessary to use large amounts of groundwater to maintain the normal growth of crops and meet the needs of agricultural production, which also causes problems such as groundwater level decline and land subsidence [4], threatening the safety of human water use and the benign operation of the natural water cycle. There is an urgent need to balance water supply and use, in order to eliminate the increasingly fierce competition of water use and deterioration of the water ecological environment [5,6].

Water resources allocation refers to the planning for macro-control of water resources according to historical water inflow data and user water demand. It is considered to be one

**Citation:** Zhang, D.; Xie, X.; Wang, T.; Wang, B.; Pei, S. Research on Water Resources Allocation System Based on Rational Utilization of Brackish Water. *Water* **2022**, *14*, 948. https:// doi.org/10.3390/w14060948

Academic Editors: Yuanzheng Zhai, Jin Wu, Huaqing Wang and Pankaj Kumar

Received: 30 December 2021 Accepted: 10 March 2022 Published: 17 March 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/).

of the most effective methods to coordinate the balance of social, economic, and ecological water use [7]. With researchers globally conducting in depth research, the water resources allocation model has been continuously developed [8,9]. Minsker, et al. [10] constructed a multi-objective allocation model of water resources based on genetic algorithm through the uncertainty of hydrological factors, and characterized the uncertain factors in the process of water resources allocation. Rosegrant, et al. [11] coupled the hydrological model with the economic development model to evaluate social and economic benefits brought by the improvement of water resources allocation and utilization efficiency. Chen, et al. [12] formulated the interval multi-level classified water resources allocation model to optimize water resources allocation to the municipal scale, and the allocation objects are industry, hydropower, and agriculture. In addition, with the development of computer technology, optimization algorithms, artificial intelligence, and space technology are increasingly used in the field of water resources allocation, such as artificial fish school algorithm [13,14], particle swarm optimization algorithm [15,16], and neuro-fuzzy reinforcement learning method [17], which greatly optimizes the scientific nature and accuracy of the water resources allocation model. In addition, it also drives the water resources allocation model to gradually transit to multi-disciplinary integration, so as to serve the solution of multiobjective problems. In the study of optimal allocation of unconventional water resources, Li, et al. [18] developed an interval number hierarchical planning model to integrate water into the water allocation system, taking into account the interests of decision-makers at different levels and the uncertainty of the water allocation system, and to guide the conventional and unconventional water supply sectors to supply water rationally and efficiently. Mooselu, et al. [19] proposed a method for optimal allocation of wastewater reuse water with the objective of minimizing the cost of water supply and regional water shortage. Gao, et al. [20] proposed an optimal allocation method for wastewater reuse water with the objective of minimizing the consumption of freshwater resources and the energy consumption of the water supply process by incorporating medium and desalinated water into the water allocation study, and used NSGA-II. The optimal water and energy-saving scheme is given on the Pareto boundary by solving this multi-objective problem using NSGA-II.

Brackish water refers to the water resources buried underground with a salinity of 2–5 g/L [21]. With the increasingly severe shortage of conventional water resources, it is urgent to bring brackish water into the optimal water resources allocation. China is rich in brackish water resources. According to statistics, the amount of brackish water resources in China is 27.7 billion m3, of which 13 billion m3 can be utilized. Underground shallow brackish water is mainly distributed in the Northern China Plain, northwest arid area and eastern coastal area, with a buried depth of about 10–100 m, which is easy to be exploited and utilized [22]. As one of the main grain production areas in China, the use of brackish water for agricultural irrigation in the Northern China Plain has very important practical significance [23]. The use of brackish water for farmland irrigation has a long history worldwide. In the southwest of the United States, the use of brackish water for irrigation of cotton, sugar beet and other crops basically does not affect the crop yields [24]. Israel, Australia, Spain, Tunisia, and other countries all have carried out practices related to brackish water irrigation and achieved good results [25]. In addition to the relevant research on the impact of brackish water irrigation on crop yield, other scholars have also carried out relevant research on the impact of brackish water salinity on soil infiltration characteristics [26], soil water and salt transport model [27], soil water and salt regulation measures [28–30], further expanding the depth and breadth of the research on the internal mechanism of the impact of brackish water irrigation on soil and crop growth.

In summary, this paper provides a comprehensive review of the research results on water resources allocation, non-conventional water resources allocation and brackish water utilization, etc. The limitations of the existing water resources allocation research are mainly reflected in the following two aspects: Firstly, the focus of the current research on optimal water resources allocation is to reasonably allocate the water supply from conventional water sources, medium water and desalinated water through various engineering and technical measures, so as to meet the demand of different users. Secondly, in terms of brackish water utilization, the current research mainly focuses on the impact of brackish water irrigation on crop growth, yield, and land quality by conducting irrigation experiments, while there are relatively few studies on the prediction and simulation of water level changes after the large-scale application of brackish water.

Based on the above discussion, the main objective of this study is to construct a water resources optimization allocation system, the core of which is to couple numerical groundwater simulation with water resources optimization allocation, to reasonably allocate brackish water as an independent water source, to realize the efficient utilization of brackish water resources in agricultural irrigation and industrial water use, and to input the allocation results into the numerical groundwater model to predict the changing trend of brackish water level. This study has certain reference value for reducing the water supply pressure of conventional water sources and inter-basin water transfer in agricultural irrigation and industrial water use, improving the development and utilization rate of brackish water resources, and realizing the synergistic allocation of conventional water resources and brackish water resources in water-scarce areas.

#### **2. Material and Method**

#### *2.1. Study Area*

Guantao County is located in the south of the Northern China Plain (longitude 115◦06 –115◦40 , latitude 36◦27 –36◦47 ), with a total area of 456.3 km2. It borders Shandong Province with Wei Canal in the east, Guangping, Quzhou and Qiu Counties in the west, Daming County in the south and Linxi County in the north. The main river in the territory is Wei Canal. Its annual average surface water resource is 4.645 million m3, its groundwater resource is 57.654 million m3, and the total amount of water resources is 62.299 million m3. By the end of 2019, the total population was 350,900, including 74,100 urban residents and 276,800 rural residents, with an urbanization rate of 21.1%. The gross national product is CNY 7.095 billion, and the proportion of the three industries is 6.2:6.3:11.1. The geographical location of Guantao County is shown in Figure 1.

**Figure 1.** Geographical location of Guantao County.

Guantao County is rich in shallow brackish water reserves. According to the relevant results of Evaluation of Water Resources of Guantao County, Hebei Province, the salinity of Guantao County is 2 g/L < M ≤ 3 g/L, and the multi-year average groundwater resources of brackish water is 15.15 million m3. The salinity is 3 g/L < M ≤ 5 g/L, and the multi-year average groundwater resources of brackish water is 4.05 million m3. The distribution of brackish water in Guantao County is shown in Figure 2. It can be seen from the figure that groundwater with salinity M ≤ 2 g/L is mainly distributed in Shoushansi Township and Chaibao Town in the south-central part of Guantao County. The groundwater with salinity 3 g/L < M ≤ 5 g/L is mainly distributed in Luqiao Township and Weisengzhai Town in the north of Guantao County, Chaibao Town in the middle and Fangzhai Town and Wangqiao Township in the south. The salinity of other areas is 2 g/L < M ≤ 3 g/L.

**Figure 2.** Distribution of brackish water in Guantao County.

#### *2.2. Optimal Water Resources Allocation System*

The optimal water resources allocation refers to the spatial and temporal allocation of limited water resources through various engineering or non-engineering measures in the basin or specific area, in accordance with the principle of natural sustainable development, so as to meet the water demand of each area to the maximum extent and coordinate the contradictions of water users without affecting the ecological environment. This aims to maximize the social and economic benefits of water resources, and to promote the sustainable and stable development of the river basin, regional economy, and the health and stability of ecological environment. In order to achieve the above objectives, this study constructs the optimal water resources allocation system on the premise of fully considering the water resources endowment conditions in the study area. The main work modules of the system include model data parameter database, water resources supply and demand prediction module, groundwater numerical simulation module, and water resources allocation module. The working relationship between modules is shown in Figure 3.

**Figure 3.** Working relationship of modules of optimal water resources allocation system.

#### 2.2.1. Model Data Parameter Database

The model data parameter database is the data source of the optimal water resources allocation system as well as the basis for the normal operation of other modules. It mainly includes the following three contents:


$$Q\_u - Q\_s - Q\_l - Q\_d = \Delta Q \tag{1}$$

where *Qu* is the upstream water inflow of the node, 10,000 m3; *Qs* is the water supply quantity of the node, 10-thousand m3; *Ql* is the water loss of the node, 10-thousand m3; *Qd* is the downstream water transmission of the node, 10,000 m3; and Δ*Q* is the storage variable of the node, 10,000 m3.

The water resources allocation system network diagram is the generalization and abstraction of water resources allocation. The river system is drawn by approximating the actual river distribution in the study area, river basin and administrative area sections, reservoir nodes, water diversion and lifting nodes, calculation units, lakes, etc., are marked in the river system. Reservoir nodes, water diversion and lifting nodes, calculation units, lakes and other nodes are used to connect local surface water sources, groundwater sources, unconventional water sources and external diverted water sources in the study area. The drainage line of the calculation unit is connected with the river to form a drainage system. The study area is partitioned to ensure that there are water transmission lines in each partition. On the basis of fully collecting the data of the study area and comprehensively considering the brackish water distribution and water use of various departments in Guantao County, Guantao County is generalized into 17 calculation units, five water lifting nodes, 12 control section nodes and one water plant node. The water resources allocation system network diagram in Guantao County is shown in Figure 4.

**Figure 4.** Water resources allocation system network diagram in Guantao County.

2.2.2. Water Resources Supply and Demand Prediction Module

The water resources supply and demand prediction module mainly includes a water demand prediction part and a water supply prediction part, wherein, the water demand prediction part adopts the combination of time series method and index quota method, and comprehensively considers the requirements of social and economic development and water conservation, so as to formulate the water demand schemes for different industries under different social development situations.

Urban and rural residents' domestic water demand forecast is calculated by the formula:

$$\mathcal{W}\_t = p\_t q\_t \tag{2}$$

where *Wt* is the projected water demand of urban and rural residents in level *t*, *pt* is the projected population in level *t*, and *qt* is the projected water quota in level *t*.

Industrial water demand forecast is calculated by the formula:

$$I\_t = \varepsilon\_t r\_t \tag{3}$$

where *It* is the predicted industrial water demand at level *t*, *et* is the predicted industrial value-added water consumption per million yuan at level *t*, and *rt* is the predicted industrial value-added at level *t*.

Agricultural irrigation water demand forecast is calculated by the formula:

$$\mathcal{W}\_{\rm mt} = \frac{\sum\_{j=1}^{m} A\_{tj} \omega\_j}{\eta} \tag{4}$$

where *Wmt* is the gross irrigation water demand in level *t*, *Atj* is the planted area of the *j*th crop in level *t*, *ω<sup>j</sup>* is the irrigation quota of the *j*th crop, *m* is the crop type, and *η* is the integrated irrigation water effective utilization coefficient.

The water supply prediction part is based on the comprehensive investigation and analysis of the engineering layout, water supply capacity, operation status, development and utilization of water resources and existing problems of existing water supply facilities, and analyzes the prospect and potential of water resources development and utilization, so as to formulate water supply planning schemes for different target years, and predict the available water supply quantity of each planning scheme. The accuracy of water demand prediction and water supply prediction has a great impact on the results of water resources allocation. Therefore, it is necessary to closely combine the relevant requirements of social and economic development and water conservancy project planning in the study area so as to avoid the excessive deviation between the prediction results and development objectives.

#### 2.2.3. Water Resources Allocation Module

The water resources allocation module is based on the binary water cycle theory, which takes the water resources system, the social economic system, and the ecological environment system as an organic whole, and efficiently promotes social and economic construction under the premise of ensuring the sustainable and sound development of the ecological environment. To promote the conservation and protection of water resources with the advancement of economy and technology, and the high-quality water resources saved can be distributed to various water users, realizing a virtuous circle of the organic whole. The water resources allocation module adopts the system analysis method based on linear programming, and establishes the balance equation and constraint equation of each control node, reservoir and calculation unit in the water resources allocation system network diagram according to the principle of water balance. It aims at the comprehensive optimization of social benefit goal, economic benefit goal and ecological environment benefit goal to construct the objective function of water resources allocation, and programs by using the general algebraic modeling system (GAMS). The optimal configuration scheme is solved by repeated iteration. The explanations of relevant variables are shown in Table 1, as the constructed objective Equations (2)–(6).


**Table 1.** Variable definition.

#### 1. Social benefit goal

The social benefit goal is mainly reflected by the safety level of water supply. Each water allocation department shall carry out multi-source joint water supply to meet the requirements of each water user for water quantity and quality. Therefore, in the process of allocation, the water supply weight of each water source shall be reasonably set to determine the optimal safety level of water supply.

$$\begin{aligned} \text{Max } F\_1 &= \sum\_j \mathbf{a\_{sur}} \times \left( \mathbf{SC}\_j + SI\_j + SE\_j + SA\_j + SR\_j \right) \\ &+ \sum\_j \mathbf{a\_{grd}} \times \left( \mathbf{GC}\_j + GI\_j + GE\_j + GA\_j + GR\_j \right) \\ &+ \sum\_j \mathbf{a\_{dir}} \times \left( DC\_j + DI\_j + DE\_j + DA\_j + DR\_j \right) \\ &+ \sum\_j \mathbf{a\_{luv}} \times \left( BI\_j + BA\_j \right) \\ &+ \sum\_j \mathbf{a\_{rcc}} \times \left( RI\_j + RE\_j \right) \end{aligned} \tag{5}$$

where *αsur* is surface water supply coefficient, *αgrd* is groundwater supply coefficient, *αdiv* is external diverted water supply coefficient, *αbw* is brackish water supply coefficient, and *αrec* is sewage disposal recycled water supply coefficient.

#### 2. Economic benefit goal

The economic benefit goal is mainly reflected by water shortage rate. For industry and agriculture, the lower the water shortage rate is, the higher the economic benefit of water supply will be. The corresponding water shortage weight is allocated to each industry according to the industrial demand and the economic benefits generated by unilateral water, so as to minimize the regional water shortage and realize the coordinated development of the regional economy.

$$\text{Min}\,F\_2 = \sum\_{j} \left( a\_j^c A C\_j + a\_j^l A I\_j + a\_j^E A E\_j + a\_j^A A A\_j + a\_j^R A R\_j \right) \tag{6}$$

where *α<sup>c</sup> <sup>j</sup>* is urban domestic water shortage coefficient; *<sup>α</sup><sup>I</sup> <sup>j</sup>* is industrial water shortage coefficient; *α<sup>E</sup> <sup>j</sup>* is ecological environmental water shortage coefficient; *<sup>α</sup><sup>A</sup> <sup>j</sup>* is agricultural water shortage coefficient; *α<sup>R</sup> <sup>j</sup>* is rural water shortage coefficient.

#### 3. Ecological environment benefit goal

The ecological environment benefit goal is mainly reflected by lake water shortage. At present, Princess Lake in the study area has shrunk to a certain extent. Therefore, it is very necessary to protect the existing lake area and give full play to its due ecological function.

$$\text{Min } F\_3 = \sum\_{\text{wl}=1}^n a\_{\text{wl}} \times WL\_{\text{wl}} \tag{7}$$

where *αwl* is the water use coefficient of lake and wetland.

4. Final goal

The final goal is to achieve the optimal comprehensive benefits of water resources allocation. In the process of water resources allocation, it comprehensively considers the carrying capacity of water resources and water environment, and gives consideration to the water supply and consumption coordination between multiple water sources and multiple users, so as to meet the water requirements for the healthy development of social economy and ecological environment.

$$F = c\_1 F\_1 + c\_2 F\_2 + c\_3 F\_3 \tag{8}$$

where *F* is the final goal of water resources allocation, and *c1*, *c2*, and *c3* are the coefficients of social benefit goal, economic benefit goal, and ecological environment benefit goal, respectively.

5. Constraints

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩

The constraints mainly include node water balance, supply and demand balance of calculation unit, water supply balance of various water sources (surface water, groundwater, external diverted water, brackish water, water for sewage treatment, etc.), river and canal overflow capacity, river network channel storage, return water balance of calculation unit, lake and wetland constraints, etc. Limited by the length of the article, see [32] for the specific constraint equation expression.

#### 2.2.4. Groundwater Numerical Simulation Module

Numerical simulation is one of the most commonly used methods to solve the problem of groundwater flow. It mainly uses the approximation principle to divide the study area, and transforms the nonlinear partial differential equation into a difference equation through discretization, which can solve the complex groundwater evaluation problem on the premise of ensuring a certain accuracy. The differential equation of groundwater flow is as follows:

$$\begin{array}{ll} \frac{\partial}{\partial x} \left[ K(H - Z) \frac{\partial H}{\partial x} \right] + \frac{\partial}{\partial y} \left[ K(H - Z) \frac{\partial H}{\partial y} \right] + \frac{\partial}{\partial z} \left[ K(H - Z) \frac{\partial H}{\partial z} \right] + \varepsilon = \mu \frac{\partial H}{\partial t} \\ \frac{\partial}{\partial x} \left[ K M \frac{\partial H}{\partial x} \right] + \frac{\partial}{\partial y} \left( K M \frac{\partial H}{\partial y} \right) + \frac{\partial}{\partial z} \left( K M \frac{\partial H}{\partial z} \right) + \mathcal{W} + p = S M \frac{\partial H}{\partial t} \\ H(x, y, z)|\_{t = 0} = H\_0(x, y, z) \\ H(x, y, z, t)|\_{\Gamma\_1} = H\_1(x, y, z, t) & x, y, z \in \Gamma\_1 \quad t > 0 \\ K M \frac{\partial H}{\partial x} \Big|\_{\Gamma\_2} = q(x, y, t) & x, y, z \in \Gamma\_2 \quad t > 0 \end{array} \tag{9}$$

where *H* is water level, m; Z is floor elevation of the first phreatic aquifer, m; *K* is hydraulic conductivity of aquifer, m/d; *ε* is rainfall infiltration and agricultural regression intensity, m/d; *μ* is water yield of the first phreatic aquifer; *W* is overflow strength, L/d; *P* is mining intensity per unit volume of aquifer, L/d; *S* is water storage rate of confined aquifer, L/m; *H*<sup>0</sup> is initial head, m; Г<sup>1</sup> is class I head boundary; *H*<sup>1</sup> is class I boundary water level, m; Г<sup>2</sup> is class II flow boundary; and *q* is boundary flow, m2/d.

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

#### *3.1. Establishment of Groundwater Numerical Model*

According to the hydrogeological conditions of Guantao County, combining with the long-term observation data of groundwater and the actual exploitation of groundwater, the phreatic layer in Guantao County is generalized into a groundwater system with unified hydraulic connection by using the 3D groundwater flow numerical simulation software GMS, and the system is generalized into a heterogeneous and unstable groundwater system whose parameters change with space.

#### 3.1.1. Stratigraphic Structure Generalization

The aquifer in the study area is a loose rock pore aquifer. Its lithology is mainly interbedded with fine sand, silty sand and silty clay, and the formation type is alluvial and lacustrine. The stratum thickness changes uniformly, and there is hydraulic connection between aquifers in the horizontal direction.

According to the spatial distribution of stratigraphic lithologic structure within the study area and combining with the horizon of existing monitoring wells in Guantao County, the strata in the study area are generalized into a three-layer structure composed of two aquifers and one aquitard which are interbedded. Referring to the deep mining depth of the groundwater mining layer in the study area, the bottom boundary depth of the model is set as 160 m.

#### 3.1.2. Generalization of Aquifer Boundary Conditions

1. Lateral boundary conditions.

The southern part of the study area is mainly recharged by lateral runoff. In terms of boundary conditions, it is regarded as class II flow inflow boundary. There is almost no water exchange between the north of the study area and the surrounding counties, and the boundary is perpendicular to the water level line, which can be generalized as a water barrier boundary. Since the overall flow direction of shallow groundwater in Guantao County is from northeast to southwest, the west and the east of the study area can be generalized as constant head boundaries.

2. Vertical boundary generalization.

The upper part of the model is the phreatic surface boundary, which is recharged by atmospheric precipitation, river leakage and agricultural irrigation return, and discharged by phreatic surface evaporation and manual mining. Therefore, the upper boundary is generalized as the water exchange boundary. The bottom of the model has relatively impermeable clay and bedrock, so it is generalized as an impermeable boundary.

#### 3.1.3. Determination of Hydrogeological Parameter Partitions and Initial Values

The hydrogeological parameters of Guantao County are determined according to the relevant hydrogeological maps and hydrogeological exploration data of Guantao County. The phreatic aquifer in Guantao County is a loose rock pore aquifer. Its lithology is mainly composed of fine sand, with NE trending silty sand bands scattered. The study area is divided into four areas. The initial values of hydraulic conductivity and water yields are shown in Table 2, and the hydrogeological parameter partitions are shown in Figure 5.


**Table 2.** Hydrogeological parameters.

**Figure 5.** Permeability coefficient partition (**a**); hydraulic conductivity zoning and observation well locations in the study area (**b**).

#### 3.1.4. Construction and Solution of Numerical Model

1. Spatial subdivision.

The study area is divided into rectangular grid elements by GMS software. After spatial subdivision, the study area is divided into rectangular grid elements, of which about 11,000 rectangular grid elements are in each layer. See Figure 6a for model element subdivision results.

2. Initial flow field.

In this study, the initial time of 1 January 2018 was taken. According to the groundwater level of groundwater monitoring wells, Kriging interpolation method is used to draw the initial water level contour of the study area, as shown in Figure 6b. Indicated by Figure 6b, due to the influence of topography, the groundwater level is generally high in the southwest and low in the northeast. Due to the long-term exploitation of groundwater, an obvious depression cone has been formed in Shoushansi area in the west and Fangzhai Town in the south of the study area.

3. Processing of source and sink items.

The main recharge items include: precipitation infiltration recharge, lateral recharge, irrigation recharge, etc. The main discharge items include: artificial mining, evaporation and lateral discharge. Darcy's Law is used to determine the lateral supply and discharge at the boundary. The precipitation is multiplied by the precipitation infiltration recharge coefficient to obtain the precipitation infiltration recharge. The determination method of irrigation infiltration recharge and evaporation discharge is the same as that of precipitation infiltration recharge. The mining yield is input into the model in a well form.

**Figure 6.** Schematic diagram of network segmentation (**a**); the initial flow field in the study area (**b**).

#### 3.1.5. Model Identification and Verification

The identification period of the model is from May 2018 to May 2019. The identification period of the model is an inversion process. According to the formation lithology characteristics of each parameter partition, and combining with the obtained initial hydrogeological parameters, the calculated water level at the monitoring well is consistent with the measured water level through appropriate adjustment of hydrogeological parameters within a reasonable parameter value range, and the water level change trend is generally consistent with the measured water level change trend. The verification period of the model is from May 2019 to December 2019. The model verification period is a forward process, that is, the parameters are kept unchanged, and observe whether the calculated water level of the monitoring well and the trend are consistent. If the two are consistent and the trend is consistent, it indicates that the model can objectively reflect the actual situation of the study area.

After identification and verification, the comparison of historical changes between the calculated water level and the measured water level of some groundwater monitoring wells in the study area is shown in Figure 7. It can be seen from the Figure that the calculated water level of the monitoring well is well fitted with the measured water level, and the constructed groundwater seepage numerical model can objectively reflect the actual situation of the study area. The model simulation error statistics are shown in Table 3. From Table 3, it can be seen that the model MAE and RMSE are both lower than 0.5 m, with high accuracy, and can be used to predict the groundwater level situation in this region.

**Table 3.** Groundwater model simulation error statistics results.


**Figure 7.** Comparison between measured and simulated groundwater level in partial typical monitoring wells.

#### *3.2. Prediction Results of Social and Economic Development and Water Demands*

This study analyzes the water supply source and water supply quantity of Guantao County in the base year, and predicts the water supply capacity of the study area in the planning target year according to the Surface Water Allocation and Utilization Plan of Hebei Province. Comprehensively considering the local social and economic development plan, total water consumption control index and water efficiency control index, this study adopts the scenario analysis method to formulate two schemes of high growth and moderate growth, respectively, according to the two scenarios of enhanced water saving and moderate water saving. A total of four water demand scenarios are combined to analyze and predict the social and economic development and water demands of Guantao County. See Table 4 for the setting of water demand prediction scheme.

**Table 4.** Setting of water demand forecasting scheme.


#### 3.2.1. Analysis and Prediction of Water Supply

Water supply quantity refers to the sum of gross water supply quantity provided by various water source projects to users, including water transmission loss. The total water supply quantity of the existing water supply facilities in Guantao County in 2019 is 84.7 million m3. Among them, the water supply quantities of surface water, groundwater, external diverted water, brackish water and reclaimed water are 38.25 million m3, 32.5 million m3, 4.95 million m3, 7 million m3, and 2 million m3, respectively, accounting for 45.16%, 38.37%, 5.84%, 8.26%, and 2.36% of the total water supply quantity, respectively. The surface water supply source mainly comes from the Wei Canal and the Yellow River Diversion Project in the east of the study area, the groundwater source comes from 9019 motor-pumped wells in the study area, the external diverted water source comes from the South-to-North Water Transfer Middle Route Project, and the reclaimed water source comes from the reclaimed water of Guantao County sewage treatment plant. The water supply capacity of the current water supply project in the study area is 103 million m3. According to the Surface Water Allocation and Utilization Plan of Hebei Province, by the end of 2025, the Weixi New Canal Connection Project and Shennong Canal Project will be added in Guantao County, and the Matou, Xiaocun and Shenjie pump stations will be started. It is estimated that the water supply capacity of 5.6 million m<sup>3</sup> will be increased. Therefore, the estimated water supply capacity of Guantao County in the planning target year is 108.6 million m3.

#### 3.2.2. Prediction of Social and Economic Development

The main contents of prediction of social and economic development include population prediction, agricultural irrigation area prediction and industrial added value prediction. The short-term (2021–2025), medium-term (2026–2030), and long-term (2031–2035) social and economic development prediction results of Guantao County are shown in Figure 8 under the scenarios of high growth and moderate growth. The results show that in terms of population, Guantao Town, Shoushansi and other areas with relatively high urbanization rates have a large population and a relatively large growth rate; in terms of irrigation area, the agricultural irrigation area in Guantao County is widely distributed on both sides of the Weixi Main Canal, of which the length is long in Luqiao Township and Chaibao Town, so the irrigation area is significantly higher; and in terms of industrial added value, the main industrial enterprises in Guantao County are distributed in Guantao Town and Shoushansi, and the industrial enterprises in the other townships are relatively few.

#### 3.2.3. Prediction of Water Demand

Based on the prediction results of social and economic development, this study adopts the quota method to predict the water demands of different industries in Guantao County. The results are shown in Figure 9, which indicates that the agricultural water demand under each scenario is significantly higher than that of other industries, followed by the urban domestic water demand. With the continuous improvement of China's attention to the ecological environment in recent years, the proportion of water demand of ecological environment is also relatively high. With the development of urbanization and industrialization, the water demands of urban life, industry, and ecological environment are gradually increasing. At the same time, the popularization of agricultural water-saving technology makes the water demand in agriculture and rural areas gradually decline. In addition, due to the relatively small scope of the study area, the water demand in the high growth scenario increases by 2–3% compared with that of the moderate growth scenario, and the water demand in the moderate water-saving scenario increases by 2–4% compared with that of the enhanced water-saving scenario.

**Figure 8.** Social and economic development prediction: (**a**) Population prediction; (**b**) Agricultural irrigation area prediction; (**c**) Industrial added value prediction.

**Figure 9.** Prediction results of water demand in different industries under different scenarios: (**a**) High growth and moderate water-saving scenario; (**b**) Moderate growth and moderate watersaving scenario; (**c**) High growth intensifies water-saving scenario; (**d**) Moderate growth intensifies water-saving scenario.

#### *3.3. Analysis of Optimal Water Resource Allocation Results*

According to the proposed three water demand schemes in different planning target years, which are respectively the high scheme (high growth and moderate water saving), medium scheme (high growth and enhanced water saving) and low scheme (moderate growth and moderate water saving), the water resources allocation results of different schemes in the planning target year of the study area are obtained through the calculation with the water resources allocation module, as shown in Table 5. The relationship between water supply and demand between different water sources and different users is shown in Figure 10.

The overall supply–demand structure of the low scheme is basically reasonable through the reasonable allocation of water sources. The depth of water shortage is the least of the three schemes. It belongs to the way of "determining production by water and taking connotative development". In the planning year, it inhibits the agricultural development of the county so as to control water shortage, and it is a better scheme from the perspective of water resource supply–demand balance. However, the scheme limits the social and economic development of Guantao County to a certain extent and hinders the development strategy of "striving to build Guantao County into a sub-central city with great radiation in the border area of Hebei, Shandong and Henan and the most influential rural tourism destination in the central south plain of Hebei", thus, it is not the preferred allocation scheme.

**Figure 10.** The relationship between water supply and demand between different water sources and different users in 2030 under the second scenario.

The high scheme adopts the mode of "determining water consumption by production and taking epitaxial development". Under the development goal of fully supporting social and economic development and building a sub-central city and ecotourism destination, it controls agricultural water consumption and increases industrial water consumption. Therefore, it is also the scheme with the highest water shortage. Although the quantity of water diverted from Wei Canal and the Yellow River and the water network coverage can be increased by building the water system connection project in the planning target year, considering the difficulty and complexity of the external water transfer project, there are still a lot of feasibility studies and preliminary demonstration work to be invested. Therefore, it is difficult to realize the large surface water supply quantity and scale by 2035. Thus, this scheme is not recommended as a priority scheme but a reservation scheme.

The medium scheme adopts the mode of "supply and demand coordination and steady and healthy development". In the planning target year, it improves the water use efficiency of various industries and comprehensively realizes the balance of regional groundwater exploitation and recharge by controlling groundwater overexploitation. In this way, under the water supply conditions of the established South-to-North water transfer project, it can not only meet the new social and economic, industrial and agricultural water demands, but also save costs and make the rational use of various resources. Therefore, it can be used as a recommended scheme. The results of water resources allocation of each calculation unit in 2030 under the medium scheme are shown in Table 6. It can be seen from Table 6 that by increasing the development and utilization of brackish water, the brackish water agricultural area of Nanxu Village, the brackish water agricultural Chaibao area 2 and the industrial area of Guantao Town will achieve the balance between supply and demand of water resources. For Shoushansi water-saving agricultural area which does not have brackish water resources, the local conventional water resources cannot support the regional agricultural economic development, so it is also the area with the highest water shortage.


**Table 5.**



**Table 6.** The water resources allocation results of each computing unit in the plan in 2030 (106 m3).

Total

 106.454

 33.467

 27.911

 10.000

 0.134

 4.770

 76.282

 4.578

 1.622

 2.413

 64.105

 3.564

 30.172

#### **4. Discussion**

#### *4.1. Results of Brackish Water Allocated According to Industry under the Recommended Scheme*

The water resources allocation model constructed in this paper can calculate the quantity of water supplied by different water sources to water users in various industries in different calculation units. Taking the brackish water allocation according to industry of the recommended scheme as an example, the brackish water allocation results in 2030 are shown in Figure 11 which indicate that only the brackish water in Guantao Town is supplied to industry, and the brackish water in other areas is only supplied to agricultural irrigation. The main reason is that the main industrial enterprises in Guantao County are located in Guantao Town and Shoushansi, and there are no brackish water resources in Shoushansi. Therefore, only brackish water in Guantao Town is used for industrial production. At the same time, brackish water consumption is positively correlated with crop yield in saline water areas. In areas with high crop yields (such as Luqiao Township, Weisengzhai Town, Chaibao Town, Wangqiao Township, etc.), the brackish water consumption is also high. The exploitation of brackish water not only relieves the water supply pressure of shallow groundwater and Wei Canal to a certain extent and alleviates the contradiction between supply and demand of water resources, but also reduces the water level of brackish water

**Figure 11.** Brackish water configuration in 2030 under the second scenario, reduces soil salinization and improves the soil ecological environment.

#### *4.2. Analysis on Variation Characteristics of Brackish Water Level under Recommended Scheme*

The brackish water allocation results of the medium scheme are substituted into the groundwater model for water level simulation, and the variation results of brackish water level in each planning target year are shown in Figure 12. Figure 12 shows that the brackish water level in the study area decreases from −19.57 m in 2019 to −32.26 m in 2035. In 2035, the highest brackish water level in the study area is located in Weisengzhai Town, with a water level of −24.21 m, and the lowest water level is located in Nanxu Village, with a water level of −53.36 m, forming three brackish water cone areas: Chaibao Town, Guantao Town, and Wangqiao Township. The decline of brackish water level can enhance the regulation and storage capacity of soil aeration zone, promote the vertical alternating

movement with shallow groundwater, and realize the desalination of brackish water, so as to increase the quantity of fresh water resources [33]. Therefore, the development and utilization of brackish water under the recommended scheme appropriately reduces the brackish water level, which is conducive to alleviating the demand for fresh water resources for agricultural irrigation and industrial production.

**Figure 12.** Groundwater flow field under high growth and enhanced water saving scheme: (**a**) at the end of 2019; (**b**) at the end of 2025; (**c**) at the end of 2030; (**d**) at the end of 2035.

#### **5. Conclusions**


to agricultural irrigation and industrial production of each calculation unit in the planning target year, so as to realize accurate water resources allocation.

(3) This paper applies the water resources allocation system to Guantao County, Handan City. Through scheme comparison, the medium scheme is finally selected as the optimal allocation scheme. Under this scheme, in 2030, the brackish water agricultural area of Nanxu Village, the brackish water agricultural Chaibao area 2 and the industrial area of Guantao Town will achieve the balance between supply and demand of water resources. Compared with the current year, the overall water shortage in the study area is decreased by 4.493 × 106 m3. Meanwhile, under the recommended scheme, the brackish water level in the study area will drop by 12.69 m in 2035. Therefore, while alleviating the tension of water supply of local conventional water sources, it also reduces soil salinization and realizes the coordinated development of society, economy, and ecology in the study area.

**Author Contributions:** D.Z., methodology, data curation, software; X.X., conceptualization, formal analysis, supervision; T.W., supervision, software; B.W., data curation; S.P., methodology, data curation. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was funded by Key R&D Program of Hebei Province (Grant No: 21374201D), National Key Research and Development Program of China (2021YFC3001005) and Hebei Province Water Conservancy Research Plan (2018-73).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available upon request due to privacy and ethical restrictions.

**Acknowledgments:** The authors would like to express their gratitude to the editors and anonymous experts for their constructive comments.

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

#### **References**


### *Article* **Groundwater Recharge Modeling under Water Diversion Engineering: A Case Study in Beijing**

**Mingyan Zhao, Xiangbo Meng \*, Boxin Wang, Dasheng Zhang, Yafeng Zhao and Ruyi Li**

Hebei Institute of Water Science, Shijiazhuang 050051, China; mingyanzhao99@126.com (M.Z.); wangboxin293@163.com (B.W.); 18501317443@163.com (D.Z.); sjz791979@126.com (Y.Z.); liruyi0117@126.com (R.L.)

**\*** Correspondence: mengxb1994@163.com; Tel.: +86-130-3103-8799

**Abstract:** The influence of surface water resource exploitation and utilization projects on groundwater has been widely studied. Surface water diversion projects lead to a reduction in river discharge, which changes the recharge of groundwater systems. In this study, the numerical simulation method is used to predict the variation in groundwater level under different diversion scale scenarios. The Zhangfang water diversion project in Beijing, China, was chosen for the case study. The downstream plain area of the Zhangfang water diversion project is modeled by MODFLOW to predict the influence of reducing water diversion on the dynamic change in the groundwater level in the downstream plain area. The model results show that the difference in groundwater recharge and discharge on the downstream plain of Zhangfang is 9,991,900 m3/a, which is in a negative water balance state, and the groundwater level continues to decrease. Reducing the amount of water diverted by the Zhangfang water diversion project to replenish groundwater is beneficial to the rise of the groundwater level in the downstream plain area. The results indicate that the groundwater flow model in the downstream plain area of Zhangfang performed well in the influence assessment of surface water resource exploitation and utilization projects on groundwater. This study also provides a good example of how to coordinate the relationship between surface water resources and groundwater resources.

**Keywords:** water diversion project; groundwater; numerical simulation; MODFLOW

#### **1. Introduction**

Groundwater plays an important role in social development and in maintaining the ecological environment [1]. River leakage recharge is an important part of groundwater recharge, so the interaction between the river and aquifer has become a research hotspot [2–4]. The interaction between surface water and groundwater is constantly changing due to the different spatial and temporal distribution of the exchange between surface water and groundwater [5,6]. A water diversion project will change the interaction between the surface water and groundwater, and directly reduce the discharge of the river channel, which is one of the man-made factors causing the change of regional groundwater level and quality [7,8]. Among them, there are many studies on the impact of external water transfer in the water receiving area, including the impact of the South-to-North water transfer project on the groundwater systems of the North China plain [9] and Beijing [10], the impact of watershed ecological water transfer on the hydrological situation in the lower reaches of the Heihe River basin [11], and the impact of the ecological water supplement on the section of Yongding River, Beijing plain [12]. Research shows that water transfer is conducive to the recovery of the groundwater level in the water receiving area. However, there have been few studies on the impact in the water transfer area. The water transfer project directly leads to the reduction of groundwater discharge in the water transfer area [13], resulting in a decline of the downstream groundwater level [14].

The Zhangfang water diversion project is located in the Fangshan district of Beijing, supplying water from the Juma River to Beijing through the Shengtian Canal. In September

**Citation:** Zhao, M.; Meng, X.; Wang, B.; Zhang, D.; Zhao, Y.; Li, R. Groundwater Recharge Modeling under Water Diversion Engineering: A Case Study in Beijing. *Water* **2022**, *14*, 985. https://doi.org/10.3390/ w14060985

Academic Editor: Juan José Durán

Received: 29 January 2022 Accepted: 17 March 2022 Published: 21 March 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/).

1981, Beijing decided to divert water from the Juma River to relieve the pressure of a water shortage. In 2004, the Zhangfang water diversion project was completed, and supplied water to Beijing. Since the construction of the project, the Zhangfang water diversion project has diverted approximately 100 million m<sup>3</sup> of water from the Juma River every year. The large amount of water diversion has reduced the surface runoff of the Nanjuma River and the Beijuma River downstream of the Zhangfang water diversion project, which has caused a serious impact on downstream water resources and the water environment. Therefore, research on the influence of the Zhangfang water diversion project on downstream water resources is conducive to rational water diversion projects and can enrich the research on the influence of water diversion on the water diversion area.

The water diversion project inevitably leads to changes in hydrogeological conditions along the river and changes the supplement and discharge relationship between downstream groundwater and the surface. A variety of methods can be used to quantitatively describe the process: the water balance method [15], Darcy law [16], the groundwater level dynamic method [17], the isotope method [18] and the numerical simulation method [19]. Numerical simulation achieves the purpose of simulating the actual groundwater system by solving the mathematical model. It is one of the most popular methods for solving groundwater-related problems and has been widely used in groundwater resource evaluation [20], the evaluation of the impact of groundwater development [21], the impact of water conservancy projects on groundwater [22], the mutual conversion of groundwater and surface water, etc. Popular groundwater flow numerical simulation codes include the finite element groundwater flow system FEFLOW [23] and modular 3D groundwater flow model MODFLOW based on the finite difference method [24]. MODFLOW-2005 [25] is a recent popular core version of MODFLOW, which uses a block-centered finite difference method to simulate confined and unconfined groundwater flow, as well as external stresses that include regional recharge, pumping wells, evapotranspiration, rivers, etc.

In this study, we use MODFLOW-2005 as a solver to construct a transient threedimensional groundwater flow numerical model in the downstream plain of Zhangfang and study the current groundwater flow of this area. Based on this model, different water diversion amounts from the Zhangfang water diversion project are designed to evaluate the impact of the Zhangfang water diversion project on the groundwater level in the downstream plain area. Additionally, this study reveals the importance of river infiltration for groundwater recharge and proposes suggestions for the rational development and utilization of downstream water resources.

#### **2. Basic information of the Study Area**

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

The study area is located between 39◦35 57" and 39◦4 50" east longitude and 115◦30 40" and 116◦15 26" north latitude. The total area is approximately 2603 km2. The research scope is shown in Figure 1. The climate in this area is a temperate continental monsoon climate. The perennial mean temperature is 13.4 °C, the winter is severely cold from December to January, and the summer is hot from July to August. The mean annual precipitation is 489.9 mm.

#### *2.2. Hydrogeological Conditions in the Study Area*

The study area is located on multistage geological faults in front of the alluvialproluvial fan of the Juma River. The main fault structures are the Huairou-Laishui deep fault and Donglutou-Qiaoliufan fault. West of the Huairou-Laishui deep fault, quaternary deposits directly cover Proterozoic carbonate rock, and the buried depth of the bedrock is generally 10–40 m. East of the fault to the line of the Donglutou-Qiaoliufan fault, underlying the quaternary strata are tertiary cemented conglomerate and glutenite, and the buried depth of the bedrock is generally 20–60 m. To the east of the Donglutou-Qiaoliufan fault line, quaternary deposits suddenly thicken to 170–550 m. The hydrogeological profile in the study area is shown in Figure 2b.

**Figure 1.** (**a,b**) Location and terrain distribution of the study area. (**c**) The location of the river system and hydrological stations in the study area.

The aquifer is mainly recharged by means of lateral runoff in the piedmont zone, and precipitation infiltration and water infiltration of the Nanjuma River and Beijuma River; it is discharged by means of artificial mining and lateral outflow of the aquifer. Affected by topography, groundwater flow in the study area mainly flows from northwest to southeast. In recent years, as a result of the large scale of water through the Zhangfang water diversion project and human exploitation of groundwater, the aquifer has almost drained.

#### *2.3. The Water Diversion Information of the Zhangfang Water Diversion Project*

The Zhangfang water diversion project supplies water to Beijing through the Shengtian Canal, with an annual flow of approximately 100 million m3 from the Juma River. The water diversion amounts can be estimated through the observation data of hydrologic stations. The location of three hydrological stations on the Juma River is shown in Figure 1. Zhangfang hydrological station is located downstream of the Duya hydrological station, with Shengtian Canal between them. Figure 3 shows the monthly average flow data of each hydrological station. According to Figure 3, the discharge data of the Duya Hydrological station is significantly higher than that of the Zhangfang hydrological station, which is caused by the interception of a large amount of water from Juma River by the Shengtian Canal. The Luobaotan hydrological station is located on the Juma River and downstream of the Zhangfang hydrological station. Except for flood season, the flow of Luobaotan hydrological station is almost 0, indicating that the Nanjuma River is in a state of perennial flow failure. The observation data of the hydrological stations show that the Zhangfang diversion project has caused a serious impact on water resources and aquatic ecology in the downstream.

**Figure 3.** Monthly average flow of hydrological station from 2008 to 2016.

#### **3. Methods**

*3.1. Groundwater Flow Numerical Model*

The western boundaries of the study area are the junction of a mountainous area and plain area, which mainly receives lateral recharge in front of the mountain and is generalized as prescribed flow boundaries. The northern boundaries of the study area are perpendicular to the multiyear average groundwater level contour, which can be generalized as no-flow boundaries. The eastern and southern boundaries are mainly the administrative boundaries of cities and counties. Observation wells are distributed near the boundaries, and the variation range of the groundwater level is small; therefore, the boundary can be generalized as prescribed head boundaries. The generalization

of boundary conditions is shown in Figure 2a. Based on the study objectives and the distribution characteristics of the aquifer, we generalized the model as a single-layer model. The upper boundary is the phreatic surface, which accepts external replenishment. Due to the continuous increase in groundwater exploitation, the depth of the groundwater table is greater than 10 m; therefore, the evaporative water loss of the diving surface can be ignored. The bottom boundary of the phreatic aquifer is selected as the model boundary and regarded as a nonflow boundary.

In this study, the flow field on 1 January 2017 was taken as the initial condition of the model. Due to topography, groundwater mainly flows from northwest to southeast. Figure 4a shows the initial flow field in the study area. The study area is located on the alluvial fan of the Juma River, and the groundwater type is quaternary loose rock pore water. According to the survey data, the northwest of the study area is mainly gravel, medium sand, and coarse sand, while the southeast is interbedded with medium sand, fine sand, and silty clay. According to formation lithology and the geological structure, the study area is divided into five regions of hydrogeological parameters. Figure 4b shows the five regions of the hydrogeological parameters. Table 1 shows the value range of hydrogeological parameters in different partitions.

**Figure 4.** (**a**) The initial flow field in the study area. (**b**) The hydrogeological parameter partition of the study area.

**Table 1.** Value range of hydrogeological parameters of the study area.


Groundwater recharge in the study area includes precipitation infiltration, agricultural irrigation regression, lateral inflow, and river leakage. The precipitation infiltration

coefficient in the study area is about 0.11–0.22 and the agricultural irrigation coefficient is about 0.1–0.15. According to the data of the Water Resources Bulletin, the average annual agricultural irrigation water consumption in the study area is about 3.74 million m3/a, and the agricultural irrigation regression water amount is about 56 million m3/a. The lateral inflow is calculated by Darcy s law according to the contour map of the groundwater level. The lateral inflow of piedmont and other inflow boundaries is calculated to be 85 million m3/a. According to the investigation data, the river leakage coefficient in the study area is 0.45. Groundwater discharge in the study area includes exploitation and lateral outflow. According to the Water Resources Bulletin, the average groundwater exploitation in the study area is 461.56 million m3/a. The lateral outflow is calculated by Darcy's law. The distribution of precipitation infiltration, agricultural irrigation and groundwater exploitation is divided according to the administrative boundaries of each county. This is mainly determined by the data of the Water Resources Bulletin that we collected. Figure 1 shows the administrative divisions of Zhuozhou, Yuyuan, Gaobeidian and Dingxing.

According to the variation in the groundwater flow field, the model is extended to a three-dimensional transient groundwater flow model with heterogeneous anisotropy. To simulate the change in groundwater level in the study area, MODFLOW-2005 is used to solve the three-dimensional numerical groundwater flow model in the study area. MODFLOW-2005 uses the finite difference method to calculate the hydraulic conduction between cells. First, the asymmetric conduction matrix is generated, and then the conjugate gradient method is used to solve the problem. MODFLOW-2005 is realized through the GMS platform developed by the Environmental Modeling Research Laboratory of Brigham Young University in the United States. In this study, we set the basic grid size of the model to 500 m×500 m. In the vertical direction, the region is discretized into one layer by using a deformed mesh. A total of 2592 grids exists in the model.

The simulation period is from January 2017 to December 2023, which is divided into 29 stress periods. Each stress period from 2017 to 2018 is one month, with 24 stress periods in total and a time step of 10 days. This time period is used as the identification and verification period of the model. From 2019 to 2023, each stress period is one year, and the time step is 30 days.

#### *3.2. Model Calibration and Validation*

Model calibration is the process of adjusting and improving the model parameter structure and parameter values. The model in this study is corrected by the trial and error method to continuously adjust the hydraulic conductivity and recharge input parameters [26,27] to minimize the error between the simulated value and the observed value [28,29]. The measured groundwater flow field on 26 May 2017 is used as the identification and verification flow field. Figure 5a shows that the simulated flow field and the measured flow field have the same trend and flow patterns. Except on the southwestern side of the study area, which has a poor simulation effect, the other parts all reflect the actual groundwater flow trend.

A total of 27 monitoring wells are used in this calibration process, and their distribution locations are shown in Figure 2a. In this study, Yihezhuang, Xiyi an, Nangaoguanzhuang, Chenzhuang, Liyuzhuang and Chencunying are selected as six typical monitoring holes for display (the rest are not shown). Figure 5b shows the process fitting curve between the calculated water levels and the measured water levels of the selected monitoring wells. The simulated groundwater level process line of the typical monitoring wells is consistent with the measured groundwater level, which accurately reflects the process of groundwater before and after water replenishment.

The parameter partition after identification and correction is shown in Table 2. The hydraulic conductivity changes show a clear trend of gradually decreasing from the upper part of the alluvial fan to the downstream plain region, and the corrected parameter is within its value range. Therefore, the established numerical model can reflect the variation

characteristics of the groundwater in the plain area downstream of Zhangfang. The model can be used to simulate and predict the influence of the Zhangfang water diversion project on the groundwater level in the downstream plain.

**Figure 5.** (**a**) Flow-field fitting diagram on 26 May 2017. (**b**) The fitting curve of typical monitoring wells in the study area.

**Table 2.** Calibrated hydrogeology parameters with the groundwater flow model.


#### **4. Results and Discussion**

*4.1. Groundwater Budget Calculation*

The groundwater budget calculation in the study area from January 2017 to December 2018 is obtained according to the operation results of the groundwater flow numerical model, as shown in Table 3. During the simulation period, the total recharge of groundwater in the study area is 466.8 million m3, the total discharge is 476.79 million m3, and the supplementary discharge difference is −9.99 million m3. The groundwater system is in a negative water balance state. Among the recharge items, the precipitation infiltration is 257.35 million m3, followed by river leakage of 68.77 million m3. They account for 69.86% of the total groundwater recharge in the study area. The discharge is mainly assembled by artificial mining, which is 461.56 million m3, which accounts for 96.81% of the total discharge.


**Table 3.** Mean groundwater budget of the study area during the simulation period (2017–2018).

#### *4.2. Schemes of Reducing Water Diversion in the Zhangfang water diversion project*

4.2.1. Consideration of Different Diversion Schemes of Zhangfang Water Diversion Project

Through model correction and validation, the groundwater flow model can be used to predict the dynamic changes of groundwater level in the downstream plain of Zhangfang in the future. According to the monitoring data of the hydrology station, the Zhangfang water diversion project diverted a large amount of water, which caused a serious impact to the downstream water resources and river ecology. Water diversion projects directly lead to a significant reduction in river discharge, coupled with a large degree of groundwater exploitation; surface water and groundwater interaction gradually developed into a single form of river recharge groundwater. Therefore, the Zhangfang water diversion project should consider the shortage of water resources in the downstream as a whole and make a reasonable allocation of water resources.

In this study, we considered three schemes to predict groundwater level changes in the study area from 2019 to 2023, we set the stress period as one year and the time step as one month in the prediction period.

Scheme 1 maintains the current water diversion volume of the Zhangfang water diversion project at 100 million m3. This scheme can be used to evaluate the change of groundwater resources of the downstream plain of Zhangfang.

Scheme 2 considers reducing the amount of water diverted from the Zhangfang water diversion project by 50%, that is, to discharge 50 million m3/a of water into the South Juma River and North Juma River. The discharge is divided between the South Juma River and North Juma River in a natural water ratio of 7:3. The direct utilization coefficient of surface water is 0.5, and the infiltration coefficient of the river is 0.45. Therefore, the river infiltration replenishment in the downstream plain area increases by 11.25 million m3/a.

Scheme 3 assumes that the water diversion volume of the Zhangfang water diversion project is 0 and analyzes the dynamic influence of groundwater recharge on the groundwater level of the downstream plain after the water volume of the Beijuma River and Nanjuma River increases. In Scheme 3, groundwater infiltration recharge is increased to 22.5 million m3/a.

4.2.2. Prediction of Groundwater Level of the Downstream Plain of Zhangfang in 2023 in Scheme 1

A groundwater flow field in the study area at the end of 2018 is shown in Figure 6a. The prediction of the groundwater level in 2023 is shown in Figure 6b. Figure 6c shows the simulation of groundwater level variation from 2019 to 2023 in Scheme 1. According to Figure 6c, if the amount of diversion of the Zhangfang water diversion project remains unchanged, the groundwater level in the downstream plain of Zhangfang will continue to decline. The water diversion of the Zhangfang water diversion project is about 100 million m3/a, accounting for 1/5 of the total replenishment of the downstream plain area. We think that the amount of water diverted is catastrophic. The region with the largest drop in groundwater level is the northwest of the study area, where the aquifer is

thin and the hydraulic conductivity is large. In the case of unfavorable upstream water inflow conditions, the groundwater in this region is very easy to lose. The maximum drop of water level in this area is more than 10 m.

**Figure 6.** (**a**) contour map of groundwater level at the end of 2018. (**b**) contour map of groundwater level of Scheme 1 at the end of 2023. (**c**) Simulation of groundwater level variation from 2018 to 2023 of Scheme 1.

The current development and utilization pattern of water resources in the downstream plain of Zhangfang is extremely unreasonable. The diversion of the Zhangfang water diversion project reduces the recharge source of groundwater and intensifies the drop of groundwater level. At the same time, the Zhangfang water diversion project leads to the interruption of downstream river flow, which causes serious harm to water ecological security. Therefore, we should consider reducing the amount of water diverted by the Zhangfang water diversion project and rationally allocate the water resources.

4.2.3. Prediction of Groundwater Level of the Downstream Plain of Zhangfang in 2023 in Schemes 2 and 3

Figure 7a,b shows the contour maps of the groundwater levels at the end of 2023 of the Schemes 2 and 3. In order to quantitatively describe the water level changes from 2018 to 2023., Figure 8a,b shows the groundwater level variation from 2018 to 2023 of Schemes 2 and 3.

**Figure 7.** (**a**) contour map of groundwater level of Scheme 2 at the end of 2023. (**b**) contour map of groundwater level of Scheme 3 at the end of 2023.

**Figure 8.** The simulation of groundwater level variation from 2018 to 2023: (**a**) Scheme 2; (**b**) Scheme 3.

According to Figure 8, compared with Scheme 1, reducing the water diversion amount of the Zhangfang water diversion project to recharge groundwater can alleviate the decline rate of the groundwater level in the whole study area. Although the water level in the eastern part of the research area still continues to decline, the maximum decline range is reduced to 2–4 m. The water level in the northwest of the study area increased greatly, and the groundwater level in the water receiving areas of the South Juma River and North Juma River also increased to a certain extent. Figure 8 shows that the unreasonable allocation of water resources in the Zhangfang water diversion project has caused serious impacts on the downstream water resources. The pressures resulting from the water shortage in the downstream can be alleviated by reasonably reducing the water diversion.

According to the operation results of Schemes 2 and 3, the calculation results of the groundwater budget in the study area are obtained. According to Table 4, the river leakage in Scheme 2 increases by 11.25 m3/a, and the supplementary discharge difference is 0.13 million m3/a. The groundwater system is in a positive water balance state. In Scheme 3, the river leakage increases by 22.5 million m3/a, and the supplementary discharge difference is 8.31 million m3/a. The groundwater system is in a positive water balance state. Therefore, in order to ensure the sustainable utilization of groundwater in the study area, we propose to reduce at the diversion water amount of the Zhangfang water diversion project by at least 50%to recharge the groundwater.

**Table 4.** Groundwater budget of the study area of Schemes 2 and 3.


According to Figure 8b, when the water diversion amount of Zhangfang water diversion project is 0, the groundwater level of the South Juma River and North Juma River receiving area rises significantly, but the water level near Baigou River still drops. The reason for this is that we only recharge groundwater though natural rivers. Therefore, when a large number of surface water sources can be used to recharge the groundwater, the replenishment amount of the surface water increased by connecting channels.

#### **5. Conclusions**

By means of groundwater monitoring and numerical simulation, this paper studied the impact of the Zhangfang water diversion project on the downstream groundwater level. Through research, the following conclusions are obtained.

The groundwater flow model in the downstream plain area of Zhangfang is established. With the flow model, it is concluded that the Zhangfang water diversion project reduces the river discharge and the river infiltration recharge in the study area, which places the supplement and discharge difference in a negative water balance state.

If the Zhangfang water diversion project maintains the current water diversion amount of 100 million m3/a, the water level of the downstream plain of Zhangfang will continue to decline. According to the prediction results, the groundwater level of the South Juma River and the North Juma River will rise to a certain extent by reducing the amount of the Zhangfang water diversion project to recharge the groundwater. The water level of the Baigou River will continue to decline, but the decreasing rate will be effectively alleviated. For the sustainable development of the region, it is recommended to reduce the amount of water diverted by the Zhangfang water diversion project by at least 50%.

Through numerical simulation, the dynamic influence of the Zhangfang water diversion project on the change of groundwater levels in the downstream is discussed from the perspective of water resources allocation, which has practical significance. The prediction scheme of the model can provide suggestions for decision- making.

**Author Contributions:** M.Z., conceptualization, formal analysis, investigation, methodology, software, and writing–original draft; X.M., conceptualization, formal analysis, and supervision; B.W., data curation and methodology; D.Z., data curation and software; Y.Z., data curation and Investigation; R.L., software. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study is funded by the Hebei Province Water Conservancy Science and Technology Project (Grant No.: JSKY2021-KY010). Research on key technology of water quality safety in underground hydraulic Mining and restoration in Miyun District, Beijing (2018000020124G037).

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data available upon request due to privacy and ethical restrictions.

**Acknowledgments:** The authors would like to express their gratitude to the editors and anonymous experts for their constructive comments.

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

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-3470-1