Next Article in Journal
Managed Aquifer Recharge as a Low-Regret Measure for Climate Change Adaptation: Insights from Los Arenales, Spain
Previous Article in Journal
The Use of Lime for Drainage of Cohesive Soils Built into Hydraulic Engineering Embankments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Application of Particle Trace Morphology and Sensitivity Analysis in Delineation of Drinking Water Protection Zone in the Luan River, North China

1
Institute of Hydrogeology and Environmental Geology, Chinese Academy of Geological Sciences, Shijiazhuang 050061, China
2
Key Laboratory of Groundwater Sciences and Engineering, Ministry of Natural Resources, Shijiazhuang 050061, China
3
China Water Huai Planning, Design and Researh Co., Ltd., Hefei 230000, China
*
Author to whom correspondence should be addressed.
Water 2022, 14(22), 3702; https://doi.org/10.3390/w14223702
Submission received: 8 October 2022 / Revised: 11 November 2022 / Accepted: 12 November 2022 / Published: 16 November 2022
(This article belongs to the Section Urban Water Management)

Abstract

:
The appropriate division of underground drinking water source protection zones is a low-cost method of preventing water source pollution and ensuring a supply of safe drinking water. Based on FEFLOW, a groundwater flow model of large water sources was established for Luan River, North China. Trace lines of particle reverse migration for 100 and 1000 days were obtained by random walks. According to the trace morphology, the water sources in the riverside water source area were divided into four categories. The first- and second-grade protection areas were delimited by ArcGIS, with areas of 0.375 and 1.20 km2. The local and global sensitivity of the permeability coefficient (K) and effective porosity (ne) effects on the area of groundwater protection zones were calculated. The area of the protection zones was positively correlated with K and negatively correlated with ne. The variation in the protected zone caused by the simultaneous changes in K and ne is the same as that of ne alone, and the global sensitivity is closer to the local sensitivity of ne. This indicates that ne has a greater impact than K on the scope of groundwater protection zones. Moreover, global sensitivity is not simply a superposition of local sensitivity, and the interaction between parameters can reduce the effect of a parameter acting alone on the delineation of protection zones. This also shows that the global sensitivity is closer to the actual situation than the local sensitivity, thus providing a scientific basis for the delimitation and monitoring of water source protection zones.

1. Introduction

As the main water source in arid and semi-arid regions, groundwater resources play an important role in guaranteeing drinking water safety, supporting social and economic development, and maintaining ecological balance [1,2]. Riverside water intake has been used worldwide to assess water safety as it can reference both surface water and groundwater [3,4,5]. The water quality of the water source near the river is easily affected by the quality of the surface water, and the protection of groundwater resources primarily focuses on prevention, with the division of the water source protection area being an important prevention measure. Therefore, the scientific division of groundwater drinking water source protection areas is critical for ensuring the safety of drinking water for residents, preventing groundwater pollution, and protecting the ecological environment of water source areas [6,7].
Developed countries have been working on improving drinking water source protection for more than 100 years. Some representative countries, such as Germany, the UK, and the United States, have established a relatively complete classification system of groundwater source protection areas and legal and regulatory technical systems. Germany was the first country in the world to establish a drinking water source protection zone system. In the 1850s, it promulgated the “Water Law”. At the end of the 1900s, Germany established the first water source protection zone (located in the Koeln area) in the world. At present, the delimitation method involved dividing the water source protection zone into three further zones and specifying the scope and water source utilization of each zone in detail [8]. The UK uses the three-division method to delimit the water source protection area, which is divided into the inner area, outer area, and watershed area. To date, thousands of well source protection areas have been delimited [9]. The United States passed the “Federal Safe Drinking Water Act Amendment” in 1986, requiring the well source protection program (WHPP) [10]. In 1997, the US Environmental Protection Agency (EPA) further put forward very detailed guidelines for the division of water source protection areas. Common methods include arbitrary distance, calculated fixed radius, analytical methods, hydrogeologic mapping, and numerical flow or flow-and-transport models.
Research on the classification of groundwater protection zones in China started later than in developed countries. The “Law of the People’s Republic of China on the Prevention and Control of Water Pollution” promulgated in 1984 stipulated the system of establishing drinking water protection zone for the first time [11]. In 2018, the Ministry of Environmental Protection issued the “Technical Specification for the Division of Drinking Water Source Protection Zones (HJ338-2018)”, which can select different methods for the division of groundwater protection zones according to the hydrogeological characteristics and water source scale of different water sources. Overall, the classification of groundwater protection zones in various countries adopts the three-level classification method, and the corresponding classification method is selected according to the actual situation in China. Common methods include the empirical value method, the formula calculation method, the analytical model method, and the numerical simulation method.
At present, the main development direction is in the application of the numerical simulation method to delimit groundwater protection zones, especially large-scale water source areas [12]. The numerical simulation method can solve many kinds of groundwater flow problems and reflects the migration law of pollutants, allowing for a more accurate delimitation of the protection zone. Zeferino redefined the water source well protection zone located in a disturbed aquifer in the densely populated urban area of Monticho, Portugal through numerical simulation, considering the population growth and the changes in flow, drawdown, and hydraulic gradient after the increase in groundwater exploitation [13]. Amin delineated the well protection zone in Iranshahr using GMS simulations of pollutant sources and migration paths [14]. Staboultzidis simulated the hydrogeological conditions near water source wells in the Keritis River Basin in Greece through the numerical simulation method and completed the division of the protection zones in this area [15]. Gurwin used MODFLOW to simulate the movement of the groundwater flow field in a water source area in southwest Poland and delineated the scope of water source protection zones at all levels based on the distance traveled by protons in the aquifer [16].
Nevertheless, hydrogeological parameters have high spatio-temporal variability, a complex aquifer structure, and uneven spatio-temporal distribution of source and sink items. There are many uncertainties in groundwater numerical models, which often directly affect the quality of the simulation and prediction [17]. Among the many uncertain factors, hydrogeological parameters have a great influence on the model operation results and have thus become the major focus of research into the uncertainties associated with groundwater numerical models. Parameter sensitivity analysis is a common method in uncertainty research, which is important to analyze parameter sensitivity in groundwater numerical models to determine the main controlling factors affecting the groundwater flow field. Sensitivity analysis can evaluate the influence of the parameter uncertainty of groundwater numerical models on simulation results. The analysis methods include local and global sensitivity analyses [18,19]. Local sensitivity analysis is used to test the influence of the change in a single parameter on the results of a groundwater numerical simulation. Local sensitivity analysis is easy to conduct, but the influence of the interaction between parameters on the model output cannot be considered. Moreover, there are certain limitations associated with the calculation and analysis. The global sensitivity analysis considers the influence of different parameter values and the interaction between the model parameters during the calculation process. There are many methods for conducting global analysis, such as the multiple regression method [20], the Morris method [21], the Fourier Amplitude Sensitivity Test [22], and the Sobol’ method based on variance analysis [23]. The Morris method, which is different from other methods with thousands or even tens of thousands of samples, can give the initial value of the parameter, give a particular disturbance (such as 10% or 20%) on the basis of the initial value, and determine the sensitivity of the parameter through the combination of different lines of the matrix.
In light of the above, this study used FEFLOW to establish a groundwater flow model for a large groundwater source near Luan River, and preliminarily delineated the first- and second-grade protection zones of the groundwater source based on particle tracing technology. The Morris method was used to conduct a global sensitivity analysis on the influence of the permeability coefficient (K) and effective porosity (ne) on the area of the protection zone so as to determine the main controlling factors affecting the division of groundwater protection zones in the study region.

2. Materials and Methods

2.1. Geography and Hydrogeology of the Study Area

The study area is located in the western suburb of Chengde, Hebei Province, China, 160 km away from the capital Beijing (40°57′–41°00′ N, 117°40′–117°44′ E). The study area is located in the intermountain valley, and the Luan River passes through it in a nearly northwest–southeast direction. The surrounding mountains are mostly between 400 and 700 m above sea level, and the highest point is 905.4 m. The mountains are mainly composed of Neoarchean Xiaowagou gneiss (Xwgn), Zhongying gneiss (Zgn), and early Proterozoic Shuiquangou porphyritic monzogranite (SPt1). The Luan River Valley extends for about 10 km; the width of the valley is generally 200–300 m, 521 m at the widest point and 172 m at the narrowest point. The elevation of the valley bottom is 368–400 m; the terrain is relatively open and flat, and the cross-section of the valley is generally U-shaped.
The study area has a semi-arid and semi-humid continental monsoon climate. Meteorological details for the study area are shown in Figure 1. Monitoring data from the Luanping county meteorological monitoring station for the past 50 years (1973–2018) show that the annual average precipitation is 545.97 mm. Over the past 10 years (2009–2018), precipitation has been relatively stable with little inter-annual change. Rainfall is distributed unevenly and varies greatly throughout the year, mostly from June to August, which accounts for 60% to 80% of annual rainfall in most years. The evaporation in the study area is large, and the annual average evaporation is 1591.3 mm, which is 2.9 times the annual average rainfall.
The study area is a quaternary single-structure phreatic aquifer, which spreads continuously and stably along the Luan River valley in strips (Figure 2 and Figure 3). The surface of the study area is generally deposited with brownish-yellow clay sand to a depth of 5.5 m. The lithology of the aquifer is mainly sand gravel with a particle size of 2–15 cm and filled with medium-coarse sand. The thickness of Quaternary aquifers is generally 7–9 m and locally 10–14 m. The water-rich area is mainly concentrated on the side beach and floodplain. The water inflow of a single well in this area is 925–4888 m3/d, with an average of 2334 m3/d. The groundwater level in the study area is relatively shallow, with 1.5–3.7 m in the floodplain and the first terrace, and 4.5–7.5 m in the second terrace and the piedmont. The groundwater runoff has clear phreatic characteristics, which is completely controlled by the topography and is consistent with the direction of surface water flow. Good fresh water with TDS less than 500 mg/L is widely distributed in the study area, and the hydrochemistry type is relatively simple, mainly HCO3-Ca.

2.2. Methodology

In this study, the aquifer structure and hydrogeological boundary conditions were determined through boreholes and field investigation. FEFLOW was used to simulate groundwater flow. According to the flow model, using the particle tracking method, ArcGIS was used to delineate the groundwater protection zone and calculate the area. The different impact factors (K, ne) were adjusted to explore their impact on the delineation of the protection zone. Local and global sensitivity analysis methods were used to determine the influence of different impact factors on the classification of the protection zone.

2.2.1. Data Sources

Based on the data from 94 boreholes constructed and 18 boreholes collected, the Kriging interpolation method was used to obtain the contour map of the Quaternary roof and floor elevation in the study area. The groundwater flow field was determined from two water level measurements (December 2019 and June 2020). Parameters such as the aquifer permeability coefficient and rainfall infiltration coefficient were obtained through data collection and field tests such as the pumping test. The series data of rainfall and evaporation from 2010 to 2019 were obtained from the Luanping County meteorological monitoring station, and the dynamic groundwater level data were obtained from the Sandaohe national hydrology station. The groundwater exploitation data were calculated by the quota method by investigating the population, livestock, agricultural planting structure, and area in the study area.

2.2.2. Hydrogeological Conceptual Model

The range of the model is nearly NW–SE along the direction of the Luan River, with a length of about 5.21 km from north to south, a width of 0.25–1.12 km from east to west, and an area of about 4.82 km2. As shown in Figure 4, the lateral recharge of the upstream groundwater is the lateral inflow boundary (boundary “AB”), and the recharge of the gully underflow is also considered as the lateral inflow boundary (boundary “EF”, “GH”, “IJ”, “KL”, “MN”, “OP”). The downstream runoff discharge of the groundwater is the lateral outflow boundary (boundary “CD”). A total of 45 groups of riverbed riser tests were arranged in the field along the Luan River to obtain the vertical permeability coefficient of the riverbed. The results show that the vertical permeability coefficient of the riverbed has a great difference in space: the value is 0.021~8.500 m/d, and the overall variation coefficient is 1.29. The river is considered as the 3rd-kind boundary condition (Cauchy type), and the rest of the boundary is the water-impervious boundary. According to the histogram of 102 boreholes that exposed the Quaternary system, the depth of the Quaternary basement in the study area is 6.0~16.0 m, with the average and median values of 11.0 m and 9.0 m, respectively. The contours of the Quaternary depth is shown in Figure 4. According to the borehole data, Kriging interpolation was used to obtain the Quaternary top and floor elevations, which are the upper and lower boundaries of the vertical direction of the model.

2.2.3. The Governing Groundwater Flow

The recharge, runoff, discharge, and water levels of groundwater in the simulation area change with time, showing the characteristics of unsteady flow. The parameters vary spatially, which reflects the heterogeneity of the groundwater system. The groundwater flow in the simulation area was generalized into heterogeneous and unstable groundwater flow systems. Combined with boundary conditions, the following two-dimensional mathematical model was used to describe the groundwater flow:
μ h t = K ( x , y ) ( h x ) 2 + K ( x , y ) ( h y ) 2 + W
h ( x , y , z , t ) | t = 0 = h 0
K n ( H h ) n = q
where h is the water level elevation of aquifer (m); K(x, y) is the horizontal permeability coefficient (m·d−1); W represents the rainfall, evaporation, or infiltration (m·d−1); μ is the specific yield (dimensionless); h0 is the initial water level distribution of aquifer (m); Kn is the permeability cot of boundary surface in the normal direction (m·d−1); n is the normal direction of the boundary surface; H is the boundary water level elevation (m); and q is the flow rate at the boundary, where inflow is positive, the outflow is negative, and the water impervious boundary is 0 (m·d−1).

2.2.4. Model Input

The model adopts a triangular mesh for discretization. The size of the base grid is 40 m, and it is refined along the river and at the boundary. The size of the refined grid is 20 m. The whole simulation area is dispersed into 7235 plane nodes and 13,526 triangular units. Hydrogeological parameters involved in the model, such as hydraulic conductivity, porosity, rainfall infiltration coefficient, and field infiltration coefficient, were obtained through a borehole pumping test and data collection, and the values varied from 30 to 459 m·d−1, 0.15–0.18, 0.1–0.24, and 0.25, respectively. It was assumed that horizontal isotropy and vertical conductivity were 0.1 times Kx. The division of aquifer permeability coefficient and riverbed vertical permeability coefficient were shown in Figure 5. Precipitation, irrigation return water, and evaporation adopted a unified non-point source sink item. Precipitation and evaporation were based on the measured values of meteorological stations. Irrigation in the study area is mainly concentrated from May to August, and irrigation infiltration was apportioned throughout this period to 10%, 40%, 40%, and 10%, and it was taken as zero in other months. The infiltration from piedmont gullies was distributed according to the monthly precipitation ratio. The groundwater extraction was summarized as the point discharge of 11 populated areas in the study area. The flow field measured in December 2018 was taken as the initial head, and the simulation period was from January to December 2019. The flow field measured in June 2019 and the monitoring data of the long-observation well were used as the basis for model verification.

2.2.5. Model Output

The groundwater flow field measured in June 2019 was used as the basis for model verification. The model error was determined by comparing the simulated water level of the observation well with the measured water level. The residual error of the model fitting was 0.12 and the error was 2.86%. The flow field fitting is shown in Figure 6, and the correlation between the simulated water level of the observation well and the measured water level is shown in Figure 7. The simulation results show that the flow field calculated in the study area conforms to the flow law of groundwater in the area, and there was high consistency between the measured and simulated values. The model can therefore reflect the groundwater flow in the water source area and meet the accuracy requirements for prediction.
The mass balance results provide important information on the quantity and reliability of the groundwater models [24]. The groundwater balance situation in the study area in 2019 (Table 1) was obtained according to the operation results of the groundwater numerical model, and the equilibrium difference in the study area in 2019 was 110.5 × 104 m3, indicating a positive equilibrium. Based on the long-bore dynamic monitoring data, the storage variable of the study area in 2019 was 31.82 × 104 m3. According to the provisions of Water Resources Investigation and Evaluation and Supplementary Rules for Groundwater Resources and Exploitation, the accuracy of equilibrium results was tested based on the relative equilibrium difference, determined as follows:
δ = X/Qtr × 100
where δ is the relative equilibrium difference (%); X is the equilibrium absolute difference, X = QtrQtd ± ΔW; ΔW is the groundwater storage variable (104 m3), which has a negative value for positive equilibrium and positive value for negative equilibrium; Qtr is the total groundwater recharge (104 m3); and Qtd is the total groundwater discharge (104 m3). According to the regulation, ∣δ∣ ≤ 10% means that the equilibrium calculation accuracy meets the requirements.
According to the calculation results of region equalization and storage variables, the relative equalization difference in the study area is 7.22%, which meets ∣δ∣ ≤ 10%, indicating that the calculation results of model equalization meet the accuracy requirements.

2.2.6. Demarcation of Groundwater Protection Zones

The classification of protected areas based on numerical methods is mainly based on the concept of a “hydraulic intercept zone”. Without considering the vertical velocity component, the hydraulic intercept zone refers to the plane area delineated by the horizontal trace of particle migration from the pollution point to the water source well within a given time [25]. Based on the results of stress periods and water exchange between grids obtained by running FEFLOW, the trace lines of water source wells are reversely tracked by random walks. If the trace length is less than two times the well spacing, each well is categorized into the first level protection zone. If the trace length is more than two times the well spacing, the wells are regarded as well groups to delimit the protection area [26]. By delineating the track range of each water source well (group) after 100 and 1000 d (protection zone classification standard), the scope of primary and secondary protection zones was delineated for the numerical simulation method.
The mass balance equation implemented for the tracing was as follows:
x ( n v x ) + y ( n v y ) = W
where n is the effective porosity of rock; vx and vy are the components of the average linear groundwater flow vector on the coordinate axis (m·d−1); and W represents the inflow and outflow amount of source and sink item per unit volume of the aquifer (1·d−1).

2.2.7. Sensitivity Analysis

In this study, local sensitivity and global sensitivity analyses were carried out on the groundwater source protection zone determined by the particle tracking method. The factors most sensitive to the delineation of the groundwater protection zone were identified.
The function f(x) has several parameters. The local sensitivity is the change in the function f(x) caused by a slight change in one of the parameters. The mathematical expression of the local sensitivity analysis is as follows:
S = f ( x 1 , x 2 , x i , , x n ) / x x = x i
where S is the sensitivity value of the function f(x) at the parameter xi. The larger the value of ∣S∣, the more sensitive the parameter xi and thus the greater the influence on the function f(x).
The Morris method is a widely used global sensitivity analysis method, which was proposed in 1991 and improved by Campolonge [27]. The Morris method can quickly screen out sensitive parameters. This method consists of a single random “one change method” experiment, in which the initial value of a given parameter is provided and a certain disturbance (such as 10% or 20%) is then derived based on the initial value, and the influence of changing the value of each parameter is successively evaluated through the combination of different rows of the matrix [28].
The Morris method calculation formula is:
S i ( x 1 , x 2 , · · · , x n , Δ x ) = [ y ( ( x 1 , x 2 , · · · , x i 1 , x i + Δ x , · · · , x n ) y ( x 1 , x 2 , · · · , x n ) ] / Δ x
where Si (X, Δx) is the sensitivity index of parameter i; y(x) is the output result of the model; X = (x1, x2,…, xn) is the n-dimensional vector of parameters; and Δx is the change in x.

3. Results and Discussion

3.1. Particle Trace Characteristics

All actions that could cause direct pollution to the water source are prohibited within the first-grade water source protection zone. Within the second-grade protection zone, it must be ensured that the main pollutants in the area gradually decrease with the flow content in the process of migration to the water source well to ensure that the water quality meets the specified requirements, and sufficient buffer distance must be reserved to deal with sudden water source pollution. According to the requirements in the Technical Specification for the Division of Drinking Water Source Protection Areas, the particle tracking time should be set at 100 and 1000 d when delimiting the first- and second-grade protection zones, respectively. Based on the simulation results of the groundwater flow field and the location of water source wells in the study area, the particle reverse trace line was obtained by using the random walks function in FEFLOW. The simulation results show that the trace length of the first-grade protection zone is 35–523 m, and the trace length of the second-grade protection zone is 50–1800 m.
The hydraulic connection between the groundwater and river changes under the exploitation of riverside source areas, causing the shape of the trace lines to become more complex. In this study, it was found that the trace pattern of the riverside source area can be divided into four types: I exploited water from lateral runoff, II exploited water from the river, III exploited water from lateral runoff and the river, and IV exploited water from lateral runoff, river, and regional runoff on the other side of the river. As shown in Figure 8, the exploited water from the water source wells J06, J14, J15, J23, J24, and J26 come from the lateral runoff of the nearby area (upstream or valley boundary), which belongs to Class I. The water source wells J10, J12, J13, J16, J17, J18, J19, J20, J21, J22, J27, J29, and J30 are located close to the river and the exploited water comes from the river, which belongs to Class II. Water source wells J04, J07, J08, and J32 are exploited from lateral runoff and rivers in the surrounding area, which belongs to Class III. The trace of water source wells J01, J02, J03, J05, J09, J11, J25, J28, and J31 cross the river and belong to Class IV.
The particle trace morphology of 1000 d is similar to that of 100 d, most of which are only the trace length changes, such as wells J08 and J09 (Figure 9). For the exploitation wells with a river supply source, especially the section with a large amount of river leakage, the particle traces of 1000 and 100 d are basically identical, such as the exploitation wells J16–J22 (Figure 9). The trace form can represent the groundwater source of exploitation wells, which has an important guiding role in the delineation of protection zones.

3.2. Protection Area Division and Sensitivity Analysis

3.2.1. Groundwater Protection Zones

When the trace length is less than two times the well spacing, each well is respectively designated as a protection zone by the particle trace range; otherwise, the wells are considered as a well group, and the outermost boundary of all traces is designated as a well group protection zone. Accordingly, the distribution of the groundwater protection zone in the riverside source area is shown in Figure 10. The area of the first-grade protection area of the water source area is designated as 0.375 km2, and the area of the second-grade protection area is 1.20 km2.

3.2.2. Sensitivity Analysis

The permeability coefficient (K) and effective porosity (ne) were selected as parameters of sensitivity analysis in the delineation of the protection area. The optimal value of a group of parameters obtained by manual parameter adjustment was used as the initial value. Under the condition of keeping other parameters unchanged, K and ne were each increased or decreased by 5%, 10%, 20%, and 30% simultaneously, and the particle trace range of 100 and 1000 d of the simulation was delineated. Through field data, it is found that K and ne in the study area are logarithmic, and the fitting formula is n e = 0.027 ln ( K ) + 0.0513 (R2 = 0.98). In about 90% of the study area, the value of K is 100–250 m/d. In this section, K and ne are positively correlated, with a correlation of 0.98. The lithology of the water source is sand gravel, and previous studies have found that K has a positive correlation with ne [29]. Therefore, in the global sensitivity calculation, K and ne increase or decrease at the same time, and thus the numerical simulation conforms to the actual situation.
Figure 11 and Figure 12 show that, under the same recharge and exploitation conditions, the areas of the first-grade and second-grade protection zones are positively correlated with K and negatively correlated with ne. Under the same amplitude, the local sensitivity ∣Sne∣ > ∣SK∣, that is, ne is more sensitive than K. When K and ne change at the same time, the change trend of the area of the protection zone caused by the joint action of K and ne is the same as that caused by ne alone; this indicates that ne plays a dominant role in the delineation of the protection zone. Therefore, the impact of ne should be taken into account in the model construction and the actual delineation of the water source protection zone.
To more intuitively compare the influence of K and ne on the particle trace range, the local sensitivity when the two parameters have a certain amplitude and the global sensitivity when the two parameters change together were calculated. Figure 13 and Figure 14 show that the global sensitivities of the first- and second-grade protection zones calculated by the Morris method are between the local sensitivities of K and ne, and the global sensitivity is closer to the local sensitivity of ne. Moreover, when the values of K and ne are changed together, they restrict each other in the global sensitivity analysis because of their mutual influence. The global sensitivity is not equal to the direct sum of the local sensitivities of the two parameters in terms of value. The change range of the global sensitivity with the change of the parameters is smaller than the local sensitivity, which indicates that the interaction between the parameters reduces the impact of a single parameter on the delimitation of the protection zone. This also shows that the global sensitivity is closer to the actual situation than the local sensitivity in practical application.

4. Conclusions

In this study, the numerical simulation method was used to delineate the groundwater protection zones of large riverside source areas through particle trace lines. The trace length (100 d) of the first-grade protection zone is 35–523 m, and the trace length (1000 d) of the second-grade protection zone is 50–1800 m. The trace pattern of the riverside source can be divided into four categories: exploited water from lateral runoff; exploited water from rivers; exploited water from lateral runoff and rivers; and exploited water from lateral runoff, rivers, and regional runoff on the other side of the river. The particle trace morphology of 1000 d is similar to that of 100 d, but only the trace length changes. The water source protection area is delimited based on the trace form, and the areas of the first-grade and second-grade protection zones are 0.375 and 1.20 km2.
The area of the protection zone was positively correlated with K and negatively correlated with ne. The change trend of the area of the protected area caused by the simultaneous change in K and ne is the same as that when ne acts alone. The global sensitivity is closer to the local sensitivity of ne, indicating that ne plays a dominant role in the delineation of the protection zone. Therefore, the influence of ne should be taken into account in the model construction and the actual delineation of the water source protection zone. At the same time, the global sensitivity is not equal to the direct sum of the local sensitivities of the two parameters, and the range of variation in the global sensitivity with the change of the parameters is smaller than that of the local sensitivity. This indicates that the interaction between the parameters reduces the impact of a single parameter on the delimitation of the protection zone and that the global sensitivity is closer to the actual situation than the local sensitivity.
Numerical simulation is an important method for groundwater protection zone division. The uncertainty of the parameters in the model seriously affects the running results of the model and the accuracy of the model prediction. Different hydrogeological conditions in the study area lead to different parameters affecting the division of the protected zone. Therefore, the precision of the main sensitive parameters should be improved as much as possible in the hydrogeological exploration of water source areas to reduce the impact of parameter uncertainty on the reliability of the division results of water source protection areas.

Author Contributions

Conceptualization, X.L. and J.L.; methodology, W.W.; software, X.L.; investigation, Z.C.; writing—original draft preparation, X.L.; writing—review and editing, W.W.; project administration, J.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Fundamental Research Funds for Central Public Welfare Research Institutes, CAGS (Grant No. SK202009) and China Geological Survey project (Grant No. DD20221752).

Acknowledgments

The authors gratefully acknowledge many important contributions from the researchers of all reports cited in our paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, X.; Wang, W.; Zhang, C. An analysis of the exploitation scheme and groundwater resources of a well field at the riverside based on the Combination of Analytical and numerical Method. China Rural. Water Hydropower 2021, 2, 95–101. [Google Scholar]
  2. Wen, X.; Cheng, Y.; Zhang, J.; Hua, D. Ecological function zoning and protection of groundwater in Asia. J. Groundw. Sci. Eng. 2021, 9, 359–368. [Google Scholar]
  3. Zhu, Y.; Zhai, Y.; Teng, Y.; Wang, G.; Du, Q.; Wang, J.; Yang, G. Water supply safety of riverbank filtration wells under the impact of surface water-groundwater interaction: Evidence from long-term field pumping tests. Sci. Total Environ. 2020, 711, 135141. [Google Scholar] [CrossRef] [PubMed]
  4. Ray, C. Worldwide potential of riverbank filtration. Clean Technol. Environ. Policy 2008, 10, 223–225. [Google Scholar] [CrossRef]
  5. Liu, P. Uncertainty on Numerical Simulation of Groundwater Flow in the Riverside Well Field. J. Jilin Univ. Earth Sci. Ed. 2008, 106, 751–752. [Google Scholar]
  6. Zhou, Y.; Parvez, S.; Nico, M. Analysis of travel time; sources of water and well protection zones with groundwater models. J. Groundw. Sci. Eng. 2015, 3, 363–374. [Google Scholar]
  7. Expósito, J.L.; Esteller, M.V.; Paredes, J.; Rico, C.; Franco, R. Groundwater protection using vulnerability maps and Wellhead protection area (WHPA): A Case Study in Mexico. Water Resour. Manag. 2010, 24, 4219–4236. [Google Scholar] [CrossRef] [Green Version]
  8. Li, J. Establishment and Protection of Drinking Water Source Protection Zone in Germany. Prog. Geogr. 1988, 17, 88–97. [Google Scholar]
  9. Jiang, G.; Teng, M. Partitioning Method of The Groundwater Wellhead Protection Zone. Saf. Environ. Eng. 2016, 23, 36–39. [Google Scholar]
  10. EPA440/5-93-001; Guidelines for Delineation of Wellhead Protection Areas. US Environmental Protection Agency: Washington, DC, USA, 1993.
  11. Fifth Session of the Standing Committee of the Sixth National People’s Congress of the People’s Republic of China. In Water Pollution Prevention and Control Law of the People’s Republic of China; Law Press: Beijing, China, 1984.
  12. Dušan, P.; Dragoljub, B.; Bojan, H.; Dragan, P. Numerical modeling and simulation of the effectiveness of groundwater source protection management plans: Riverbank filtration case study in Serbia. Water 2022, 14, 1993. [Google Scholar]
  13. Zeferino, J.; Paiva, M.; Carvalho, M.D.R.; Carvalho, J.M.; Almeida, C. Long Term Effectiveness of Wellhead protective Areas. Water 2022, 14, 1063. [Google Scholar] [CrossRef]
  14. Amin, M.; Mohammad, N.; Mahdi, L.; Tafreshi, G.M. Determination of the travel time and path of pollution in Iranshahr aquifer by particle-tracking model. SN Appl. Sci. 2019, 1, 1616. [Google Scholar]
  15. Staboultzidis, A.; Dokou, Z.; Karatzas, G. Capture Zone Delineation and Protection Area Mapping in the Aquifer of Agia; Crete; Greece. Environ. Process. 2017, 4, 95–112. [Google Scholar] [CrossRef]
  16. Gurwin, J. Integration of numerical models and geoinformatic techniques in the delimitation of a protection zone for the Study on the characteristics of aquifer system in China. Geologos 2015, 21, 169–177. [Google Scholar] [CrossRef] [Green Version]
  17. Zhang, L.; Su, X.; Meng, X.; Du, S.; Meng, J. Global Sensitivity Analysis of Groundwater Flow Numerical Simulation Parameters Based on the Uniform Design. China Rural. Water Hydropower 2014, 8, 92–97. [Google Scholar]
  18. Tomer, T.; Katyal, D.; Joshi, V. Sensitivity analysis of groundwater vulnerability using DRASTIC method: A case study of National Capital Territory; Delhi; India-ScienceDirect. Groundw. Sustain. Dev. 2019, 9, 100271. [Google Scholar] [CrossRef]
  19. Khonok, A.; Tabrizi, M.; Babazadeh, H. Sensitivity analysis of water quality parameters related to flow changes in regulated rivers. Int. J. Environ. Sci. Technol. 2021, 7, 3001–3014. [Google Scholar] [CrossRef]
  20. Mckay, M.; Beckman, R.; Conover, W. A comparison of three methods for selecting values of input variables in the analysis of output from a Computer Code. J. Comput. Sci. Technol. 1999, 21, 239–245. [Google Scholar]
  21. Dennis, D. The Investigation of Groundwater Contamination in Wicomico County’s Morris Mill Community. J. Environ. Health 2016, 78, 16–25. [Google Scholar]
  22. Uys, L.; Hofmeyr, J.; Rohwer, J. Coupling Kinetic models and Advection-Diffusion equations. 2. Sensitivity analysis of an advection-diffusion-reaction model. J. Silico Plants 2021, 3, diab014. [Google Scholar] [CrossRef]
  23. Karunanidhi, D.; Aravinthasamy, P.; Subramani, T.; Kumar, D.; Setia, R. Investigation of health risks related with multipath entry of groundwater nitrate using Sobol sensitivity indicators in an urban-industrial sector of south India. Environ. Res. 2021, 200, 111726–111738. [Google Scholar] [CrossRef] [PubMed]
  24. Saravanan, R.; Balamurugan, R.; Karthikeyan, M.; Rajkumar, R.; Anuthaman, N.; Gopalakrishnan, A.N. Groundwater modeling and demarcation of groundwater protection zones for Tirupur Basin—A case study. J. Hydro-Environ. Res. 2011, 5, 197–212. [Google Scholar] [CrossRef]
  25. Li, X. Study on Wellhead Protection Area along the Yellow River in Zhengzhou by Numerical Simulation Method. Master’s Thesis, China University of Geosciences, Beijing, China, 2014. [Google Scholar]
  26. Ministry of Environmental Protection of the People’s Republic of China. Technical Guideline for Delineating Source Water Protection Areas: HJ338-2018; China Environmental Science Press: Beijing, China, 2018.
  27. Morris, M. Factorial sampling plans for preliminary computational experiments. Technometrics 1991, 33, 161–174. [Google Scholar] [CrossRef]
  28. Jin, M.; Tian, Y.; Michele, L.; Yu, J.; Zheng, C.-M. Global sensitivity analysis of LID facility model parameters based on Morris; Sobol and EFAST. China Rural Water Hydropower 2022, 6, 104–110. [Google Scholar]
  29. Zhang, H.; Wang, Z.; Wei, Q. Sensitivity analysis of particle tracking trace length and division of groundwater source protection area—An example of Ganhezi Well Field. South-North Water Transf. Water Sci. Technol. 2020, 18, 62–69. [Google Scholar]
Figure 1. Annual average meteorological information from Luanping meteorological station (1973–2018).
Figure 1. Annual average meteorological information from Luanping meteorological station (1973–2018).
Water 14 03702 g001
Figure 2. Water abundance zoning map of the study area.
Figure 2. Water abundance zoning map of the study area.
Water 14 03702 g002
Figure 3. Hydrogeological section perpendicular to the river (Position A–A′ is between Liudaohe Village and Qidaohe Village in Figure 2).
Figure 3. Hydrogeological section perpendicular to the river (Position A–A′ is between Liudaohe Village and Qidaohe Village in Figure 2).
Water 14 03702 g003
Figure 4. Model boundary condition map and Quaternary basement buried depth contour map.
Figure 4. Model boundary condition map and Quaternary basement buried depth contour map.
Water 14 03702 g004
Figure 5. Division of aquifer permeability coefficient and riverbed vertical permeability coefficient.
Figure 5. Division of aquifer permeability coefficient and riverbed vertical permeability coefficient.
Water 14 03702 g005
Figure 6. Flow field fitting diagram for the water source in June 2019.
Figure 6. Flow field fitting diagram for the water source in June 2019.
Water 14 03702 g006
Figure 7. Water level fitting curve of the groundwater dynamic monitoring well.
Figure 7. Water level fitting curve of the groundwater dynamic monitoring well.
Water 14 03702 g007
Figure 8. Particle trace pattern for the 100 d period in each water source well. The red dots are water source wells, and the blue dots form particle traces, indicating the source of water in the exploitation wells. According to the source of groundwater, water source wells can be divided into four categories: (c,i,k) are Class I, (f,h,j,m,o) are Class II, (b,d,q) are Class III, and (a,e,g,l,n,p) are Class IV.
Figure 8. Particle trace pattern for the 100 d period in each water source well. The red dots are water source wells, and the blue dots form particle traces, indicating the source of water in the exploitation wells. According to the source of groundwater, water source wells can be divided into four categories: (c,i,k) are Class I, (f,h,j,m,o) are Class II, (b,d,q) are Class III, and (a,e,g,l,n,p) are Class IV.
Water 14 03702 g008aWater 14 03702 g008bWater 14 03702 g008cWater 14 03702 g008d
Figure 9. Comparison of particle trace patterns of water source wells for 100 and 1000 d, where blue represents the 100 d particle trace and the green represents the 1000 d particle trace. (a) means that when the exploitation well is mainly supplied laterally, the trace morphology of 1000 days and 100 days are basically the same, but the length is increased. (b) means that when the exploitation well is mainly supplied by the river, the trace morphology of 1000 days and 100 days are are basically identical.
Figure 9. Comparison of particle trace patterns of water source wells for 100 and 1000 d, where blue represents the 100 d particle trace and the green represents the 1000 d particle trace. (a) means that when the exploitation well is mainly supplied laterally, the trace morphology of 1000 days and 100 days are basically the same, but the length is increased. (b) means that when the exploitation well is mainly supplied by the river, the trace morphology of 1000 days and 100 days are are basically identical.
Water 14 03702 g009
Figure 10. First-grade and second-grade groundwater source protection areas in the study area.
Figure 10. First-grade and second-grade groundwater source protection areas in the study area.
Water 14 03702 g010
Figure 11. Variation in the area of the first−grade protection zone with K and ne.
Figure 11. Variation in the area of the first−grade protection zone with K and ne.
Water 14 03702 g011
Figure 12. Variation in the area of the second−grade protection zone with K and ne.
Figure 12. Variation in the area of the second−grade protection zone with K and ne.
Water 14 03702 g012
Figure 13. Sensitivity of K and ne to the area of the first−grade protection zone.
Figure 13. Sensitivity of K and ne to the area of the first−grade protection zone.
Water 14 03702 g013
Figure 14. Sensitivity of K and ne to the area of the second−grade protection zone.
Figure 14. Sensitivity of K and ne to the area of the second−grade protection zone.
Water 14 03702 g014
Table 1. Groundwater equilibrium table for the study area in 2019.
Table 1. Groundwater equilibrium table for the study area in 2019.
Equilibrium TermsBudget (104 m3)
Recharge itemsUpstream lateral recharge98.5
River recharge820
Irrigation + Precipitation—Evaporation26
Gully recharge145
Discharge itemsRiver discharge186
Artificial exploitation183
Downstream lateral discharge610
Equalization difference110.5
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, X.; Li, J.; Wang, W.; Cheng, Z. Application of Particle Trace Morphology and Sensitivity Analysis in Delineation of Drinking Water Protection Zone in the Luan River, North China. Water 2022, 14, 3702. https://doi.org/10.3390/w14223702

AMA Style

Li X, Li J, Wang W, Cheng Z. Application of Particle Trace Morphology and Sensitivity Analysis in Delineation of Drinking Water Protection Zone in the Luan River, North China. Water. 2022; 14(22):3702. https://doi.org/10.3390/w14223702

Chicago/Turabian Style

Li, Xiaoyuan, Jianxiu Li, Wenzhong Wang, and Zhongshuang Cheng. 2022. "Application of Particle Trace Morphology and Sensitivity Analysis in Delineation of Drinking Water Protection Zone in the Luan River, North China" Water 14, no. 22: 3702. https://doi.org/10.3390/w14223702

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop