Comparative Study of Methods for Delineating the Wellhead Protection Area in an Unconﬁned Coastal Aquifer

: Various delineation methods, ranging from simple analytical solutions to complex numerical models, have been applied for wellhead protection area (WHPA) delineation. Numerical modeling is usually regarded as the most reliable method, but the uncertainty of input parameters has always been an obstacle. This study aims at examining the results from di ﬀ erent WHPA delineation methods and addressing the delineation uncertainty of numerical modeling due to the uncertainty from input parameters. A comparison and uncertainty analysis were performed at two pumping sites—a single well and a wellﬁeld consisting of eight wells in an unconﬁned coastal aquifer in Israel. By appointing numerical modeling as the reference method, a comparison between di ﬀ erent methods showed that a semi-analytical method best ﬁts the reference WHPA, and that analytical solutions produced overestimated WHPAs in unconﬁned aquifers as regional groundwater ﬂow characteristics were neglected. The results from single well and wellﬁeld indicated that interferences between wells are important for WHPA delineation, and thus, that only semi-analytical and numerical modelling are recommended for WHPA delineation at wellﬁelds. Stochastic modeling was employed to analyze the uncertainty of numerical method, and the probabilistic distribution of WHPAs, rather a deterministic protection area, was generated with considering the uncertain input hydrogeological parameters.


Introduction
Worldwide, groundwater provides essential and valuable water resources for drinking water production, agricultural irrigation [1] and industrial processes [2]. Particularly in arid and semi-arid areas where ample surface water sources do not exist, groundwater is often the major, or even the sole, drinking water resource. However, human activities, due to increasing livelihood and development demands, have contributed to the deterioration of groundwater quality. Experience in the last decades has shown that once groundwater is contaminated by chemical, biological or radiological agents, it is always difficult to clean up these contaminants, and the remediation involves high costs [3]. Therefore, protecting groundwater from pollution is a major task for sustainable water management. A central prerequisite for this endeavor is the delineation of wellhead protection area (WHPA) [4].
WHPA delineation is a concept that refers to the adequate management of the areas surrounding groundwater pumping wells where contamination could potentially exist [5,6]. Once WHPA delineation has been conducted, subsequent land-use restrictions and remediation techniques can be applied. WHPA delineation is an effective means of protecting groundwater resources; hence, it is necessary and important to obtain thorough knowledge on how to delineate such an area accurately and precisely. The very early groundwater resource protection with safety area being defined around pumping wells should be dated back to the 1950s in Europe [5], in the purpose of eliminating pathogenic bacteria and viruses from entering groundwater [7]. Later, in the United States, the Amendments to the Safe Drinking Water Act (SDWA) established the first nationwide program to protect groundwater resources [4]. A series of different delineation solutions were then proposed by the U.S. Environmental Protection Agency (EPA) [8]. These methods can be summarized into five major categories, namely, arbitrary and simplified variable shapes, analytical solutions, hydrogeological mapping, semi-analytical method, and numerical modeling [9,10]. Different methods vary greatly in their degrees of complexity and precision, as well as producing different types of protection zones such as the zone of contribution (ZOC), which is the land surface area of the aquifer volume containing all the groundwater that may eventually flow toward the pumping well, or the zone of travel (ZOT), which is defined with a given time of travel (TOT) value for the water to reach the well [8].
Traditional practices determine WHPA using the arbitrary or easily calculated analytical solutions, which ignore the temporal and spatial variability of the local hydrogeological conditions. On the other hand, applying the three-dimensional numerical model gives rise to the possibility of integrating more hydrogeologic characteristics to suit the highly anisotropic and heterogeneous settings [11,12]. The commonly used model in delineating WHPA is MODFLOW [13,14]. Coupled with MODPATH, it provides numerical solution to tracking the pathlines of the particles moving within groundwater flow in a forward or backward manner [15]; thus, the area in which particles could move to can be delineated as the WHPA. However, despite their powerful simulation efficiency, the application of numerical models for WHPA delineation requires a great deal of hydrogeological data from field investigation, through which a complete knowledge of aquifer properties has never been fully achieved. Therefore, the uncertainty of required input data has consistently been an obstacle for obtaining the numerical modeling results [16,17]. Some studies have been conducted to assess the influence of hydrogeological properties on WHPA delineation results, and hydraulic conductivity was identified as the most significant parameter in determining the dimensions of WHPA [18,19]. The principle idea of analyzing the uncertainty of WHPA simulation is to select the input parameters randomly, then run the stochastic modeling over the numerical model to delineate WHPA [6,19].
Reviewing all of the WHPA delineation methods, the simple analytical solution allows us to obtain the protection area through quick calculations, but the results are usually the least accurate, since the real local aquifer conditions are often oversimplified. Numerical modeling has been regarded as the most accurate WHPA delineation method to represent complex aquifer conditions, but the model construction requires the performance of a considerable intensive field investigation and the integration of a large number of hydrogeological parameters, necessitating careful consideration of the modeling uncertainty. Many countries (e.g., the United States, the United Kingdom, Australia, Israel) had listed different methods for defining WHPA, and most of them adopted analytical solutions to define different levels of WHPAs by using ZOTs with different travel times for the national-wide groundwater protection programs [20]. However, a detailed comparative study of using different delineation methods and the comparison between delineation results are still missing in most of the countries. From a practical point of view, the most appropriate method for WHPA delineation should simplify the groundwater flow systems as much as possible for the ease of nationwide program implementation, while still preserving the essential geological and hydrologic characteristics for effective groundwater protection [21].
Therefore, a comparison between the common types of WHPA delineation methods is necessary to ensure a proper protection zone. It is also particularly important for proposing recommendations on the choice of a WHPA delineation method based on the comparison. The major objectives of this study are: (1) to present a comparative study on different methods for delineating WHPAs to examine their performances and to reveal the essential factors that impact WHPA delineation results; (2) to address WHPA delineation uncertainty with numerical modeling by stochastic simulation; and (3) to suggest a general WHPA delineation method selection guideline for conducting groundwater protection program. This study was conducted using two pumping sites in the unconfined coastal aquifer in Israel as an example.

Delineation Methods Selection
Various delineation methods have been introduced and been applied based on different delineation criteria: distance to the water source, drawdown range of pumping well, time of travel (TOT) of contaminant with flow, assimilative capacity of contaminant, and boundaries which control groundwater flow [8,20,22]. Given the fact that delineating WHPA based on contaminant travel time allows defining WHPA into different vulnerability levels, most of the countries incorporated different TOT values into delineation methods to legislate groundwater protection programs [20,[23][24][25][26]. Therefore, in this study, TOT was selected as the delineation criterion. Five commonly used delineation methods were compared: the calculated fixed radius (CFR) method [8], the uniform flow equation method [8,27,28], the WhAEM2000 model [29], the HYBRID method [30] and the MODFLOW-MODPATH model [13,15]. All of the selected methods have been thoroughly explained in the literature, so only brief introductions are provided here. Table 1 summarized the equations for calculating the protection areas by each of the analytical methods.  Uniform flow equation HYBRID The CFR method is based on the assumption that water pumped out from the well comes from the porous aquifer media around the well screen, which is distributed as a cylinder, and the infiltrated water is contributed by precipitation. By using CFR method, a fixed radius around the pumping well will be calculated in terms of a specified TOT value. The uniform flow equation calculates the down-gradient flow boundary X L , the maximum width of the up-gradient zone Y L and the distance to up-gradient boundary r x (with regard to the specified TOT) to delineate the envelope-shaped protection zone. WhAEM2000 is a semi-analytical model developed by U.S. EPA for facilitating capture zone and protection area mapping. WhAEM2000 includes the CFR and uniform flow methods, but also has its own semi-analytical solution to model steady pumping wells and can solve the influence of hydrological boundaries. HYBRID method combines some of the analytical equations from the CFR and uniform flow equation method. By following the five steps illustrated by Paradis and Martel [30], this method is relatively straight forward. MODFLOW-MODPATH is a widely used groundwater numerical model. All of the WHPA delineation analysis was based on the data compiled from past annual average values which covered both the dry and wet seasons.
Application of the selected delineation methods requires that TOT values be specified for different dimensions of protection areas. In Israel, WHPAs are classified into three protection levels, as in many other countries. According to the regulations from the Israel Ministry of Health [31], a Level-A protection zone is a circle closely located to the pumping well with a radius of 10 m, a Level-B protection zone is delineated by a 100-day TOT, and a Level-C protection zone is delineated by a 400-day TOT. The Level-B and Level-C protection zone results from selected methods were compared, since all the Level-A zones are identical.

Reference WHPA
In order to compare the protection zones produced by different methods, a reference WHPA is necessary. We chose the WHPA delineated by MODFLOW-MODPATH as the reference WHPA because numerical modeling most completely accounts for hydrological information of the study area. Referring to the work from Paradis [21] and Miller [32] who tried to assess the differences between WHPAs, we used a comparative index C i to quantify the difference between the reference WHPA and the WHPAs produced by other methods. C i index represents how the WHPA from the tested method fits the WHPA from the reference method. The principle of C i can be illustrated by Figure 1 and C i is calculated by Equation (1): (1) where S CA is the common area between the reference and tested WHPA, S NPA is the non-protected area by the tested method, and S OPA is the over-protected area by the tested method.
Water 2019, 11, x FOR PEER REVIEW 4 of 17 has its own semi-analytical solution to model steady pumping wells and can solve the influence of hydrological boundaries. HYBRID method combines some of the analytical equations from the CFR and uniform flow equation method. By following the five steps illustrated by Paradis and Martel [30], this method is relatively straight forward. MODFLOW-MODPATH is a widely used groundwater numerical model. All of the WHPA delineation analysis was based on the data compiled from past annual average values which covered both the dry and wet seasons. Application of the selected delineation methods requires that TOT values be specified for different dimensions of protection areas. In Israel, WHPAs are classified into three protection levels, as in many other countries. According to the regulations from the Israel Ministry of Health [31], a Level-A protection zone is a circle closely located to the pumping well with a radius of 10 m, a Level-B protection zone is delineated by a 100-day TOT, and a Level-C protection zone is delineated by a 400-day TOT. The Level-B and Level-C protection zone results from selected methods were compared, since all the Level-A zones are identical.

Reference WHPA
In order to compare the protection zones produced by different methods, a reference WHPA is necessary. We chose the WHPA delineated by MODFLOW-MODPATH as the reference WHPA because numerical modeling most completely accounts for hydrological information of the study area. Referring to the work from Paradis [21] and Miller [32] who tried to assess the differences between WHPAs, we used a comparative index Ci to quantify the difference between the reference WHPA and the WHPAs produced by other methods. Ci index represents how the WHPA from the tested method fits the WHPA from the reference method. The principle of Ci can be illustrated by Figure 1 and Ci is calculated by Equation (1): where S CA is the common area between the reference and tested WHPA, S NPA is the non-protected area by the tested method, and S OPA is the over-protected area by the tested method.  [21]).
According to the definition of the comparative index, a higher C i value indicates a better fitting to the WHPA delineated by reference numerical modeling method, and vice versa. According to the definition of the comparative index, a higher C i value indicates a better fitting to the WHPA delineated by reference numerical modeling method, and vice versa.

Stochastic Modeling for Uncertainty Analysis
In order to address the impacts of uncertain input parameters of the MODFLOW-MODPATH numerical model on the results of delineation, an uncertainty analysis was conducted by Latin hypercube sampling (LHS) based stochastic modeling. An uncertainty analysis is an evaluation that takes a set of randomly chosen input values and passes them through the testing model to obtain the distribution of the outputs [33]. LHS, rather the common Monte Carlo sampling, was used to avoid clustering samples and to produce fewer iterations of random input number series. By using LHS method, the probability distribution of the sample parameter firstly being stratified into equal probability intervals on the cumulative probability scale, then a random value is taken from each interval. The number of intervals equals the number of iterations. Hydraulic conductivity (K) was chosen as the randomized parameter as it has been examined by many studies on the importance of determining the size of WHPA and its large uncertainty [6,16,34]. Other parameters, such as aquifer effective porosity, can also affect the WHPA delineation. However, the effects of its spatial variability can be smaller than that of the hydraulic conductivity [19]. By giving the minimum, mean, and maximum values of K, as well as the distribution pattern of K, a series of random K values were produced. The random K values were then introduced into the numerical WHPA delineation model for each study area to get the series of WHPAs with different probabilities.

Study Area
The study was conducted at two pumping sites in Israel: a single groundwater pumping well in the Nahariya city and a multi-well pumping field in the Rehovot area ( Figure 2a). Both pumping sites are municipal drinking water supplying wells located on the unconfined coastal aquifer. The overall relief of the coastal plain aquifer increases gradually from east to west, since a series of mountains are located to the east and in the west is the Mediterranean Sea. The coastal plain aquifer of Israel can be represented by a typical cross-section in Figure 2b [35]. The coastal sediments have a thickness up to 200 m in the west and pinch-out to the east. Calcareous sandstone is the major component of the sediments, and is also the major component of the aquifer storing groundwater [36,37]. Several silt and clay aquicludes exist within the sandstone, which divide the overall aquifer into distinct sub-aquifers, significantly improving the heterogeneity of the local aquifer [38]. Due to the local topography, groundwater moves from the east and is discharged into the sea.

Stochastic Modeling for Uncertainty Analysis
In order to address the impacts of uncertain input parameters of the MODFLOW-MODPATH numerical model on the results of delineation, an uncertainty analysis was conducted by Latin hypercube sampling (LHS) based stochastic modeling. An uncertainty analysis is an evaluation that takes a set of randomly chosen input values and passes them through the testing model to obtain the distribution of the outputs [33]. LHS, rather the common Monte Carlo sampling, was used to avoid clustering samples and to produce fewer iterations of random input number series. By using LHS method, the probability distribution of the sample parameter firstly being stratified into equal probability intervals on the cumulative probability scale, then a random value is taken from each interval. The number of intervals equals the number of iterations. Hydraulic conductivity (K) was chosen as the randomized parameter as it has been examined by many studies on the importance of determining the size of WHPA and its large uncertainty [6,16,34]. Other parameters, such as aquifer effective porosity, can also affect the WHPA delineation. However, the effects of its spatial variability can be smaller than that of the hydraulic conductivity [19]. By giving the minimum, mean, and maximum values of K, as well as the distribution pattern of K, a series of random K values were produced. The random K values were then introduced into the numerical WHPA delineation model for each study area to get the series of WHPAs with different probabilities.

Study Area
The study was conducted at two pumping sites in Israel: a single groundwater pumping well in the Nahariya city and a multi-well pumping field in the Rehovot area ( Figure 2a). Both pumping sites are municipal drinking water supplying wells located on the unconfined coastal aquifer. The overall relief of the coastal plain aquifer increases gradually from east to west, since a series of mountains are located to the east and in the west is the Mediterranean Sea. The coastal plain aquifer of Israel can be represented by a typical cross-section in Figure 2b [35]. The coastal sediments have a thickness up to 200 m in the west and pinch-out to the east. Calcareous sandstone is the major component of the sediments, and is also the major component of the aquifer storing groundwater [36,37]. Several silt and clay aquicludes exist within the sandstone, which divide the overall aquifer into distinct subaquifers, significantly improving the heterogeneity of the local aquifer [38]. Due to the local topography, groundwater moves from the east and is discharged into the sea. The aquifer properties of two sites were compiled from past studies. In the single well pumping site (the Nahariya Site), aquifer recharge due to natural precipitation is approximately 250 mm/year  The aquifer properties of two sites were compiled from past studies. In the single well pumping site (the Nahariya Site), aquifer recharge due to natural precipitation is approximately 250 mm/year [39], local hydraulic gradient was estimated to be ranging between 0.002 and 0.0032, transmissivity is between 775 and 4650 m 2 /day, and hydraulic conductivity is 50-300 m/day. The porosity of the supplying aquifer is estimated as 30% [40]. The Nahariya pumping well has a nearly constant pumping rate of 6600 m 3 /day [41]. In the multi-well pumping field (the Rehovot Site), the local aquifer recharge rate is around 208 mm/year with a hydraulic conductivity of 10.5 m/day [35]. Eight pumping wells in the wellfield are distributed as shown in Figure 2c. The average pumping rate of Well 1 is 3000 m 3 /day, with 1500 m 3 /day of Well 2 and 2500 m 3 /day of Well 6; the remaining wells have the same pumping rate, i.e., 1000 m 3 /day.

Results
For the tested analytical methods (CFR method, uniform flow equation method and HYBRID method), the input parameters are the average values of the data compiled from past field records listed in Section 2. For the semi-analytical model and the reference numerical modeling method, input hydrogeological parameters may have been calibrated to fit the model. Table 2 shows the hydrogeological parameters used by different WHPA delineation methods on the Nahariya and the Rehovot pumping sites. The application of CFR method produced circular protection zones around pumping wells (Figure 3). With a pumping rate of 6600 m 3 /day and the saturated thickness of well is 15.5 m, the radii of Level-B (TOT = 100 days) and Level-C (TOT = 400 days) protection zone was 211 m and 421 m, respectively. The surface area of Level-B zone was 0.14 km 2 and the area of the Level-C zone was 0.56 km 2 . Applying the uniform flow equation method, an average hydraulic gradient of 0.0025 and hydraulic conductivity of 250 m/day [41] were used. The calculated distances from well to the down-gradient boundary (X L ) and the maximum width (Y L ) were 108 m and 341 m. X L and Y L values were the same for both Level-B and Level-C protection zone. The size of WHPA depends on distance to the up-gradient boundary (r x ), which is the travel distance of contaminants with the up-gradient groundwater flow to enter pumping well. r x of Level-B zone was 365 m and was 1098 m for Level-C zone. Thus, two envelope-shaped WHPAs with different dimensions that extend to the up-gradient direction were produced by uniform flow method. The WHPAs delineated by HYBRID method were represented by two ellipses with different horizontal dimensions (d/2) and vertical dimensions (w/2). Surface area of Level-B and Level-C protection zones were 0.14 km 2 and 0.55 km 2 . the ground surface, the hydrogeological WhAEM2000 model was constructed as a simple unconfined sandy aquifer with an average thickness of 15.5 m. The groundwater flow direction and magnitude were assigned by model calibration with the groundwater heads data from three observation wells [41]. As shown in Figure 3, a smaller Level-B WHPA with the surface area of 0.097 km 2 and a Level-C WHPA with the surface area of 0.37 km 2 were delineated. The reference WHPAs delineation was based on the simulation conducted by the Groundwater Modeling System (GMS) 10.2 based on MODFLOW and MODPATH. A 2D single sandy layer model was constructed as shown in Figure 4. The no-flow boundary condition was set for the two boundaries parallel to the groundwater flow, and specific head boundary condition was set for the other two boundaries at the west and east. The same recharge rate of 250 mm/year was used as the analytical methods. WHPAs were outlined with regard to the simulation range of backward particle tracking. The reference Level-B protection zone was 0.0073 km 2 and the Level-C protection zone was 0.29 km 2 . Reference WHPAs aligned with the groundwater flow direction.  Since the aquifer at the Nahariya pumping well is less heterogeneous than that at the Rehovot site, and groundwater was extracted from the substantial sandstone underneath a thin loam layer on the ground surface, the hydrogeological WhAEM2000 model was constructed as a simple unconfined sandy aquifer with an average thickness of 15.5 m. The groundwater flow direction and magnitude were assigned by model calibration with the groundwater heads data from three observation wells [41]. As shown in Figure 3, a smaller Level-B WHPA with the surface area of 0.097 km 2 and a Level-C WHPA with the surface area of 0.37 km 2 were delineated.

WHPA Delineation of the Nahariya Site (Single Well)
The reference WHPAs delineation was based on the simulation conducted by the Groundwater Modeling System (GMS) 10.2 based on MODFLOW and MODPATH. A 2D single sandy layer model was constructed as shown in Figure 4. The no-flow boundary condition was set for the two boundaries parallel to the groundwater flow, and specific head boundary condition was set for the other two boundaries at the west and east. The same recharge rate of 250 mm/year was used as the analytical methods. WHPAs were outlined with regard to the simulation range of backward particle tracking. The reference Level-B protection zone was 0.0073 km 2 and the Level-C protection zone was 0.29 km 2 . Reference WHPAs aligned with the groundwater flow direction. Since the aquifer at the Nahariya pumping well is less heterogeneous than that at the Rehovot site, and groundwater was extracted from the substantial sandstone underneath a thin loam layer on the ground surface, the hydrogeological WhAEM2000 model was constructed as a simple unconfined sandy aquifer with an average thickness of 15.5 m. The groundwater flow direction and magnitude were assigned by model calibration with the groundwater heads data from three observation wells [41]. As shown in Figure 3, a smaller Level-B WHPA with the surface area of 0.097 km 2 and a Level-C WHPA with the surface area of 0.37 km 2 were delineated.  Figure 4. The no-flow boundary condition was set for the two boundaries parallel to the groundwater flow, and specific head boundary condition was set for the other two boundaries at the west and east. The same recharge rate of 250 mm/year was used as the analytical methods. WHPAs were outlined with regard to the simulation range of backward particle tracking. The reference Level-B protection zone was 0.0073 km 2 and the Level-C protection zone was 0.29 km 2 . Reference WHPAs aligned with the groundwater flow direction.

WHPA Delineation of the Rehovot Site (Multi-Well Field)
The same delineation methods were applied for the Rehovot multi-well pumping site. For the application of analytical methods, the WHPA of each pumping well was calculated separately. For the WhAEM2000 and MODFLOW-MODPATH modeling methods, groundwater flow model was constructed for the whole site. Hydrogeological parameter values used for each method were shown in Table 2. The pumping rate of each well in the Rehovot pumping site is much smaller than the pumping rate of the Nahariya well. Consequently, the WHPA of each pumping well would be very small if the TOT values were the same as the Israeli regulation. In order to project all of the Level-B and Level-C WHPAs of each pumping well on the same map to compare the differences, different time scale for delineating Level-B and Level-C protection areas was used for the Rehovot Site. A new TOT of 400 days was used for the Level-B protection zone, and a TOT of 15 years was used for the Level-C protection zone.
The calculation results of the protection area lateral dimensions and surface area of the analytical methods are listed in Table 3. The results of uniform flow equation methods showed great differences between X L , Y L and r x . Therefore, the WHPAs of all pumping wells delineated by the uniform equation method expanded significantly in the direction that was perpendicular to the groundwater flow direction, which already deviated from the realistic protection areas for protecting groundwater resources. Owning to this, the uniform flow equation was not used for delineating WHPAs in the Rehovot wellfield and the surface area of WHPAs were not calculated. Since WhAEM2000 is incapable of simulating the aquicludes within aquifer, a homogeneous aquifer similar to the Nahariya site was built and eight pumping wells were distributed with different saturation thickness and pumping rates.  For the numerical simulations, a three-dimensional model ( Figure 5) of saturated flow and particle transport was used [42]. Considering the more complex hydrogeological conditions at the Rehovot pumping site, five major materials composing the aquifer were defined: sand, silt, sandstone, clay, and gravel as shown by the cross-section A-A' in Figure 5b. The model was calibrated to the observed groundwater levels in observation wells (grey dots in Figure 5a) from 1975-1985 [42], and it was found that a significantly improved fit was obtained when the major water-bearing sandstone was divided into two sub-materials with different hydraulic properties which were located in the western and eastern domain of the model (Figure 5b). The recharge rate (Table 2) was produced through calibration as well and varied in space. The MODPATH module was operated by releasing particles at the groundwater table of each pumping well and running backward tracking. All the WHPAs delineated by different methods can be found in Figure 6. Level-C WHPAs delineated with a longer TOT got overlapped between a few of the wells.
Water 2019, 11, x FOR PEER REVIEW 9 of 17 sandstone was divided into two sub-materials with different hydraulic properties which were located in the western and eastern domain of the model (Figure 5b). The recharge rate (Table 2) was produced through calibration as well and varied in space. The MODPATH module was operated by releasing particles at the groundwater table of each pumping well and running backward tracking. All the WHPAs delineated by different methods can be found in Figure 6. Level-C WHPAs delineated with a longer TOT got overlapped between a few of the wells.

Uncertainty Analysis Results of Numerical Modeling (MODFLOW-MODPATH) Method
Hydraulic conductivities of two pumping sites were randomized to produce a series of conductivity values. For the Nahariya pumping site, the numerical model was built as a homogeneous aquifer; thus only one hydraulic conductivity value, which represents the entire site, was randomized. For the Rehovot pumping site, the hydraulic conductivities of two major sandstone sub-aquifers were randomized. The statistical characteristics of hydraulic conductivities are listed in Table 4. For each pumping site, 64 iterations of parameter sampling were created. sandstone was divided into two sub-materials with different hydraulic properties which were located in the western and eastern domain of the model (Figure 5b). The recharge rate (Table 2) was produced through calibration as well and varied in space. The MODPATH module was operated by releasing particles at the groundwater table of each pumping well and running backward tracking. All the WHPAs delineated by different methods can be found in Figure 6. Level-C WHPAs delineated with a longer TOT got overlapped between a few of the wells.

Uncertainty Analysis Results of Numerical Modeling (MODFLOW-MODPATH) Method
Hydraulic conductivities of two pumping sites were randomized to produce a series of conductivity values. For the Nahariya pumping site, the numerical model was built as a homogeneous aquifer; thus only one hydraulic conductivity value, which represents the entire site, was randomized. For the Rehovot pumping site, the hydraulic conductivities of two major sandstone sub-aquifers were randomized. The statistical characteristics of hydraulic conductivities are listed in Table 4. For each pumping site, 64 iterations of parameter sampling were created.

Uncertainty Analysis Results of Numerical Modeling (MODFLOW-MODPATH) Method
Hydraulic conductivities of two pumping sites were randomized to produce a series of conductivity values. For the Nahariya pumping site, the numerical model was built as a homogeneous aquifer; thus only one hydraulic conductivity value, which represents the entire site, was randomized. For the Rehovot pumping site, the hydraulic conductivities of two major sandstone sub-aquifers were randomized. The statistical characteristics of hydraulic conductivities are listed in Table 4. For each pumping site, 64 iterations of parameter sampling were created. The objective of the stochastic simulation was to generate the stochastic map of the possible areas to be protected and to evaluate the probability distribution of WHPAs that are caused by the uncertainty of input hydraulic conductivities. Here, only the stochastic map produced at the Rehovot multi-well pumping site is shown in Figure 7 since the Nahariya site is a single well and the stochastic modeling results can be represented by the results of the multi-well site. As shown in Figure 7, delineating WHPAs with considering uncertainty of hydraulic conductivity produced a series of protection zones with different possibilities around each pumping well. The darker color area with a higher probability represents a higher priority that this particular location should be protected, and the lighter color refers to a lower probability and a lower protection demand. From Figure 7a,b, it was found that more uncertainties of the generated protection zones, which were indicated by the wider distributed different probability zones, appeared on the southeast of the pumping site. The reason for this was that the groundwater was contributed to by the mountainous area in the southeast side. The heterogeneity of this direction was more influential for the entire groundwater system, as well as for the particle movements. Thus, more uncertainty was encountered in this direction.  The objective of the stochastic simulation was to generate the stochastic map of the possible areas to be protected and to evaluate the probability distribution of WHPAs that are caused by the uncertainty of input hydraulic conductivities. Here, only the stochastic map produced at the Rehovot multi-well pumping site is shown in Figure 7 since the Nahariya site is a single well and the stochastic modeling results can be represented by the results of the multi-well site. As shown in Figure 7, delineating WHPAs with considering uncertainty of hydraulic conductivity produced a series of protection zones with different possibilities around each pumping well. The darker color area with a higher probability represents a higher priority that this particular location should be protected, and the lighter color refers to a lower probability and a lower protection demand. From Figure 7a, b, it was found that more uncertainties of the generated protection zones, which were indicated by the wider distributed different probability zones, appeared on the southeast of the pumping site. The reason for this was that the groundwater was contributed to by the mountainous area in the southeast side. The heterogeneity of this direction was more influential for the entire groundwater system, as well as for the particle movements. Thus, more uncertainty was encountered in this direction. Comparing the uncertainty of Level-B WHPAs (Figure 7a) with the uncertainty of Level-C WHPAs (Figure 7b), more uncertainties were found in the larger protection zone delineation which was indicated by the more expanded protection areas with different probabilities around each well in Figure 7b. This is because confidently determining a larger area is more difficult than determining a relatively smaller area. Spatial variables can increase greatly with expanding area. The red shapes in Figure 7 represent the WHPAs for Rehovot pumping site generated by the deterministic Comparing the uncertainty of Level-B WHPAs (Figure 7a) with the uncertainty of Level-C WHPAs (Figure 7b), more uncertainties were found in the larger protection zone delineation which was indicated by the more expanded protection areas with different probabilities around each well in Figure 7b. This is because confidently determining a larger area is more difficult than determining a relatively smaller area. Spatial variables can increase greatly with expanding area. The red shapes in Figure 7 represent the WHPAs for Rehovot pumping site generated by the deterministic MODFLOW-MODPATH numerical model without considering the uncertainty of the input hydraulic conductivities. By comparing the deterministic results with the stochastic modeling results, it was easy to find that the probabilistic maps basically covered all of the 100% to be protected zones and most of the 80% to be protected zones, which indicates that WHPAs delineated by the deterministic numerical modeling method can provide the most necessary protection for groundwater resources. However, there are still the uncovered protection areas due to the uncertainty of input parameters and this could result in possible risks on groundwater resources contamination.

Comparison of WHPA Delineation Results
The results of WHPA delineation at the two study areas are shown in Figures 3 and 6. In order to compare the application of five different methods, the Level-B and Level-C WHPA results were analyzed in terms of the surface area, geometric property, and the C i value.

Comparison of Level-B Protection Zones
The Level-B protection zones of the Nahariya and the Rehovot pumping site are shown in Figures 3a  and 6a, respectively. WHPAs delineated by different methods were represented by different colors, and the red zones are the reference protection areas delineated by numerical model. The comparative index C i values and the surface area of protection zone of each method are listed in Table 5.
The level-B protection zone delineation results from the two pumping sites showed that the differences among the results were significant when different methods were compared. In general, the protection area from the semi-analytical method-WhAEM2000 model, best fit the reference results in both pumping sites (C i = 58.6% in the Nahariya site, Average C i = 34.1% in the Rehovot site), while the delineation results from the uniform flow equation method (C i = 25.0% in the Nahariya site) and the HYBRID method (C i = 27.0% in the Rehovot site) least fit the reference results. From the perspective of the surface area (S) of the WHPAs being delineated, all of the analytical methods generated overestimated protection zones compared to the surface area of the WHPAs delineated by the MODFLOW-MODPATH model. Exceptionally, the WhAEM2000 model produced a similar scale of protection zones as the reference method.
The application of the CFR method was based on the water balance between the pumping extraction and the aquifer storage and recharge, although the recharge rate can be neglected, as it is usually a small value. In this study, the circular protection zones delineated by the CFR method was four times larger than the reference WHPAs at the Rehovot pumping site, and the C i values at the two pumping sites were relatively small (C i = 43.1% at the Nahariya site, and the average C i = 28.8% at the Rehovot site), mainly due to the unrealistic shape of the WHPAs from the CFR method. The equation of the CFR method suggests that this method is only valid when the supplying aquifer around pumping well can be regarded as a cylinder, which means that the CFR method is most suitable for confined aquifers and unconfined aquifer with a nearly flat gradient (aquifer hydraulic gradient < 0.0015). In the Nahariya site, the well was located in an unconfined aquifer with a hydraulic gradient of 0.0025 (cannot be regarded as a flat water table); thus the reference protection area ought to elongate along the direction of the groundwater flow, which significantly differed from the circular WHPA determined by the CFR method. Similar results can be observed in the Rehovot site. The hydraulic gradient in this wellfield increased gradually from west to eat (from 0.0015-0.0070). Well 1 was located in the west with a nearly flat hydraulic gradient, therefore, the reference WHPA of well 1 tended to be in a similar shape as the WHPA determined by the CFR method. The WHPAs of the other wells in the east had an elongated shape towards the groundwater flow direction as the gradient increased, which can be seen in Figure 6a.
Implementing the uniform flow equation method requires a certain degree of sloping groundwater table (hydraulic gradient). According to the equations of this method, for a nearly flat aquifer with a small hydraulic gradient value, an extremely large width for the up-gradient zone would be produced.
In this study, the uniform flow equation method was only used at the Nahariya site as the pumping well has a relatively large pumping rate (6600 m 3 /day), and it was located in an unconfined aquifer where a significant asymmetrical drawdown could be created, and the difference between X L and Y L values were appropriately small. However, the uniform flow method does not guarantee the water balance between extraction and storage as the CFR method does [30]. Therefore, the C i value (25.0%) of the uniform flow method was still small at the Nahariya site. The HYBRID method integrates the water balance factor and the boundary limits of the protection area. An ellipse-shaped WHPA is usually produced by this method. The C i value of the HYBRID method delineated WHPA was 51.7%, indicating that this method was well suited to the reference method in this setting. However, the HYBRID method performed poorly in delineating WHPAs in the Rehovot site (average C i = 27.0%). Several WHPAs determined by the HYBRID method at the Rehovot wellfield were similar to the WHPAs determined by the CFR method, particularly for the wells located in small hydraulic gradient zones.
The semi-analytical WhAEM2000 model had the highest accuracy (C i = 58.6%) compared with the other analytical methods. At the Rehovot wellfield, although the WhAEM2000 model was still more accurate (C i = 34.1%) than the other methods, the WHPAs from the WhAEM2000 model were much smaller than the reference WHPAs and the WHPAs from other methods. This was because the WhAEM2000 model cannot represent the saturated thickness of an individual well. In the WhAEM2000 model, only a constant thickness (in the Rehovot case, saturated thickness = 85 m) can be given to all the well, which was much greater than the real saturated thickness of each well in the numerical model. Hence, under similar hydraulic conductivity values and hydraulic gradients, the WHPAs determined by the WhAEM2000 method had to be much narrower than the WHPAs determined by the numerical model, in order to ensure that the water flux to the well was equivalent to the pumping rate.

Comparison of Level-C Protection Zones
Level-C protection zones were delineated with a longer TOT value. The Level-C WHPAs of two pumping sites can be found in Figures 3b and 6b. The comparison data are listed in Table 6. Since the Level-C WHPA of the individual wells at the Rehovot wellfield overlapped each other, the protection zone information of the entire wellfield is presented. The comparison results of the two sites showed that the WhAEM2000 model still provided the best fitting WHPAs with the highest C i values, followed by the HYBRID method then the CFR method. With the CFR methods, a much smaller C i value (25.1%) was obtained in the Level-C WHPA at the Nahariya site. This was likely due to the fact that the analytical CFR method does not account for the heterogeneity of local aquifers. When the CFR method was implemented to delineate a relatively small area (Level-B zone), the local heterogeneity could be ignored. However, with the expansion of delineating the higher level protection zone, aquifer heterogeneity could no longer be ignored, and the uncertainty that accumulated with the heterogeneity finally resulted in the deviation of the WHPA. At the Rehovot site, the C i value of the entire Level-C WHPAs was basically the same as the average C i value of the Level-B protection zones. However, CFR method didn't produce the WHPAs that reflect the local groundwater flow condition. For instance, by the reference method, the WHPA of well 1 in the Rehovot site should be nearly circular and not overlapping other protection zones of other wells.
The HYBRID method was still more accurate than CFR method, since it considered the groundwater flow property. The WhAEM2000 model generally fits the reference WHPA well. Despite the fact that the WHPA of well 1 delineated by this model at the Rehovot site was more narrowed and elongated to the hydraulic gradient direction. This is because WhAEM2000 model can only handle the same hydrogeological conditions for all pumping wells, which means that the individual properties of the aquifer where each well was located were disregarded. This further proves that WhAEM2000 model is less appropriate for simulating the multi-well field, though this method provides a relatively accurate approximation and offers an easier alternative compared with building a numerical model.
The WHPAs generated by different methods at the Rehovot multi-well pumping site indicated that WHPAs in the wellfield can be overlapped when there are a few wells nearby (Figure 6), or when delineating WHPAs to a larger range. Comparing the C i values (Tables 5 and 6) of both Level-B and Level-C WHPAs at the Nahariya site with the average C i value at the Rehovot wellfield, smaller C i values were obtained in the Rehovot wellfield, indicating that the same methods performed less accurately in the multi-well field. This is likely due to the influence of well-interference. In a wellfield consisting of several pumping wells, the groundwater flow, as well as the particle transport with the groundwater flow system, can be greatly affected by adjacent wells. Other methods, except WhAEM2000 model and numerical modeling, do not account for the adjacent well interferences.

Recommendation for Selecting Delineation Method
The comparison between WHPAs, determined by different types of methods, showed significant differences. The WhAEM2000 model fits the reference method best, while all the analytical methods overestimated the protection zones. Analysis in detail indicated that it is important to apply the analytical methods to the appropriate hydrogeological setting.
The principle idea of CFR method is to ensure the water balance, so it can only be used for a confined aquifer or an unconfined aquifer with a nearly flat water table (local hydraulic gradient < 0.0015), and it can only provide a very rough estimate for WHPA delineation. The uniform flow equation is only suitable for a sloping aquifer with a hydraulic gradient higher than 0.0015, but this method aligns well to the groundwater flow direction which reflects the regional flow characteristics. The HYBRID method had better fitting compared to other analytical methods, which indicated that water balance and regional flow characteristics are two essential factors for delineating a proper WHPA. Water balance controls the general size of the protection zone, and the consideration of flow characteristics determines the shape of the protection zone. Furthermore, well interferences existed in the multi-well field could have significant impact on the accuracy of WHPAs delineation. Since the analytical methods are not able to consider the influences from adjacent wells, only WhAEM2000 model and numerical modeling methods are recommended. But the WhAEM2000 model is limited to handle the uniform aquifer condition, which means it's not appropriate to be used for WHPA delineation while the aquifer parameters vary greatly at each well location.
The uncertainty analysis through stochastic modeling produced the probabilistic distribution map of the protection zones (Figure 7). This map shows the different necessities of protecting the areas around the pumping well, and improved the deterministic delineation method which excludes the possibly threated zones. It is useful for land-use risk management around pumping wells to prevent potential contaminating activities. A general WHPA delineation method selection flow chart is depicted in Figure 8. This scheme is recommended to choose a suitable protection zone delineation method based on data availability and delineation expectations. interference. In a wellfield consisting of several pumping wells, the groundwater flow, as well as the particle transport with the groundwater flow system, can be greatly affected by adjacent wells. Other methods, except WhAEM2000 model and numerical modeling, do not account for the adjacent well interferences.

Recommendation for Selecting Delineation Method
The comparison between WHPAs, determined by different types of methods, showed significant differences. The WhAEM2000 model fits the reference method best, while all the analytical methods overestimated the protection zones. Analysis in detail indicated that it is important to apply the analytical methods to the appropriate hydrogeological setting.
The principle idea of CFR method is to ensure the water balance, so it can only be used for a confined aquifer or an unconfined aquifer with a nearly flat water table (local hydraulic gradient < 0.0015), and it can only provide a very rough estimate for WHPA delineation. The uniform flow equation is only suitable for a sloping aquifer with a hydraulic gradient higher than 0.0015, but this method aligns well to the groundwater flow direction which reflects the regional flow characteristics. The HYBRID method had better fitting compared to other analytical methods, which indicated that water balance and regional flow characteristics are two essential factors for delineating a proper WHPA. Water balance controls the general size of the protection zone, and the consideration of flow characteristics determines the shape of the protection zone. Furthermore, well interferences existed in the multi-well field could have significant impact on the accuracy of WHPAs delineation. Since the analytical methods are not able to consider the influences from adjacent wells, only WhAEM2000 model and numerical modeling methods are recommended. But the WhAEM2000 model is limited to handle the uniform aquifer condition, which means it's not appropriate to be used for WHPA delineation while the aquifer parameters vary greatly at each well location. The uncertainty analysis through stochastic modeling produced the probabilistic distribution map of the protection zones (Figure 7). This map shows the different necessities of protecting the areas around the pumping well, and improved the deterministic delineation method which excludes the possibly threated zones. It is useful for land-use risk management around pumping wells to prevent potential contaminating activities. A general WHPA delineation method selection flow chart

Conclusions
Since the WHPA cannot be measured but only be calculated or simulated, it is necessary to justify the performance of the commonly applied WHPA delineation methods. In this comparative study, five different WHPA delineation methods were employed to define the protection zones for two pumping sites in the unconfined coastal aquifer in Israel. The MODFLOW-MODPATH numerical modeling method was selected as the reference method. By comparing the WHPA results with the reference WHPA, the semi-analytical WhAEM2000 model provides the best fitting, while all of the analytical methods overestimated the protection zones. The CFR method, which is based on water balance between aquifer extraction and storage, provided fast WHPA approximation. But it can only be used for the confined or the unconfined aquifer with a nearly flat water table. On the other hand, the uniform flow equation method preserved groundwater flow property and can only be applied at sloping aquifers. However, the latter does not consider the water balance as the CFR method does; therefore, it could produce unrealistic WHPAs. The better fitting to the reference method of the analytical HYBIRD method indicated that it is essential to ensure water balance calculation as well as to reflecting regional flow characteristics while defining the WHPA. Dealing with a complex hydrological condition and a multi-well field, only semi-analytical or numerical modeling should be applied. Stochastic modeling of both pumping sites generated a probability distribution map for WHPAs, which suggested the different necessities of protecting various areas around the pumping well and offered risk management into WHPA delineation. Therefore, selecting a WHPA delineation method should be based on data availability and delineation expectations. The recommended procedure for choosing the proper delineation method could facilitate this process and resolve the puzzle of finding the "best" method.
Author Contributions: This paper was completed under the supervision of A.Y. and N.W.; Y.L. conducted the calculation and modeling, analyzed the results and wrote the paper. N.W. contributed to the analysis discussion and helped with reviewing the manuscript. A.Y. guided building the concepts of this research, helped with analyzing the results and building up the numerical models, as well as reviewing this manuscript.
Funding: This research was funded by Maccabi Carasso, grant number 87382911.