Spatial Assessment of Groundwater Quality Monitoring Wells Using Indicator Kriging and Risk Mapping, Amol-babol Plain, Iran

The main aim of monitoring wells is to assess the conditions of groundwater quality in the aquifer system. An inappropriate distribution of sampling wells could produce insufficient or redundant data concerning groundwater quality. An optimal selection of representative monitoring well locations can be obtained by considering the natural and anthropogenic potential of pollution sources; the hydrogeological setting; and assessment of any existing data regarding monitoring networks. The main objective of this paper was to develop a new approach to identifying areas with a high risk of nitrate pollution for the Amol-Babol Plain, Iran. The indicator kriging method was applied to identify regions with a high probability of nitrate contamination using data obtained from 147 monitoring wells. The US-EPA DRASTIC method was then used in a GIS environment to assess groundwater vulnerability to nitrate contamination, and combined with data concerning the distribution of sources to produce a risk map. In the study area, around 3% of the total area has a strong probability of exceeding the nitrate threshold and a high–moderate risk of pollution, but is not covered adequately by sampling wells. However, the number of monitoring wells could be reduced in most parts of the study area to minimize redundant data and the cost of monitoring.


Introduction
Groundwater is one of the most significant natural sources [1,2], and can be used as an alternative to surface water for drinking, irrigation and industry usage.Poor drinking water quality, high cost of water purification, human health problems, and loss of water supply are attributable to groundwater contamination.The monitoring of the chemical, physical and biological conditions of groundwater is considered to be critical for the planning strategy for the protection of groundwater quality [1].The data obtained from a groundwater monitoring network are valuable for understanding, identifying, and describing modifications in the condition of the groundwater [3,4].A good monitoring network should be indicative of both adequate and appropriate information concerning the groundwater quality as well as be effective in terms of cost [5].Although the amount of information collected from a monitoring network could be increased by more sampling wells, it is costly and probably provides redundant information.Therefore, the optimal monitoring network should provide sufficient data concerning the groundwater quality using the minimum number of monitoring wells [5].Some of the monitoring network designs are inefficient due to the shortage or redundancy of information [6].A new technique can be developed from the probability estimation of the groundwater contaminant concentrations, hydrogeological approaches and evaluation of the pollution risk from anthropogenic activities to assess the groundwater quality monitoring network and evaluate the risky zones of the aquifers.
Geostatistics is a spatial statistical technique that can be used to assess and represent the distribution of concentration over space and time [7].This technique predicts the estimated values based on the relationship between the sample points and estimates the uncertainty of that prediction [8][9][10].Kriging is a linear interpolation procedure that is used to create probabilistic models of uncertainty relating to the values of the attributes.Indicator kriging (IK) is an efficient non-parametric geostatistical method with no assumption regarding the distribution of variables [11], and has the ability to take the data uncertainty into account and predict the conditional probability of certain data for an unsampled location [12].Indicator kriging is also used to identify areas of high probability as potential sites for monitoring based on the current monitoring wells.However, this method alone is not sufficient for the optimal design of monitoring wells, without considering the potential risk from anthropogenic activities and the vulnerable hydrogeological characteristics.The vulnerability of groundwater is characterized by the hydrogeological and geological attributes of the aquifer [13] to specific areas that are more prone to contamination.The DRASTIC model is the most commonly applied vulnerability model based on the physical environmental aquifer parameters to assess groundwater vulnerability [14][15][16][17][18].The existence of potential contamination activities should be considered as a risk for groundwater pollution since the vulnerability only represents the intrinsic characteristic of the aquifer.As the population of an area grows, intensive agricultural activities, inappropriate placement of commercial and industrial regions and high intensity residential areas can potentially cause pollution of the groundwater.Therefore, land use is an additional parameter that can be integrated into the DRASTIC method to evaluate the potential risk in different areas [13,15].
Few researchers have applied the integration of geostatistical techniques and vulnerability assessment as a new approach for redesigning the groundwater monitoring networks [5,[19][20][21].The density of monitoring wells was considered together with vulnerability assessments by Dawoud [22].Yeh et al. [21] applied a genetic algorithm and the factorial kriging method for nine variables-electrical conductivity (EC), total dissolved solids (TDS), Cl − , Na + , Ca 2+ , Mg 2+ , SO 4 2− , Mn, and Fe-for optimal selection of monitoring wells in Pingting Plain, Taiwan.To reach a similar objective, Baalousha [5] developed a new methodology by combining vulnerability and ordinary kriging maps based on the nitrate concentration from groundwater in the Heretaunga Plain, New Zealand.Preziosi et al. [20] developed a GIS based procedure to select the most appropriate monitoring points by combining the actual contamination data, the attributes of water points, and the vulnerability conditions for a variable density network design.
In this study, the data collected from 147 monitoring wells were applied to estimate the potential contamination risk of nitrate in drinking water using indicator kriging on the Amol-Babol Plain, in northern Iran.Moreover, the potential risk zones of groundwater to pollution are specified by integrating the vulnerability and risk mapping in the study area.The main purpose of this study was to develop a new approach to identify areas with high potential pollution and assess the efficiency of the current monitoring wells by probability risk assessment method.

Study Area
The Amol-Babol Plain, which is situated in Mazandaran Province in the northern part of Iran, covers about 1822 km 2 .The plain has a subtropical and humid climate with a hot summer.The annual mean temperature is 17.9 °C [23].The average temperature decreases from the Caspian coastal area towards the Alborz Highlands in the southern region (Figure 1).The annual precipitation is about 880 mm with the majority of precipitation occurring during the rainy season (November and December).
Geomorphologically, the transportation and sediment deposition by the Haraz River, Talar River and Babol River (1b) have evolved into alluvial fans, marine deposits and a flood plain.
The agricultural land, which consists of irrigated fields, dry farming and orchards, covers around 51%, 12.8% and 8% of the Amol-Babol Plain, respectively (Figure 1) [23].Groundwater constitutes 63% of the total water supply in the study area.The annual water abstraction is about 44 million m 3 for domestic consumption and 390 million m 3 for agricultural usage [23].Water in the study area is abstracted from 61,496 shallow wells and 6634 deep wells.Huge amounts of fertilizer, mainly nitrates and phosphates, are applied by farmers to the agricultural land, especially in the first quarter of the year for the agricultural activities (Table 1).
Livestock is the second most important economic activity in the rural parts of the study area.Nitrogenous compounds are the main polluting components present in livestock waste [24].Around 56,885 kg/day of nitrogenous compounds are produced by livestock in the area [23].Runoff, containing nitrate, phosphate and coliforms from the livestock yards or irrigated land may infiltrate the groundwater, especially in areas with a shallow water table (less than 30 m [26], such as the Amol-Babol Plain.

Groundwater Sampling and Analysis
Groundwater nitrate content is highly related to soil type and geological conditions.Generally, the natural value of nitrate in groundwater should only be a few milligrams per liter [27].Inorganic nitrogen in fertilizers, such as nitrate (NO 3 − ), ammonium (NH 4 ) and organic nitrogen in waste, break down to ammonia in the soil and then oxidize into nitrate.Nitrate is essential for plant growth and is used in the synthesis of organic nitrogen compounds [27,28].When not all the nitrate is used by the plants, it may accumulate in the soil and infiltrate into the aquifer.Under anaerobic conditions in the aquifer, nitrate could be completely denitrified or degraded into nitrogen.In order to analyze the nitrate concentration in the groundwater of the Amol-Babol Plain, 147 water samples were collected from the groundwater monitoring wells in 2009.The initial selection of spatial distribution monitoring wells locations was mostly obtained from local experts based on their simple professional experience or empiric-type criteria of the possible sources of nitrate contamination in the study area.The groundwater was pumped out for 10 to 15 min to flush away the non-representative samples of polluted water.The groundwater samples were stored in polyethylene bottles [29] and kept at less than 4 °C in a refrigerator and were analyzed within 24 hours.Nitrate concentration was performed on groundwater samples according to APHA [29] using the colorimetric method.

Vulnerability Assessment
The US Environmental Protection Agency (EPA) developed the DRASTIC model to standardize the methods for evaluating the vulnerability of groundwater to contamination [30].Intrinsic groundwater vulnerability is an index and overlay method, which is dependent on the different hydrogeological parameters of the aquifer system [14].Seven parameters-depth to water table, net recharge, aquifer media, soil media, topography, impact of vadose zone, and hydraulic conductivity (DRASTIC)-are applied by the model for assessing groundwater vulnerability.
A specific rate from 1 to 10 is assigned to each parameter, based on the hydrogeological characteristics (Table 2), with the higher values representing greater vulnerability potential.Then, weights are allocated to each of the seven hydrogeological settings from 1 to 5 (Table 2).[30].The most important parameters have a weight of 5 and the least important have a weight of 1.The DRASTIC index is calculated by multiplying each factor's rating by the assigned weights, as follows:

Table 2. DRASTIC parameters
where D, R, A, S, T, I, and C represent the seven hydrogeological factors, r is the rating and w the weight.The DRASTIC index represents a relative measure of groundwater vulnerability and can range from 26 (very low vulnerability) to 226 (extremely high vulnerability) [31].
Although areas with a low DRASTIC index are less susceptible to pollution compared to areas with a high DRASTIC index, it does not mean that these areas are completely free from groundwater contamination.The DRASTIC model is based on four assumptions: (a) the contaminant is introduced at the ground surface; (b) the contaminant enters the groundwater by precipitation; (c) the contaminant has mobility; and (d) the area should be 400 m 2 or larger [15].

Risk Assessment
The risk of groundwater contamination is not only determined by the intrinsic groundwater vulnerability map, but is also related to the existence of contamination sources.Human activities, which mainly occur at the land surface, could be the source of groundwater contamination.To assess the potential risk of groundwater, an extra factor (land use) should be added to the study (Table 3).The risk map was created based on the following equation [17] where L refers to land use, r to the rating and w to weight.Secunda et al. [32] and Adamat et al. [17] specified three dominant classes of land use activity-built-up area, irrigated field crops and uncultivated land-that could affect groundwater quality (Table 3).Similar to the vulnerability map, the groundwater contamination risk map could be classified into eight classes, from very low (<145) to extremely high (>270) contamination risk (Table 4) [33].

Variogram Analysis
Geostatistics have been defined by Matheron [34] as "the application of probabilistic methods to regionalized variables," indicating that any variable in an area has both random and spatial properties [35].This technique was developed to create mathematical models for a spatial correlation structure [34][35][36] with a variogram that quantifies the spatial variability of random variables between two points [37].The empirical semivariogram, γ (h) , is calculated as half the average quadratic difference between data points separated by the distance vector h [35]: where n (h) is the total number of the variable pairs separated by this distance; and z (x) is the value of the variable.The experimental variogram is fitted into a theoretical model, which comprises eleven different functions: circular, spherical, tetra spherical, pentaspherical, exponential, Gaussian, rational quadratic, Hole effect, K-Bessel, J-Bessel and stable.The cross validation estimation is applied as the goodness-of-fit method for selecting the best variogram model.For accurate prediction, the mean error (ME) and kriged reduced mean squared error (KRMSE) are computed as follows: where Z 0,i is the observed value at location I; Z p,i is the predicted value at the location i; and N is the number of observations and predicted value; S is the standard deviation of the observed value.The corresponding sill (C 0 + C), nugget (C 0 ), and range values of the best-fitting theoretical model are observed.The nugget-sill ratio is utilized in the classification of the spatial dependency of groundwater quality parameters [38].
The variogram can be computed in different directions to detect any anisotropy of the spatial variability.An anisotropic model generally includes geometric anisotropy and zonal anisotropy [39].

Indicator Kriging
Indicator kriging (IK) is applied as a non-parametric geostatistical method to approximate the conditional cumulative distribution function at an unsampled point based on the correlation structure of indicator-transformed data points [11].The indicator kriging function of the observation Z(u) at point u related to the threshold value Z is formulated as follows [11,40]: where Z k is the threshold level.The expected value of I (u;z k ), conditional on n surrounding data, is written as: where ; ( ) ( ) is the conditional cumulative distribution function of Z (u) ≤ Z (k) .Indicator kriging is a form of estimation methodology, in which the method is based on an estimator represented as: where ( ; ) expresses the mount of the indicator at the measured point, U j , j = 1, 2… n, and λ j is a weighting factor of ( ; ) used in estimating ( ; ).
The geostatistical extension module of ArcGIS 9.3 was used for the indicator kriging estimation in this study [10].

Groundwater Vulnerability and Risk Map
In order to assess a vulnerability map, the parameters of the DRASTIC model were prepared using ArcGIS software [10].
-Depth to water (D) The depth to water is one of the significant parameters for determining the distance between the land surface and the groundwater level, through which the contaminants percolate to the groundwater.This layer was obtained using the kriging method for the water depth values measured in 2008.The groundwater depth on the Amol-Babol Plain varies between 2 m and 59.8 m.The shallow wells are mostly located near the coastal area at the northern side and the deep wells are found along the western and southern parts of the study area.Shallow wells have a rating of about 9 (D varies between 1.5 and 4.6) because the groundwater could easily be contaminated by surface runoff and contaminants (Table 2).

-Net Recharge (R)
The net recharge is the amount of water available from precipitation and artificial sources that can penetrate to the groundwater levels [41].The net recharge layer was obtained by calculating the minimum, maximum, and mean value of net recharge using the Guttman Equation [42], as follows: = 0.15 , < 300 = .534 ( − 216), 300 ≤ ≤ 650 = 0.8 ( − 360), > 650 (9) where, R is annual recharge and P is precipitation, both in millimeters.
The average annual rainfall is around 800 mm in 2010.Therefore the net recharge was estimated to be about 350 mm, which is classified in rate 10 (Table 2), for each weather station in the Amol-Babol Plain.The final layer was obtained using the kriging method and the raster map was classified using the rating in the ArcGIS environment.
-Aquifer media (A) Aquifer media refers to the lithology of the saturated zone [43], which was obtained in the study area from the database of the boreholes and wells.The groundwater flow system is influenced by aquifer media, due to the role of rocks or sediment in a formation or group of formations to transfer quantities of water to wells and springs [30].Sandstone, gravel, and silt are the main lithological components of the aquifers from the Amol-Babol Plain.In this layer, the highest rating is 8, which represents the sand and gravel layer (Table 2) along the western side of the study area.The rating is about 6 for the bedded sandstone (Table 2), which is found in the central and along the southern side of the area.The lowest rating is 3, which indicates a clay layer in the limited area of the north and northeastern side of the Amol-Babol Plain.A raster map of the aquifer was obtained by interpolating between the rate values and reclassified into three groups, including the gravel layer, bedded sandstone and clay layer.
-Soil media (S) Soil media forms the uppermost part of the vadose zone, and influences the potential contamination [44].Soil data were driven from the area soil reports [23].The predominant soil in the Amol-Babol Plain was found to be river alluvium.The soil varies gradually from coarse-grained soils, mainly gravel and sand in the highlands, to fine-grained soils, mainly silts and clays, in the coastal area on the northern side of the study area.Therefore, classes 4, 5 and 6 were the identified rates according to silty loam, loam, and sandy loam, respectively, and 9 for sandy soil (Table 2).
-Topography (T) The topography represents the slope of the land surface.This layer was created using digital elevation model (DEM) maps of the Amol-Babol Plain.In low slope areas there is increased infiltration from runoff, which shows greater potential for groundwater contamination, while high slope areas, cannot retain runoff for a long time reducing infiltration to the groundwater.The topography rating changed from 10 for most parts of the plain to 1 for the highlands on the southern side of the study area.

-Impact of the vadose zone
The vadose zone is the ground portion found between the aquifer and the soil cover in which pores or joints are unsaturated, depending on its permeability and the attenuation characteristics of the media [41].The duration of the movement of contaminants into the groundwater is influenced by the zone's components.Lithological and geophysical data from boreholes [23] indicates the presence of silt, clay, shale, and gravel in the vadose zone of the study area.A lower rating (rating 3) is assigned for clay and silt, as they present a less penetrable layer.Larger grain sizes, such as gravels and shales, show higher permeability and lower ability to filter contamination [41].Therefore, a rating of 8 is assigned to the shale and gravel layer (Table 2), which covers the southern part of the Amol-Babol Plain.
-Hydraulic conductivity (C) This parameter represents the ability of the aquifer to transmit water and also controls the rate of contaminant migration from the source to the aquifer [30].Hydraulic conductivity values were calculated based on the transmissibility data from pumping wells and also the thickness of the aquifer from geophysical data.This parameter was computed based on the following equation: where k is the hydraulic conductivity (m/d), T is the transmissivity (m 2 /d), and b represents the thickness of the aquifer (m).
The transmissivity values vary from 100 m 2 /d on the northern and northeastern parts to more than 2000 m 2 /d in the southern and southwestern parts of the Amol-Babol Plain.The aquifer thickness also changes from 25 m in the north to more than 300 m in the west and southwestern side.Therefore, the hydraulic conductivity in the plain varies from 0.04 m/d to 16.4 m/d.The interpolation method was used to create a map, which was classified into rates 1, 2 and 4 (Table 2) to create the conductivity map of the study area.

-Vulnerability mapping
The integration of the seven obtained maps, after multiplying each map with the respective related rating and weights (Equation ( 1)), provides the vulnerability to contamination map for the Amol-Babol Plain (Figure 2a).The minimum and maximum DRASTIC index could be varied between 27 and 240, respectively (Table 4).This range was divided into 8 equal classes, from very low to extremely high vulnerability.The resulting DRASTIC scores from the study area vary between 107 and 169, and were categorized as: 107-119 = moderate low vulnerability (representing 9.4% of the area), 120-139 = moderate vulnerability (representing 72% of the area), 140-159 = moderate high vulnerability (representing 17.06% of the area), and 160-169 = high vulnerability (representing 1.54% of the area).Most of the study area presents moderate and moderate low vulnerability (Figure 2a).The moderate high vulnerable areas occur in the western and northwestern parts, while the highest vulnerability area is located along the western side (Figure 2a).-Risk mapping Land use is an additional parameter that can be integrated into the vulnerability parameters for the potential risk of groundwater contamination.Amol-Babol Plain is mostly covered by irrigated field crops, orchards, marsh, and forests.Therefore, the elements presenting risk are mainly agricultural activities and urban development.The highest weight (5) was used for the land use layer, and a specific rate was assigned to each category (Table 3) based on the significance of the pollution potential contamination for each class.The land use layer for the Amol-Babol Plain shows that agricultural lands and urban areas are the dominant land uses covering 1760 km 2 of the 1822 km 2 total area, and, therefore, were assigned a rating of 8.The resulting risk scores vary between 141 and 211, and the areas are classified as a very low to moderately high risk index area (Table 5).Most of the area corresponds to low to moderate risk zones (Figure 2b).

Groundwater Probability Map
Indicator kriging was applied to create a probability map for nitrate concentration in the groundwater of the Amol-Babol Plain.The measured data at each sampling well were subjected to a continuous scale and converted to a discrete indicator variable with a value of either "1" or "0".Arslan and Demir [2] noted that the difference between variogram results may be due to the weather conditions, irrigation or drainage system of the study area.The best-fitting theoretical models and related semivariogram parameters were chosen to obtain the most accurate estimation, by evaluating eleven different models.The cross-validation was undertaken to determine the difference between the measured and estimated nitrate values in the groundwater, and the spherical semivariogram model was the best-fitting model (Figure 3 and Table 6).The cross-validation, which represents the accuracy of the predictions, shows that the mean error of nitrate in the groundwater is close to 0 (0.002) and that the root mean square standardized is close to one (1.006).17,27].This value was used as the threshold value for human consumption.Concentrations lower than the threshold value were assigned 0, while the higher threshold was assigned as 1.
The potential of nitrate leaching from the soil to the groundwater depends on the amount of rainwater, soil type and well depth.The probability range is classified between very weak (0.0-0.2) and very strong (0.8-1.0) (Table 7 and Figure 4).The probability map of nitrate shows two vulnerable zones-the northwestern and southern sides of the plain (Figure 4)-which cover around 116.2 km 2 (6.37) of the total area.About 12.3 km 2 in the southern part of Babol City was classified as highly vulnerable, as the nitrate concentration in the area exceeded the threshold limit (0.6 to 0.8) (Figure 4 and Table 7).About 103.9 km 2 (5.7%) shows a "moderate" probability of exceeding the threshold value, occurring around the north and northwest area of Amol City, the southern part of Babol City, and the northeastern part of the plain (Figure 4).There is no specific pattern on the area covered by a strong probability (Figure 4), which could be due to the source of nitrate contamination.Domestic sewage and industrial wastewater coupled with an excessive use of fertilizers pose a nitrate hazard to the groundwater of the area.
The amount of infiltration also depends on the groundwater depth.Nitrate concentration is directly related to the extent of irrigation near the well and inversely to the well depth [45].The shallow wells of less than 30 m deep constitute 60.3% of the sampling wells in the Amol-Babol Plain, which are the most vulnerable to nitrate contamination.However, about 1700 km 2 of the total area is suitable for drinking purposes, and the probability of exceeding nitrate concentration threshold is "weak" and "very weak" (Figure 4 and Table 7).

Monitoring Network Assessment
An approach to optimize the groundwater quality monitoring network was obtained using the hydrogeological setting and geostatistics.To achieve this purpose, we used a combination of risk mapping to identify vulnerable areas and probability maps to select areas with insufficient or maybe redundant information.The probability risk map was obtained using a combination of risk and probability maps, which were re-classified into 5 classes (from 1 to 5) to avoid any bias (Figure 5a).The risk map was re-classified from 1 (very low risk) to 5 (moderate high risk) according to Table 5.The probability map of nitrate concentration was also re-classified from 1 (very weak probability) to 5 (very strong probability) based on Table 4.In the probability risk map, the study area was classified from 2 to 7 (Table 8).The highest combined index represents areas with a very strong probability of nitrate concentration and moderate-high risk of contamination (Table 8).The areas with the highest classes should be monitored carefully and more stringently than the areas with lower indexes.However, this index was not observed in the combined map of the Amol-Babol Plain (Table 8).Otherwise, it was found that there is no monitoring well station in the combined index from 6.0 to 8.0 with moderate risk and strong probability areas in the southern part of Babol City and the northwestern side of Amol City.
Moreover, in the area with a combined index between 6.0 and 8.0, in the northwest of Amol City, the distance between two monitoring network wells was found to be more than 3.5 km.An area of 30 km 2 in the northeastern part of the plain with a moderate probability of nitrate concentration and moderate low risk has only one monitoring well.Furthermore, no monitoring well was considered for the northwestern part of the study area.Therefore, the monitoring network should be expanded in the aforementioned areas to provide accurate water quality monitoring and prevent a shortage of information.Additional monitoring wells could be chosen from those agricultural wells that have not been selected for previous monitoring networks (Figure 5b).The extensive areas located in the northern, southern, and eastern parts of the plain have a low-class value with low risk and weak probability of nitrate contamination.Near Babol City and the central part of the plain, the distance between the monitoring wells is less short than 800 m.The number of monitoring wells can be reduced in these areas to minimize redundant information and reduce monitoring costs (Figure 5b).
A monitoring network could be suggested based on 128 sampling wells, for future nitrate assessment in groundwater of the area.This can be achieved by adding 12 new wells in the high risk areas and removing 30 wells from the low risk areas.
Based on the risk map, the results show that part of the southern portion of Babol City was characterized as low risk, while in the probability map, the area was found to have a high probability of nitrate contamination.Hence, when the two maps were compared, it was found that the area has moderate to high probability risk.These results have proven that relying on a risk map alone may be misleading.Therefore, it is advisable to combine both maps to achieve reliable results.
A probability risk map, comprising the risk and probability maps of nitrate concentration from collected samples, is required to evaluate the efficiency of the initial groundwater quality monitoring networks and for redesigning the quality monitoring wells.

Conclusions
This paper developed an approach to assess and delineate areas that require groundwater quality monitoring.The proposed approach assesses the efficiency of the existing network of monitoring wells, considering not only the samples collected from the wells, but also the hydrogeological characteristics of the aquifer and the potential for anthropogenic pollution in the study area.The combination of indicator kriging to estimate the probability of nitrate concentration, together with the vulnerability and risk assessment by the DRASTIC method, represents a powerful and reliable method for identifying optimal sampling locations in the study area.The methodology was designed and successfully applied in the Amol-Babol Plain (Iran).The resultant map of the overlaid factors shows that some areas with strong nitrate probability and high risk of pollution on the south side of Babol City, the northeastern part of the plain and the northwestern part of Amol City are not covered by an adequate number of monitoring wells.However, the majority of the plain, which has a weak probability of nitrate concentration and low risk of contamination, is monitored by multiple sampling wells.Therefore, the existing monitoring wells should be reduced in the lower risk areas and increased in the areas with the highest risk of nitrate contamination.The proposed methodology is general and can be applied to any type of aquifer that is threatened by natural or anthropogenic pollution.In future studies, the proposed method could be applied to assess and redesign the monitoring wells based on various types of pollutants in the aquifer.

Figure 2 .
Figure 2. Groundwater vulnerability map (a); and risk map (b) of the Amol-Babol Plain.

Figure 3 .
Figure 3. Experimental variogram of nitrate concentration and the fitting of theoretical model.

Figure 4 .
Figure 4. Probability map of nitrate concentration in the Amol-Babol Plain.

Figure 5 .
Figure 5. Combined probability map of nitrate concentrations and risk map of pollution (a); Suggested monitoring network (b).

Table 1 .
[25]age of annual nitrate fertilizer usage in the study area[25].

Table 4 .
Classification of vulnerability and risk index values.

Table 5 .
Distribution of risk zones in the Amol-Babol Plain.

Table 6 .
Cross-validation and semivariogram model parameters for probability map of nitrate concentration.

Table 7 .
Probability ranges of area exceeding groundwater nitrate threshold by indicator kriging.

Table 8 .
Classification of probability risk map of nitrate contamination.