Next Article in Journal
Removal of Manganese(II) from Acid Mine Wastewater: A Review of the Challenges and Opportunities with Special Emphasis on Mn-Oxidizing Bacteria and Microalgae
Previous Article in Journal
Gendered Water Insecurity: A Structural Equation Approach for Female Headed Households in South Africa
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modified Index-Overlay Method to Assess Spatial–Temporal Variations of Groundwater Vulnerability and Groundwater Contamination Risk in Areas with Variable Activities of Agriculture Developments

1
Graduate Institute of Applied Geology, National Central University, Taoyuan 32000, Taiwan
2
Department of Geology, University of Science Ho Chi Minh City, Ho Chi Minh City 700000, Vietnam
*
Author to whom correspondence should be addressed.
Water 2019, 11(12), 2492; https://doi.org/10.3390/w11122492
Submission received: 28 October 2019 / Revised: 21 November 2019 / Accepted: 23 November 2019 / Published: 26 November 2019
(This article belongs to the Section Hydrology)

Abstract

:
The groundwater vulnerability (GV) assessment for contamination is an effective technique for the planning, policy, and decision-making, as well as for sustainable groundwater resource protection and management. The GV depends strongly on local hydrogeological settings and land-use conditions that may vary in response to the activities of agricultural development. In this study, a modified DRASTIC model, which employs an additional factor of land use coupled with the analytic hierarchy process (AHP) theory, was used to quantify the spatial and temporal variation of GV and groundwater contamination risk in the Pingtung groundwater basin. The results show that the GV slightly decreased due to the decrease in agricultural areas under the change of land use over two decades (1995–2017). The yearly changes or a shorter period of observations incorporated with the accurate land-use map in DRASTIC parameters could improve GV maps to obtain a better representation of site-specific conditions. Meanwhile, the maps of yearly contamination risk indicated that the counties of Jiuru and Ligang are at high risk of nitrate pollution since 2016. In other agriculture-dominated regions such as Yanpu, Changzhi, and Gaoshu in the Pingtung groundwater basin, the climate conditions influence less the temporal variations of groundwater contamination risk. The results of this study are expected to support policy-makers to adopt the strategies of sustainable development for groundwater resources in local areas.

1. Introduction

Groundwater is a crucial water resource for domestic use in arid and semi-arid regions because of its relatively stable and large storage capacity, easy extraction, and low vulnerability to contamination as compared to surface water resources [1]. However, due to the growth of population, leading to an increase in agricultural activities, sewage disposal, and industrial wastewater, the pollution of groundwater is now a vital issue for groundwater resource sustainability [2,3]. Contaminated groundwater treatment is costly, and the prevention of groundwater pollution is, thus, an effective solution to protect groundwater resources [4]. The assessment of groundwater vulnerability (GV) is recognized as an efficient approach to determine the potential zones of contamination linked to natural and social–economic conditions [5].
The concept of GV was proposed for decades, and the vulnerability maps are critical references in decision-making for land-use planning, groundwater management, monitoring, and remediation [6]. Previous investigations proposed three basic approaches to evaluate GV for specific groundwater basins. These approaches can be categorized as (1) index-overlay methods, (2) process-based physical methods, and (3) statistical methods [4]. Among these approaches, the index-overlay methods are recognized to be efficient and useful to identify GV and are especially proper for utilization with a geographical information system (GIS), which provides useful tools to overlay and integrate the different multiple maps [2,3,7,8].
The DRASTIC model is the most efficient and widely used index-overlay method for quantifying GV because of its ease of use, minimal data demand, and good definition of GV. With the development of GIS and the associated database, the simple calculations of the DRASTIC model are physically meaningful and well-suited for large-scale problems. Investigators took into account the correlation between the vulnerability index and contaminant concentration in combination with the sensitivity analysis of DRASTIC parameters for vulnerability assessment [7,8,9,10]. However, the GV maps obtained from the DRASTIC model might be unreliable because the feature ratings and weights in the model are site-specific, and the accuracy highly depends on the experience of the investigators. Recently, many studies proposed the modification of the conventional DRASTIC model for mapping groundwater vulnerability zones, including (i) geostatistics [11,12], (ii) an artificial neural network (ANN) [13,14], and (iii) the analytical hierarchy process (AHP) [7,8,15]. Among them, the AHP revealed that the modification based on the local hydrogeological setting shows better and more accurate results. Optimizing the weights of DRASTIC using AHP was more effective [15]. Additionally, the DRASTIC model with the AHP also shows a high correlation between the GV map and the contaminants as compared to field measurements [16].
The DRASTIC method is flexible, and it can be integrated with other parameters such as lineament, land use/land cover, aquifer thickness, and impact of contaminant, depending on the local hydrogeological conditions and data availability in specific areas [2,7,17]. Recent studies considered land use in a DRASTIC model to improve GV maps [18,19,20]. In general, the change in land use might influence the net recharge in the study areas [21]. Land use may directly influence the groundwater recharge process, and it plays a critical role in the prediction of the impacts of human activities on groundwater quality. Ribeiro [22] assessed the impact of the land use on groundwater pollution potential for an agricultural coastal plain using the susceptibility index (SI) method. Lima et al. [17] applied a modified DRASTIC model and Dynamic Conversion of Land-Use (Dyna-CLUE) model to predict the effect of future agricultural expansion on GV. The factors in GV models typically use fixed groundwater levels and average precipitation to quantify the impact of land use on the variations of GV maps. Li and Merchant [23] further analyzed the impact of projected land-use change on the vulnerability of groundwater to pollution. They considered that DRASTIC factors such as depth to water, recharge, and land-use conditions might change in response to future changes in climate and socio-economic conditions. The model was used to evaluate groundwater pollution risks under future climate and land-use changes in North Dakota. Results showed that high GV would expand northward and northwestward in eastern North Dakota under different scenarios. Unlike previous studies, this study is the first to consider temporal variations of GV that are connected to future climate and socio-economic conditions.
The critical implementation of the GV map is to quantify the water-related health risks that are induced by different anthropogenic activities. Groundwater pollution risk is defined as the process of estimating the possibility that a particular event may occur under a given set of circumstances [24]. The groundwater risk assessment integrates information from the aquifer vulnerability, the pollution loads, and the probability of contaminant occurrence. With recent developments in the risk analysis models, Kazakis and Voudouris [25] employed the DRASTIC model and a nitrogen loss method to assess the GV and pollution risk in porous aquifers. They replaced qualitative parameters such as aquifer type, soil, and impact of the vadose zone with quantitative parameters, including aquifer thickness, nitrogen losses from the soil, and hydraulic resistance. Results of the modified DRASTIC models enabled the estimation of the pollution risk to nitrates. Goudarzi et al. [26] combined the DRASTIC and Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) method to map the zones of groundwater pollution risk caused by agricultural practices. The results from their study efficiently identified the area with high-risk levels, which might need emergency recovery action plans. However, the data without long-term observations of climate and land-use data limit the implementation of GV and risk maps to account for temporal variations in a specific groundwater basin. The spatial and temporal evolution of the health risk to groundwater pollution is critical to support decision-making for local development plans.
Based on long-term observations of groundwater quality data and available data of land use in Taiwan, this study aimed to quantify the spatial and temporal evolution of GVs and zones of groundwater contamination in the Pingtung groundwater basin in southern Taiwan. More specifically, the DRASTIC model, coupled with the factor of land use and the AHP theory, was employed to calculate the variations of GV. The results of GV were validated with site-specific nitrate concentration values obtained from 2010–2017. The comparison of GV for different land-use maps aimed to quantify the influence of land use on the results of different DRASTIC models. The year-averaged risk maps were then obtained by integrating the validated GV maps, the pollution severity interpolated by ordinary kriging (OK), and the probability of occurrence calculated by indicator kriging (IK). Additionally, the study conducted sensitivity analysis to determine significant input variables that influence the reliability of GV mapping. The results of this study are expected to support policy-makers to adopt the strategies of sustainable development for groundwater resources such as the mitigation of groundwater pollution, the planning of land-use patterns and practices, and the allocation of groundwater resources in local areas.

2. Materials and Methods

2.1. Study Area

The Pingtung plain groundwater basin has a total area of 1210 km2, and it is located in the southwestern part of Taiwan (see Figure 1). It is bound by the low hill of the Quaternary sediments in the west, foothills and river valleys to the north, the Fengshan fault to the west, the Taiwan Strait to the south, and the Chaozhou fault to the east. The plain mainly comprises flat areas of Kaoping, Tungkang, and Linpien river catchments. These rivers form not only the natural discharge base for regional groundwater in the basin but also the largest drainage area in Taiwan [27]. The average annual rainfall and evaporation are about 2493 mm and 1723 mm per year, respectively. However, the distribution of precipitation is uneven in that the ratio of rain in the dry seasons (October to April) to that in the wet season (May to September) is 1:9 [28].
There are two unconsolidated deposits of the Late Pleistocene and the Holocene, which constitute the principal aquifer in the groundwater basin. According to the hydrogeological investigation from the Taiwan Central Geological Survey (CGS) [28] in 2002, the Pingtung plain is divided generally into a distal fan and a proximal fan. Subsurface hydrogeological analysis within a depth of 250 m was conducted between 1995 and 1998, which partitioned the plain sediments in the distal fan into eight overlapping sequences, comprising four marine sequences and four non-marine sequences [28]. Non-marine sequences with highly permeable coarse debris, arranged from medium sand to gravel, are identified as the aquifers in the basin. Meanwhile, marine sequences are composed of fine debris, classified from clay to fine sand, and they are considered as aquitards (see Figure 2). The aquitards are mainly observed in the distal fan. The aquifers are highly permeable with an averaged hydraulic conductivity greater than 100 m/day. There is no aquitard observed in the proximal fan. In the groundwater basin, the aquifer recharge is mainly from the river valleys and the proximal fan in the eastern and northern study areas.
The Pingtung plain is a crucial area for agricultural production in Taiwan, with a significant proportion of the area used for agriculture. It is not common that only relatively 45.8% of people of the Pingtung plain utilize clean groundwater that was treated at a water plant (about 92.93% of the population in Taiwan uses tap water) [11]. The substantial amount of groundwater in the Pingtung plain still needs to meet drinking and household demands. Recently, nitrate contamination was determined as a significant non-point source of groundwater contamination in Taiwan. The report of elevated levels of NO3–N in groundwater were published by Agriculture Engineering Research Center [29] based on long-term groundwater quality monitoring. Therefore, the assessment of GV is a useful approach for creating the strategy of groundwater quality protection and the development of groundwater resource sustainability.

2.2. DRASTIC Method

In the study, the DRASTIC model linked to GIS was used for assessing the vulnerability of the study area. As described by Aller et al. [1], the DRASTIC model is a composite description of all the major geological and hydrological factors that affect and control the groundwater movement into, through, and out of an area. The conventional DRASTIC model uses seven hydrogeological parameters to assess the GV, namely, depth to groundwater table (D), net recharge (R), aquifer media (A), soil media (S), topography (T), the impact of vadose zone (I), and hydraulic conductivity (C). The DRASTIC vulnerability index (DVI) is a weighted sum of the seven parameters, as shown in Equation (1).
D V I = D r D w + R r R w + A r A w + S r S w + T r T w + I r I w + C r C w ,
where D, R, A, S, T, I, and C represent the seven hydrogeological parameters. Subscriptions r and w denote the corresponding ratings and weightings for the parameters.
The minimum and maximum DVI values of the different vulnerability maps usually vary in different ranges. Therefore, the study normalized the DVI values for comparison purposes. The advantage of taking normalization is to break down the relationships of anomalies to generate smaller, well-structured relationships among different DRASTIC maps. The value of DVI maps was converted to the range of 0 to 100 using the following formula:
X = X X min X max X min × 100 ,
where X’ is the normalized index value, and X is the original index value obtained from different DRASTIC models. In the study, we added land use (LU) as the additional parameter for different DRASTIC models and compared the results influenced by the change of land use. The calculations of GV followed the same concept and procedures.

2.3. Analytical Hierarchy Process (AHP)

The AHP is one of the most commonly used multi-criteria evaluation methods developed by Saaty [31] for decision-making when confronted with conflicting and qualitative criteria. AHP constructs a series of multiple pair-wise comparisons based on a standardized comparison scale of nine levels. A comparison is performed to quantify a rating or weighting for each criterion using a number from 1–9 (extremely less important to extremely more influential). The effect of the factors or criteria for the decision-making process is quantified based on the scale combined with the experience and knowledge of specialists or users [32]. AHP can capture both subjective and objective evaluation measures and test the stability of the evaluation methods and options proposed by specialists or decision-makers, thereby decreasing the errors in decision-making [32]. For the cases in the study, there are significant steps to solve a decision-making problem using AHP [31], shown below.
Step 1: Decompose the issue into the hierarchy of elements, namely, criteria and sub-criteria. In the study, the standards were the eight parameters (i.e., D, R, A, S, T, I, C, and LU) in the DRASTIC model and the sub-criteria were the ranges within each setting.
Step 2: Develop (n × n) judgment matrices (B) by pair-wise comparison on the n criteria, in which all the diagonal elements are one as the same criterion. Afterward, the relative weight for each matrix is identified using the right eigenvector (w) that agrees with the largest eigenvalue (λmax).
B w = λ max w ,
where, as λmax approaches n, the judgments in the pair-wise comparison matrix become more consistent. Therefore, the difference, λmaxn, can be used as a test of inconsistency.
Step 3: Check the consistency of the judgment. Saaty [30] defined a consistency index (CI) as the relationship between the admission of B matrices, i.e., bij × bjk = bik. The formula can be described by
C I = λ max n n 1 .
The calculation of the consistency index is based on the random pairwise comparison applied to matrices with different sizes. Taking the average of CIs for each matrix size allows obtaining the random index (RI). The ratio of CI and RI is then defined as the consistency ratio (CR).
C R = C I R I .
If the CR value is greater than 0.1, it means that there is severe inconsistency in the judgments while creating the pairwise comparison matrix. Hence, CR ≤ 0.1 should be kept for maintaining the steadiness of the array. In this study, the CR was 0.017 for weight estimation and 0.025 for rating calculation, indicating that a logical matrix was created for the judgement. Table 1 shows the consistency ratios of the selected parameters in this study. The CR values for different parameters were less than 0.1.

2.4. Groundwater Contamination Risk

Groundwater contamination risk can be defined as the probability that groundwater in the aquifer will become contaminated to an unacceptable level due to activities on the immediately overlaying land surface. In other words, the risk is the combination of the occurrence probability and consequences of certain hazardous events. According to Morris and Foster [33], they defined the risk of aquifer contamination as the amount of damage and the probability of the contaminants exceeding the acceptable threshold caused by agricultural and human activities. These factors yield the following relationship for estimating pollution risk:
Risk = probability of an event × consequential damage.
The consequential damage includes the factors of GV and groundwater contamination. Probability is the likelihood of the occurrence of a contaminant which damages the aquifer. In this study, the criterion considered for the damage of the aquifer was the standard of irrigation water quality. Based on the proposed concept, the calculation of groundwater contamination risk in Equation (6) gives the following formula:
Risk = Vulnerability × groundwater pollution × probability.
To obtain the risk maps for the study area, we used different interpolation methods that collected data at well locations and estimated preferred variables at points without observations. The study employed the ordinary kriging (OK) method for the interpolation of the nitrate concentration in the Pingtung plain groundwater basin. For the calculation of the probability of exceeding the tolerable level of 10 mg/L, we conducted indicator kriging (IK) interpolation based on the specified threshold. Note that the tolerable nitrate concentration is the irrigation water quality standard in Taiwan. Other target contaminants and tolerable standards can be used for cases different from the one used in the study. The kriging methods are geostatistical interpolation methods that follow the concept of best linear unbiased estimation. In the OK and IK, the required geostatistical parameters must be determined based on spatial analysis of observation data, i.e., the nitrate concentration at monitoring wells in the Pingtung plain groundwater basin. In the study, the selection of the OK variogram model relied on the best fit of the experimental variograms obtained from the data. We also checked the suitability of interpolated concentration by using the Pearson coefficient and spatial stratification in the geographical detector [34].

2.5. The Research Flowchart

Figure 3 shows the flowchart of the study. To obtain the variations of GV and risk maps for the study area, the conventional DRASTIC model used the original seven parameters. The modified DRASTIC model included the land-use maps in 1995 and 2017 as the additional parameter. In the modified AHP/DRASTIC model, the study employed the AHP theory to modify the weightings and ratings in the modified DRASTIC model. Validations of the three different DRASTIC models used the data of site-specific nitrate concentration values taken from 71 monitoring wells in the Pingtung plain groundwater basin. The effect of the land-use parameter on the GV estimation was then quantified based on the change of land use in the past two decades. Coupled with the climate and the associated recharge estimations, we aimed to assess the variations of GV in the period of 2010–2017. The generation of groundwater risk maps relied on the integration of validated GV, the pollution severity, and the probability of occurrence. In the study, the preparation of the pollution severity employed the ordinary kriging (OK) interpolation, while the distribution of the probability of occurrence used the indicator kriging (IK).

2.6. Sensitivity Analysis

The advantage of sensitivity analysis is to provide beneficial facts in terms of the effect of weights and ratings given to each factor, and it helps experts in hiding the importance of subjective elements [35]. In most cases, the sensitivity analysis can be conducted based on removal sensitivity analysis introduced by Lodwick et al. [36] and a single-parameter sensitivity analysis proposed by Napolitano and Fabbri [37].
The map removal sensitivity analysis describes the sensitivity of the vulnerability index when removing one or more parameters from the suitability analysis [36]. The evaluation relies on the calculation of the variation index (%), as shown in Equation (8).
S = V i N V x i n V i × 100 ,
where S, Vi, and Vxi are the sensitivity index, total vulnerability index, and the vulnerability index excluding one map layer, respectively. Notation N denotes the number of map layers used for computing Vi, and n represents the number of maps used for calculating Vxi.
The study also employed the single-parameter sensitivity analysis [37] to assess the influence of each parameter on the vulnerability index, according to the methods used for the GVs. Therefore, the effective weight of each parameter was determined for the final GV map. The calculation of the effective weighting for each i-th sub-area used the following formula:
W x i = P r i × P w i V i × 100 ,
where Wxi is the effective weight of each parameter, and Pri and Pwi denote the rating and weight of the parameter (x) in the i-th sub-area. Notation Vi is the overall vulnerability index.

2.7. Data Preparation and Processing

In the study, there were seven fundamental parameters for the conventional DRASTIC model. We prepared the parameter maps using ArcGIS software. The modified DRASTIC considered the additional land-use factor to estimate the GV. Different types of datasets were used to create the ratings and weights of the seven factors for the appropriate DRASTIC methods. All data layers created for the study area were georeferenced within a GIS environment using the TWD_1997_TM_Taiwan projection system. All vector data were converted into a raster format with a cell size (pixel) of 1000 × 1000 m. Table 2 shows the sources of data used for creating different GVs. Below, we briefly describe the data collection and the data processing used for different DRASTIC models.

2.7.1. Depth to Water (D)

This factor is crucial because the term “depth to water” determines the potential distance of contaminant transport through the infiltration process. The D is defined as the distance (in meters) from the land surface to the groundwater table at a specified point. In the study, we collected the groundwater levels at 47 monitoring wells from 2010 to 2017 and used the averaged groundwater levels to prepare the groundwater depth map (Figure 4). There were two types of D prepared for the study. The long-term averaged D values accounted for the averaged groundwater levels in the period of 2010 to 2017. However, the yearly averaged groundwater levels in each year were used for assessing spatial and temporal variations of GV and pollution risks.

2.7.2. Net Recharge (R)

The net recharge is the total amount of annual precipitation that percolates from the land surface to the aquifer system. The percolation enables the contaminant to transport vertically to the groundwater table and allows it to be horizontally distributed in the aquifer [1]. The study determined the factor R using the water balance equations described by Bisson and Lehr [38].
R = P E T R s f ,
R s f = ( P 0.3 S ) 2 ( P + 0.7 S ) ,
I a = 0.3 S ,
S = 25400 C N 254 ,
where W, P, ET, and Rsf denote the net recharge (mm), precipitation (mm), evapotranspiration (mm), and surface runoff (mm). S, Ia, and CN present the maximum potential retention (mm), initial abstraction (mm), and curve number (dimensionless). The observation data of monthly precipitation and evapotranspiration were obtained from the Central Weather Bureau of Taiwan for 2010–2017. The runoff curve number (CN) method considers the rainfall, soil types, and land uses (or land covers) to calculate the potential runoff (Rsf). In this study, we obtained CN for different hydrologic soil groups from Hong and Adler [39], and used the land-use map of 2017 for the recharge estimation. Therefore, the distribution of factor R could be obtained based on the concept of the water balance for the aquifer system [38]. Note that the land uses and land covers represent different sources of data. The effect of using land use or land cover on the results of GV maps is discussed later. The factor R had a distribution pattern similar to the regional precipitation, indicating that the regional precipitation dominates the evaluation of net recharge in the study area (Figure 4).

2.7.3. Aquifer Media (A)

This parameter characterizes the geologic composition of the aquifer in the top layer, which was obtained using a hydrogeological map, geological cross-sections, and borehole profiles. The top layer of the aquifer controls the pollutant attenuation processes through the saturated zone material properties. The types of composition have a distinct range of permeability of the aquifer. In the study, the subsurface geology map, geological sections, and 47 drilling profiles were used to define the lithological properties for the aquifer media in the DRASTIC models (Figure 4). With the available data from the Pingtung plain groundwater basin, there were four types of alluvium aquifer material classified for the GVs (Table 3). The map of aquifer media shows that the low-hydraulic-conductivity area appears in the distal fan. The river confluence flow zone of the Kaoping river in the mid-fan also obtained a relatively low hydraulic conductivity in the top layer of the aquifer.

2.7.4. Soil Media (S)

The soil media plays a significant role in the movement of groundwater contamination from the surface. The soil media refers to the topmost portion of the vadose zone, particularly described by substantial biological activities [40]. The spread of pollutants to the aquifer can be reduced, delayed, or accelerated by fine-textured materials such as silts and clays. The distribution of soil taxonomy was obtained from digitalized soil maps of the Pingtung plain groundwater basin. The soil media could be classified into three categories, including sandy loam, loam, and clay loam (Figure 4).

2.7.5. Topography (T)

In the DRASTIC model, the evaluation of topography relies on calculating the slope variability of the land surface. The slope can be the indicator to represent the detention time of the contaminant on the land surface. The concept implies that the surface runoff on the land can carry the contaminant with the surface water. Therefore, a steep slope can lead to relatively high surface runoff and results in a low vulnerability of groundwater contamination [2,40]. In the study, a digital elevation model (DEM) with a resolution of 20 m was used, and the slope percentage was calculated for the parameter T (Figure 4). We classified the slope into five different categories (Table 3).

2.7.6. Impact of Vadose Zone (I)

The vadose zone is defined as the unsaturated layer above the groundwater table, which is an essential parameter in the assessment of vulnerability because it affects the time available for the propagation of pollutants to the aquifers. Thus, the type of vadose zone media determines the attenuation characteristics of the material below the typical soil layer and above the water table through some bio-geochemical processes like bio-degradation, mechanical filtration, sorption, volatilization, and dispersion [1]. The data of the vadose zone were collected from the 47 well logs provided by Taiwan CGS [28]. Silt to clay, fine to coarse sand, and fine to coarse gravel were the major classes of the vadose zone in the study area (Figure 4).

2.7.7. Hydraulic Conductivity (C)

Hydraulic conductivity is necessary because it refers to the ability of the aquifer material to transmit water, and it controls the rate contaminant movement [1]. In general, a high vulnerability of an aquifer comes with a high value of hydraulic conductivity [41]. In the study, the hydraulic conductivity of the first aquifer layer was obtained based on the pumping tests applied to the available 47 monitoring wells in the Pingtung plain groundwater basin [28]. The values of hydraulic conductivity varied from 1 to 507.17 m/day (Figure 4). In the DRASTIC models, the values of hydraulic conductivity were characterized based on six different classes following the guideline proposed by Aller et al. [1].

2.7.8. Land Use (LU)

Land use is one of the essential parameters for the assessment of GVs. The contaminants triggered by anthropogenic and agricultural activities can directly or indirectly influence the groundwater quality. Similar to most areas in river deposit plains, the land use in the Pingtung plain includes agricultural fields, water bodies, building and industry, mining fields, public areas, and other areas according to the specified land-use categories. Note that the land-use parameter is valid for the modified DRASTIC and modified AHP/DRASTIC models. The agricultural areas might lead to a high potential of groundwater contamination. The highest rating of 10 was specified for the agricultural areas. Table 3 details the strategy for the scoring criteria. The study further employed the land-use maps in 1995 and 2017 to address the influence of land use on the GV variations. Figure 5 shows the details of the land-use maps for 1995 and 2017. Over 10 years of land development in the Pingtung plain resulted in a slight change of land use along the Kaoping river.
In the study, we used OK to obtain the spatial distribution of D and other factor maps for different DRASTIC models. Table 3 lists the weightings and ratings for the selected parameters in the different DRASTIC models. The weightings and ratings were modified based on the required procedures involved in the different DRASTIC models. Figure 4 shows the distribution of all factors in the DRASTIC models for the Pingtung plain groundwater basin. Note that the distributions of the parameters use the assigned rating values based on the observation data.

3. Results and Discussion

3.1. Assessment of Groundwater Vulnerability with Different DRASTIC Models

The three different vulnerability maps were generated by overlaying the weightings and ratings from eight parameter maps. The range of vulnerability indices was classified into three classes (low, moderate, and high) based on the classification method of natural breaks (Jenks) that describe the relative probability of contamination of the groundwater resources [42]. Figure 6a–c show the vulnerability indices obtained from the conventional DRASTIC, modified DRASTIC, and modified AHP/DRASTIC models, respectively. The index values were between 100 and 188 for the conventional DRASTIC model (Figure 6a). Figure 6b shows a similar GV pattern but more detailed variations of GV obtained from the modified DRASTIC model. The result in Figure 6b shows that the vulnerability indices varied from 115 to 238. In Figure 6c, the GV result based on the AHP modification shows a significant change of highly vulnerable areas when comparing the GV result in Figure 6c with that in Figure 6a,b. The main difference between results in Figure 6a–c is the land-use factor used in the models. The differences of the GV maps in Figure 6 indicate that the factor LU might play an essential role in the local variations of GV maps in the Pingtung plain. Moreover, the high-vulnerability areas in the GV maps highly correlated with the highly permeable vadose zone, aquifer media, and hydraulic conductivity. Most areas in the Pingtung plain groundwater basin exhibit low slope variations, indicating the high GV for the entire basin. The result also implies that the influence of the topography on the calculations of GVs may be small because the ratings are similar. The result shows that the modified vulnerability indices varied from 0.1374 to 0.4024 (Figure 6c). Figure 6c also reports that the highly vulnerable areas were in the townships of Gaoshu, Yanpu, and Changzhi, located in the northeastern sector of the study area. In the Pingtung plain groundwater basin, the areas of high, moderate, and low vulnerability were 34.9%, 41.6%, and 23.5%, respectively.

3.2. Validation of GVs and the Influence of Long- and Short-Term Data on the GVs

The validation of the vulnerability maps relied on comparing field measurements of common groundwater pollutants with the vulnerability index created by the DRASTIC models [43]. In general, the correlation coefficient between the vulnerability indices and actual contaminant measurements at specified locations was used to validate the accuracy of GV mapping [10]. The nitrate concentration is a non-point source contaminant that is used to evaluate the GV and contamination risk [44]. This study considered NO3–N obtained from 2010 to 2017 as the target contaminant for the validation of different DRASTIC models. Field observations of groundwater quality showed that the nitrate concentration over the period of 2010–2017 in terms of minimum, maximum, and average values was 0.01 mg/L, 14.50 mg/L, and 3.93 mg/L, respectively. We then defined the averaged data in the period from 2010 to 2017 as the long-term case and the averaged data in the period from 2016 to 2017 as the short-term case. The comparison also employed the land-use maps obtained from 1995 and 2017 and assessed the influence of land-use evolution on the GV variations.
Figure 7 shows the comparison of nitrate concentration and the normalized values of GVs obtained from different DRASTIC models. In Figure 7, the data for the comparison were the long-term case, while the land-use map was the recent observation in 2017. The result of the comparison shows the promising GV map obtained from the modified AHP/DRASTIC model. With a simple linear regression applied to the scatter plot in Figure 7, the Pearson correlation coefficient (r) and coefficient of determination (R2) were also employed to account for the accuracy of the GV maps [7].
Table 4 lists the validation of GVs for different DRASTIC models. The results reveal that recent records of nitrate concentrations in groundwater were highly correlated with the land use in 2017 as compared with the GV that considered the land use in 1995. The short-term averaged (2016–2017) DRASTIC parameters and the land-use map in 2017 had the best results in different DRASTIC models. As expected, the modified AHP/DRASTIC model performed the best in terms of the GV map by evaluating the r and R2 coefficients. The variations of r and R2 for different DRASTIC models show that the sources of DRASTIC parameters and the NO3–N concentration controlled the accuracy of GV maps. With the systematical comparison, the contamination did not remain constant through the years from 2010 to 2017, since there was a higher correlation for the short-term analyses than for the long-term ones. The results indicate the intrinsic characteristics of nitrate compounds as being soluble and mobile in groundwater. Additionally, the distribution of nitrate contamination represents the evolution of human activities that are responsible for the production of nitrate sources in time and space. The comparison also implies that the GV maps need to be modified by considering yearly changes or a shorter period of observations for DRASTIC parameters.
Table 5 lists the correlation between observed nitrate concentration and the annual GV maps. The results in Table 5 also show that the modified AHP/DRASTIC model performed the best estimations of GV maps. With the same sources of data for calculating GVs, the correlation in different years for different DRASTIC models did not vary significantly.

3.3. The Variations of GV Induced by the Change of Land Use

According to Huan et al. [45], land-use types that are related to human activities represent the major factor that directly affects the migration of nitrate in the vadose zone, which should be considered as a critical parameter in groundwater vulnerability assessment. In the study, the land-use maps in 1995 and 2017 were included in the DRASTIC models to quantify the variation of vulnerability in the period from 1995 to 2017. We assumed that the variations of physical properties such as aquifer media, soil media, topography, vadose zone, and hydraulic conductivity were negligible because of the time scales for the changes in physical properties. The parameter of depth to water used the available data of groundwater monitoring from 2010 to 2017. Calculations of the net recharge employed two land-use maps in 1995 and 2017. Furthermore, we adjusted the GV in 1995 and 2017 by using a similar classification method in 2017 for comparison purposes. Such adjustment yielded slightly lower GVs for 2017 as compared with those for 1995. Figure 8 shows the results of different DRASTIC models using different land-use maps. In general, the GV maps obtained from different DRASTIC models showed similar distributions of low, moderate, and high vulnerability to contamination. However, the area extension for low, moderate, and high vulnerability was different; for example, the areas of low, moderate, and high vulnerability were 12.7%, 43.9%, and 43.4% in 1995, while they were 14.7%, 44.6%, and 40.7% in 2017. In the Pingtung plain area, the low-vulnerability area was increased from 1995 to 2017. The high-vulnerability area was decreased, indicating that the land use and human activities in the Pingtung plain improved the reduction of the high-vulnerability area. Over the past two decades, the agriculture fields decreased by 5.3% in the Pingtung area, while the areas of land-use types describing built-up areas, traffic, and mining increased by 2.7%, 2.1%, and 0.8%, respectively. Additionally, the groundwater resources in Pingtung plain are currently over-exploited due to non-uniform precipitation and lack of withdrawal management [46]. The over-pumping leads to a significant decline in hydraulic head. As the distance from the land surface to the groundwater table increased, the potential of contaminant transport decreased. Therefore, The depth to water might not be the main contributing factor to the increase in GV. The land use should be considered as a dominant parameter in the study area.

3.4. The Mapping of Groundwater Contamination Risk

The groundwater contamination risk is the probability that the groundwater in the aquifer will become contaminated to an unacceptable level. The risk represents the combination of the occurrence probability and consequences of specified hazardous events (Equation (7)). In the study, we followed the concept proposed by Morris and Foster [33] and considered the standard of irrigation water as a 10 mg/L concentration threshold. It is interesting to understand how precipitation and groundwater usage influence the variation of contamination risk. The GV maps for the years from 2010 to 2017 were then created based on the corresponding observations of nitrate concentration in different years. Note that the GV maps used the land-use map of 2017. The DRASTIC parameters of depth to water and net recharge considered yearly averaged groundwater levels and precipitation. Other DRASTIC parameters were fixed because of negligible change over time.
Figure 9 shows GV obtained from the modified AHP/DRASTIC model, in terms of nitrate concentration, probability of occurrence, and groundwater contamination risk maps for the years from 2010 to 2017. The legend of these GV maps in each year was normalized in the range of 0 to 100 for comparison purposes. The results of GV maps show that the GV patterns were consistent over the years. The variations of high, moderate, and low vulnerability were not significant. The result implies that the parameters of depth to water and net recharge might not affect the DRASTIC models to a large extent.
The nitrate concentration maps in Figure 9 show that the eastern areas of the Pingtung plain exhibited similar distributions but slightly different intensity. However, the western areas of the Pingtung plain showed significant variations of nitrate concentration in terms of both distribution and intensity. Based on the indicator kriging algorithm applied to the observations, the probability maps of nitrate occurrence in the Pingtung plain were obtained, and the probability was normalized for comparison purposes (i.e., the third column in Figure 9). In the study, the calculations of probability maps relied on the concentration threshold (i.e., 10 mg/L) defined in the study. Figure 9 shows that there was no general pattern obtained in the probability maps because the significant agricultural activities probably changed on a yearly basis, depending on the market of the products. Comparing the GV and contamination maps gave a good match for these yearly averaged maps. In general, Gaoshu, Yanpu, Changzhi, Daliao, and Xinpi counties had high-risk irrigation water in these regions based on a combination of the vulnerability and the intensity of nitration concentration.
Based on the Equation (7), we integrated the production of the probability map, pollution map, and vulnerability map, and obtained the risk map for different years (see the fourth column in Figure 9). Note that the values of normalized contamination risk were based on the ranges between maximum and minimum values of the GW contamination risk. The results clearly show that the areas in counties of Neipu, south Changzhi, north Daliao, south Wanluan, north Chaozhou and Xinpi, and south Fangliao exhibited a high risk of nitrate pollution, and the water quality in these regions could affect the health of their inhabitants. The series of risk maps enabled the monitoring of time-varied agriculture activities in the Pintung plain. In 2012, the high risk occurred in Yanpu county, south Gaoshu county, north Changzhi county, and small areas in other counties. However, the probability of nitrate in this year was relatively low as compared to the probability maps in other years. Risk maps for years 2013 to 2017 showed continuous high risks of nitrate pollution in areas of west Yanpu, south Changzhi, north Neipu, north Wandan, Xinpi, and north Fangliao. The high-risk areas excluded Fangliao county after 2016. Moreover, Jiuru and Ligang counties recorded no risk of groundwater contamination in the past; however, in recent years, these counties showed a high risk of groundwater contamination, indicating that agricultural activities increased considerably in recent years.

3.5. Sensitivity Analysis

In this study, the sensitivity analysis used the modified AHP/DRASTIC model to assess the influence of different DRASTIC parameters on the GV results. In the study, before the sensitivity analyses, fundamental statistics were conducted to obtain general insight into the relationship between DRASTIC parameters. Table 6 summarizes the statistical properties of the eight DRASTIC parameters for the modified AHP/DRASTIC model. The mean values of topography and soil media were 0.4066 and 0.3893, indicating the strong influence of vulnerability indices. However, depth to water showed the lowest mean value of 0.1313. Hydraulic conductivity and aquifer media had mean values of 0.2223 and 0.2488, which suggest a low vulnerability of aquifer contamination. In the analyses, the net recharge, vadose zone, and land use presented moderate influences on the estimations of vulnerability. Table 6 also shows that aquifer media had the highest coefficient of variance (CV) with 59.7%, followed by hydraulic conductivity with a CV value of 52.3%. Recharge had the lowest CV value of 6.2%. The depth of water (44.7%), vadose zone (44.6%), land use (41.4%), and soil media (40.7%) represented moderate variations, while topography had the lowest variation in the study. A higher CV of a DRASTIC parameter indicates a more significant contribution of the factor to the variation in GV. On the contrary, a low CV value of a DRASTIC parameter suggests less of a contribution to GV variation [47,48].
Table 7 presents Spearman’s rank-order correlations for the eight DRASTIC parameters in the modified AHP/DRASTIC model [2]. In the study, we conducted 28 rank-order correlation analyses. However, only a few significant correlations at 95% confidence level are listed for evaluation. It can be seen that the aquifer media had an essential effect on the depth of water (r = 0.70). Furthermore, net recharge and hydraulic conductivity also contributed to the fluctuation of water depth. Their values of correlation coefficient (r) were 0.64 and 0.51, respectively. Site-specific conditions in the Pingtung plain area demonstrate that the soil media and hydraulic conductivity might have virtually no relationship (r = 0.12). The correlation coefficients between net recharge and aquifer media and hydraulic conductivity show that the rates of net recharge in various areas might be different because of the moderate influences of net recharge on aquifer material and hydraulic conductivity (r = 0.44 and r = 0.42). Additionally, the aquifer media and hydraulic conductivity slightly affect each other (r = 0.31).
Table 8 lists the statistics of the removal sensitivity analysis. The vadose zone parameter showed an average single sensitivity value of 1.8%, indicating it as the most significant parameter influencing the variation of GV. This result was consistent with the conventional DRASTIC model that showed the high theoretical weighting for estimating GVs. The removal of the topography and land-use parameters obtained sensitivity values of 1.0% and 0.9% on average, respectively. These two parameters are expected to have high sensitivity for estimating GVs. Although the theoretical weight of topography only contributed to the weight of 1, the mean rating score was 0.4066. The overall contribution to the GV estimations was relatively significant. Removing the factor of aquifer media (average sensitivity value of 0.7%) also demonstrated the importance of the parameter on GV estimations. Additionally, the removal of depth to water, soil media, and hydraulic conductivity showed that they were all relatively sensitive to the variations of the covariance index (values of 0.6%). Net recharge was one of the essential parameters of vulnerability variation with a value of 0.4%. Based on the removal sensitivity analysis, the order of variation index (I > T > LU > A > D = S = C > R) was different from the order of the theoretical weights in the DRASTIC model (i.e., D = I = LU > R > A = C > S > T). Accordingly, the removal sensitivity of DRASTIC parameters depends on the cumulative influence of the extensive data range, the mean of the rating score, and the theoretical weight of each parameter [2,20,47], as well as the covariance [48]. The finding suggests that the GV varies with the removal of parameters, and it is indispensable to apply all eight parameters to evaluate GV in the study area.
Single-parameter sensitivity analysis was used to assess the impact of the effective weight of the eight parameters on the vulnerability index compared to the theoretical weight of parameters from the conventional DRASTIC model. The term “effective weight” is a function of the value of the single parameter concerning the other parameters, as well as the weight assigned to it by the DRASTIC model [2]. Table 9 reveals that the impact of the vadose zone dominated the modified AHP/DRASTIC vulnerability index because of its mean effective weight of 24.17% against the theoretical weight of 17.86%. The calculated effective weights for net recharge (15.54%), aquifer media (10.09%), and hydraulic conductivity (9.8%) corresponded to their theoretical weights (14.29%, 10.71%, and 10.71% respectively). Soil type displayed a higher effective weight value (9.55%) than its theoretical one (7.14%), while land use showed a lower effective weight (15.31%) compared to its theoretical one (17.86%). The lowest effective weight among the parameters was exhibited by the topography, which presented a slight change of 5.42% upon comparing it to its theoretical weight of 3.57% in the DRASTIC model. Additionally, the depth to water parameter showed a high difference between theoretical and effective weight. The result indicates that this parameter might have a low contribution to GV mapping because of the distribution of groundwater levels in the study area. A large depth from the land surface to the groundwater level can restrict pollutants from reaching the groundwater body in the study area. In summary, the single-parameter sensitivity analysis suggested that the importance of net recharge, vadose zone, and land-use needs to be highlighted when assessing vulnerability using the DRASTIC model. The distributions of accurate, detailed, and representative data are essential for enhancing the overall performance of the modified AHP/DRASTIC model.

4. Conclusions

With long-term observations of groundwater quality data in Taiwan, this study assessed and quantified the evolution of GV and zones of groundwater pollution risk in the Pingtung plain groundwater basin in southern Taiwan. The results suggested that the modified AHP/DRASTIC model performed the best in terms of estimating GV, which was highly consistent with the site-specific nitrate observations. The land use is one of the crucial factors that directly reflect the effect of anthropogenic activities on GVs. The newest parameter LU can enhance the reliability of GV mapping because the LU reflects the state of the land use in the study area. The estimations of the spatial–temporal variations of GV and groundwater pollution risk enabled the identification of variable agriculture activities in the study period. In recent years, the increase of specific agriculture zones (or agriculture parks) and other anthropogenic activities yielded an increase in vulnerability to groundwater pollution. Because of the water demands, all these specific agriculture zones (or agriculture parks) are mainly developed by rivers.
Over two decades (1995–2017), the groundwater vulnerability in Pingtung plain slightly decreased due to the change in land-use types such as built-up areas, traffic, and mining areas. The changes in land-use patterns, which are highly related to the land development policy in Pingtung plain, can lead to variations in GV and, therefore, influence the strategies of groundwater resource protection and management. The groundwater contamination risk showed that Jiuru and Ligang counties recorded no risk of groundwater contamination in the past; however, in recent years, these counties exhibited a high risk of groundwater contamination, indicating that agricultural activities increased considerably in recent years. The counties of Yanpu, Changzhi, Gaoshu, and several small areas in Pingtung plain recorded high pollution risk from 2010 up to now. With farmland conditions and changes in plant type, these high-risk areas excluded Fangliao county after 2016. These variations in contamination risk are essential for decision-making and policies of regional development plans.
The concept of spatial–temporal variations in GV and pollution risk shows the importance of human activities on the groundwater environment. In this study, a land-use map was created in 2017, which might not be suitable for estimating GVs from 2010 to 2016 and in the years after 2017. Taking advantage of satellite and imaging technologies, detailed and high-frequency land-use maps can be obtained for any specific area. It is possible to monitor the dynamic variations of GV or pollution risk to track the impact of real-time human activities on the groundwater environment and the influence of health for groundwater usage. Additionally, the prediction of GV and pollution risk can also be possible if the GV estimations are coupled with real-time groundwater monitoring networks and physical groundwater models. Lastly, the physical models can predict and quantify the DRASTIC parameters that correspond to the development plans or strategies for water resources.

Author Contributions

T.-D.V. and C.-F.N. conceptualized the main idea of the paper. T.-D.V. and C.-F.N. developed and tested the models. T.-D.V. and W.-C.L. prepared the figures. T.-D.V., C.-F.N., W.-C.L., and M.-H.T. wrote the paper. All authors read and approved the final draft.

Funding

This study is jointly supported by the Ministry of Science and Technology, the Republic of China under grants MOST 107-2625-M-008 -001, MOST 107-2116-M-008 -003 -MY2, and MOST 108-2625-M-008 -007, by the Soil and Groundwater Pollution Remediation Fund in 2018 and 2019, and by the Water Resources Planning Institute under grant 107705.

Acknowledgments

We thank Water Resource Bureau and Central Geological Survey, Ministry of Economic Affairs, for providing groundwater and geological data. We also thank the anonymous reviewers for their constructive criticism and helpful comments.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Aller, L.; Bennet, T.; Leher, J.H.; Petty, R.J.; Hackett, G. DRASTIC: A Standardized System for Evaluating Ground Water Pollution Potential Using Hydrogeologic Settings; U.S. Environmental Protection Agency: Chicago, IL, USA, 1987.
  2. Rahman, A. A GIS based DRASTIC model for assessing groundwater vulnerability in shallow aquifer in Aligarh, India. Appl. Geogr. 2008, 28, 32–53. [Google Scholar] [CrossRef]
  3. Kaliraj, S.; Chandrasekar, N.; Peter, T.S.; Selvakumar, S.; Magesh, N.S. Mapping of coastal aquifer vulnerable zone in the south west coast of Kanyakumari, South India, using GIS-based DRASTIC model. Environ. Monit. Assess. 2014, 187, 4073. [Google Scholar] [CrossRef] [PubMed]
  4. Nobre, R.C.M.; Rotunno Filho, O.C.; Mansur, W.J.; Nobre, M.M.M.; Cosenza, C.A.N. Groundwater vulnerability and risk mapping using GIS, modeling and a fuzzy logic tool. J. Contam. Hydrol. 2007, 94, 277–292. [Google Scholar] [CrossRef] [PubMed]
  5. Jang, W.S.; Engel, B.; Harbor, J.; Theller, L. Aquifer vulnerability assessment for sustainable groundwater management using DRASTIC. Water 2017, 9, 792. [Google Scholar] [CrossRef]
  6. Wang, J.; He, J.; Chen, H. Assessment of groundwater contamination risk using hazard quantification, a modified DRASTIC model and groundwater value, Beijing Plain, China. Sci. Total Environ. 2012, 432, 216–226. [Google Scholar] [CrossRef]
  7. Sener, E.; Davraz, A. Assessment of groundwater vulnerability based on a modified DRASTIC model, GIS and an analytic hierarchy process (AHP) method: The case of Egirdir Lake basin (Isparta, Turkey) Evaluation de la vulnérabilité de l’eau souterraine basée sur un modèle DRASTIC. Hydrogeol. J. 2012, 21, 701–714. [Google Scholar] [CrossRef]
  8. Saida, S.; Tarik, H.; Abdellah, A.; Farid, H.; Hakim, B. Assessment of Groundwater Vulnerability to Nitrate Based on the Optimised DRASTIC Models in the GIS Environment (Case of Sidi Rached Basin, Algeria). Geosciences 2017, 7, 20. [Google Scholar] [CrossRef]
  9. Javadi, S.; Kavehkar, N.; Mohammadi, K.; Khodadadi, A.; Kahawita, R. Calibrating DRASTIC using field measurements, sensitivity analysis and statistical methods to assess groundwater vulnerability. Water Int. 2011, 36, 719–732. [Google Scholar] [CrossRef]
  10. Rupert, M.G. Calibration of the DRASTIC Ground Water Mapping Method. Ground Water 2001, 39, 625–630. [Google Scholar] [CrossRef]
  11. Liang, C.P.; Jang, C.S.; Liang, C.W.; Chen, J.S. Groundwater vulnerability assessment of the Pingtung plain in Southern Taiwan. Int. J. Environ. Res. Public Health 2016, 13, 1167. [Google Scholar] [CrossRef]
  12. Jang, C.S.; Chen, S.K. Integrating indicator-based geostatistical estimation and aquifer vulnerability of nitrate-N for establishing groundwater protection zones. J. Hydrol. 2015, 523, 441–451. [Google Scholar] [CrossRef]
  13. Allah, A.; Sedghi, Z.; Khatibi, R.; Gharekhani, M. Mapping vulnerability of multiple aquifers using multiple models and fuzzy logic to objectively derive model structures Science of the Total Environment Mapping vulnerability of multiple aquifers using multiple models and fuzzy logic to objectively derive model structures. Sci. Total Environ. 2017, 593–594, 75–90. [Google Scholar]
  14. Nadiri, A.A.; Sedghi, Z.; Khatibi, R.; Sadeghfam, S. Mapping specific vulnerability of multiple confined and unconfined aquifers by using artificial intelligence to learn from multiple DRASTIC frameworks. J. Environ. Manag. 2018, 227, 415–428. [Google Scholar] [CrossRef] [PubMed]
  15. Neshat, A.; Pradhan, B.; Dadras, M. Groundwater vulnerability assessment using an improved DRASTIC method in GIS. Resour. Conserv. Recycl. 2014, 86, 74–86. [Google Scholar] [CrossRef]
  16. Thirumalaivasan, D.; Karmegam, M.; Venugopal, K. AHP-DRASTIC: Software for specific aquifer vulnerability assessment using DRASTIC model and GIS. Environ. Model. Softw. 2003, 18, 645–656. [Google Scholar] [CrossRef]
  17. Lima, M.L.; Zelaya, K.; Massone, H. Groundwater vulnerability assessment combining the drastic and Dyna-CLUE model in the Argentine Pampas. Environ. Manag. 2011, 47, 828–839. [Google Scholar] [CrossRef]
  18. Secunda, S.; Collin, M.L.; Melloul, A.J. Groundwater vulnerability assessment using a composite model combining DRASTIC with extensive agricultural land use in Israel’s Sharon region. J. Environ. Manag. 1998, 54, 39–57. [Google Scholar] [CrossRef]
  19. Panagopoulos, G.P.; Antonakos, A.K.; Lambrakis, N.J. Optimization of the DRASTIC method for groundwater vulnerability assessment via the use of simple statistical methods and GIS. Hydrogeol. J. 2006, 14, 894–911. [Google Scholar] [CrossRef]
  20. Sadat-Noori, M.; Ebrahimi, K. Groundwater vulnerability assessment in agricultural areas using a modified DRASTIC model. Environ. Monit. Assess. 2016, 188, 19. [Google Scholar] [CrossRef] [PubMed]
  21. Toews, M.W.; Allen, D.M. Evaluating different GCMs for predicting spatial recharge in an irrigated arid region. J. Hydrol. 2009, 374, 265–281. [Google Scholar] [CrossRef]
  22. Ribeiro, L. Um Novo índice de Vulnerabilidade Específico de Aquíferos. Formulação e Aplicações. (SI: A New Index of Aquifer Susceptibility to Agricultural Pollution); Internal report; ERSHA/CVRM, Instituto Superior Tecnico: Lisbon, Portugal, 2000; p. 12. [Google Scholar]
  23. Li, R.; Merchant, J.W. Modeling vulnerability of groundwater to pollution under future scenarios of climate change and biofuels-related land use change: A case study in North Dakota, USA. Sci. Total Environ. 2013, 447, 32–45. [Google Scholar] [CrossRef] [PubMed]
  24. Finizio, A.; Villa, S. Environmental risk assessment for pesticides A tool for decision making. Environ. Impact Assess. Rev. 2002, 22, 235–248. [Google Scholar] [CrossRef]
  25. Kazakis, N.; Voudouris, K.S. Groundwater vulnerability and pollution risk assessment of porous aquifers to nitrate: Modifying the DRASTIC method using quantitative parameters. J. Hydrol. 2015, 525, 13–25. [Google Scholar] [CrossRef]
  26. Goudarzi, S.; Jozi, S.A.; Monavari, S.M.; Karbasi, A.; Hasani, A.H. Assessment of groundwater vulnerability to nitrate pollution caused by agricultural practices. Water Qual. Res. 2017, 52, 64–77. [Google Scholar] [CrossRef]
  27. Zhou, Y.; Ting, C.S.; Liu, C.W. Design of Groundwater Monitoring Networks, with Case Study of the Pingtung Plain, Taiwan; Wu-Nan Books: Taipei, Taiwan, 2003. [Google Scholar]
  28. Taiwan Central Geological Survey (CGS). Hydrogelogical Survey Report of Pingtung Plain, Taiwan; Central Geological Survey, Ministry of Economic Affairs, Executive Yuan: Taipei, Taiwan, 2002.
  29. Agriculture Engineering Research Center. Survey, Analysis and Assessment of Groundwater Quality in Taiwan Areas in 2009; Water Resources Agency, Ministry of Economic Affairs, Executive Yuan: Taipei, Taiwan, 2009.
  30. Huang, P.S.; Chiu, Y.C. A simulation-optimization model for seawater intrusion management at Pingtung Coastal Area, Taiwan. Water 2018, 10, 251. [Google Scholar] [CrossRef]
  31. Saaty, T.L. The Analytical Hierarchy Process; McGraw-Hill: New York, NY, USA, 1980. [Google Scholar]
  32. Ho, W. Integrated analytic hierarchy process and its applications—A literature review. Eur. J. Oper. Res. 2008, 186, 211–228. [Google Scholar] [CrossRef]
  33. Morris, B.; Foster, S. Cryptosporidium Contamination Hazard Assessment and Risk Management for British Groundwater Sources. Water Sci. Technol. 2015, 41, 67–77. [Google Scholar] [CrossRef]
  34. Wang, J.; Hu, Y. Environmental Modelling & Software Environmental health risk detection with GeogDetector. Environ. Model. Softw. 2012, 33, 114–115. [Google Scholar]
  35. Gogu, R.C.; Hallet, V.; Dassargues, A. Comparison of aquifer vulnerability assessment techniques. Application to the Néblon river basin (Belgium). Environ. Geol. 2003, 44, 881–892. [Google Scholar] [CrossRef]
  36. Lodwick, W.A.; Monson, W.; Svoboda, L. Attribute error and sensitivity analysis of map operations in geographical informations systems: Suitability analysis. Int. J. Geogr. Inf. Syst. 1990, 4, 413–428. [Google Scholar] [CrossRef]
  37. Napolitano, P.; Fabbri, A.G. Single-parameter sensitivity analysis for aquifer vulnerability assessment using DRASTIC and SINTACS. In Proceedings of the 2nd HydroGIS Conference, Vienna, Austria, 16–19 April 1996; Volume 235, pp. 559–566. [Google Scholar]
  38. Bisson, R.; Lehr, J.H. Modern Groundwater Exploration; John Wiley & Sons: New York, NY, USA, 2004. [Google Scholar]
  39. Hong, Y.; Adler, R.F. Estimation of global SCS curve numbers using satellite remote sensing and geospatial data. Int. J. Remote Sen. 2008, 29, 471–477. [Google Scholar] [CrossRef]
  40. Anane, M.; Abidi, B.; Lachaal, F.; Limam, A.; Jellali, S. GIS-based DRASTIC, Pesticide DRASTIC and the Susceptibility Index (SI): Comparative study for evaluation of pollution potential in the Nabeul-Hammamet shallow aquifer, TunisiaDRASTIC-SIG, DRASTIC Pesticide et Indice de Sensibilité (SI): Étude comparative. Hydrogeol. J. 2013, 21, 715–731. [Google Scholar] [CrossRef]
  41. Chenini, I.; Zghibi, A.; Kouzana, L. Hydrogeological investigations and groundwater vulnerability assessment and mapping for groundwater resource protection and management: State of the art and a case study. J. Afr. Earth Sci. 2015, 109, 11–26. [Google Scholar] [CrossRef]
  42. Krishna, R.; Iqbal, J.; Gorai, A.K.; Pathak, G.; Tchounwou, P.B.; Tuluri, F. Groundwater vulnerability to pollution mapping of Ranchi district using GIS. Appl. Water Sci. 2014, 5, 345–358. [Google Scholar] [CrossRef] [Green Version]
  43. Sahoo, S.; Dhar, A.; Kar, A.; Chakraborty, D. Index-based groundwater vulnerability mapping using quantitative parameters. Environ. Earth Sci. 2016, 75, 1–13. [Google Scholar] [CrossRef]
  44. Assaf, H.; Saadeh, M. Geostatistical assessment of groundwater nitrate contamination with reflection on DRASTIC vulnerability assessment: The case of the upper litani basin, Lebanon. Water Resour. Manag. 2009, 23, 775–796. [Google Scholar] [CrossRef]
  45. Huan, H.; Wang, J.; Teng, Y. Assessment and validation of groundwater vulnerability to nitrate based on a modified DRASTIC model: A case study in Jilin City of northeast China. Sci. Total Environ. 2012, 440, 14–23. [Google Scholar] [CrossRef]
  46. Chiang, C.-Y. Hydrogeological Survey of Pingtung Plain with the Project of Groundwater Observation Network in Taiwan; Central Geological Survey. Ministry of Economic Affairs, Executive Yuan: Taipei, Taiwan, 2002.
  47. Babiker, I.; Mohamed, M.A.A.; Hiyama, T.; Kato, K. A GIS-based DRASTIC model for assessing aquifer vulnerability in Kakamigahara Heights, Gifu Prefecture, central Japan. Sci. Total Environ. 2005, 345, 127–140. [Google Scholar] [CrossRef]
  48. Ahmed, I.; Nazzal, Y.; Zaidi, F.K.; Al-Arifi, N.S.N.; Ghrefat, H.; Naeem, M. Hydrogeological vulnerability and pollution risk mapping of the Saq and overlying aquifers using the DRASTIC model and GIS techniques, NW Saudi Arabia. Environ. Earth Sci. 2015, 74, 1303–1318. [Google Scholar] [CrossRef]
Figure 1. The Pingtung plain groundwater basin and the associated river system (TWD_1997_TM_Taiwan projection system).
Figure 1. The Pingtung plain groundwater basin and the associated river system (TWD_1997_TM_Taiwan projection system).
Water 11 02492 g001
Figure 2. The hydrogeological framework for cross-section (A–A’) in Figure 1 (modified from Reference [30]).
Figure 2. The hydrogeological framework for cross-section (A–A’) in Figure 1 (modified from Reference [30]).
Water 11 02492 g002
Figure 3. The flowchart used in the study to quantify the variations of groundwater vulnerability (GV) and pollution risk.
Figure 3. The flowchart used in the study to quantify the variations of groundwater vulnerability (GV) and pollution risk.
Water 11 02492 g003
Figure 4. Spatial distribution of the rating map of different parameter maps used in the DRASTIC models.
Figure 4. Spatial distribution of the rating map of different parameter maps used in the DRASTIC models.
Water 11 02492 g004aWater 11 02492 g004b
Figure 5. Land-use maps of the Pingtung plain groundwater basin based on the data in (a) 1995 and (b) 2017.
Figure 5. Land-use maps of the Pingtung plain groundwater basin based on the data in (a) 1995 and (b) 2017.
Water 11 02492 g005
Figure 6. Normalized vulnerability maps obtained from different DRASTIC models: (a) conventional DRASTIC model, (b) modified DRASTIC model, and (c) modified analytic hierarchy process (AHP)/DRASTIC model. Note that the land-use data for (b) and (c) represent the conditions of human or agriculture activities in 2017.
Figure 6. Normalized vulnerability maps obtained from different DRASTIC models: (a) conventional DRASTIC model, (b) modified DRASTIC model, and (c) modified analytic hierarchy process (AHP)/DRASTIC model. Note that the land-use data for (b) and (c) represent the conditions of human or agriculture activities in 2017.
Water 11 02492 g006
Figure 7. The relationship between nitrate concentration and GV maps using the land-use map of 2017: (a) conventional DRASTIC, (b) modified DRASTIC, (c) modified AHP/DRASTIC. The nitrate concentration used in the analyses was the long-term averaged concentration (i.e., from 2010 to 2017).
Figure 7. The relationship between nitrate concentration and GV maps using the land-use map of 2017: (a) conventional DRASTIC, (b) modified DRASTIC, (c) modified AHP/DRASTIC. The nitrate concentration used in the analyses was the long-term averaged concentration (i.e., from 2010 to 2017).
Water 11 02492 g007aWater 11 02492 g007b
Figure 8. The comparison of GV maps obtained from different DRASTIC models and land-use maps: (a,d) conventional DRASTIC model; (b,e) modified DRASTIC model; (c,f) modified AHP/DRASTIC model. Subfigures (a,b,c) in the first row show the results based on the land-use map of 1995. Subfigures (d,e,f) employed the land-use map of 2017.
Figure 8. The comparison of GV maps obtained from different DRASTIC models and land-use maps: (a,d) conventional DRASTIC model; (b,e) modified DRASTIC model; (c,f) modified AHP/DRASTIC model. Subfigures (a,b,c) in the first row show the results based on the land-use map of 1995. Subfigures (d,e,f) employed the land-use map of 2017.
Water 11 02492 g008aWater 11 02492 g008b
Figure 9. The variation of yearly vulnerability maps, nitrate concentration, probability of occurrence, and groundwater pollution risk maps.
Figure 9. The variation of yearly vulnerability maps, nitrate concentration, probability of occurrence, and groundwater pollution risk maps.
Water 11 02492 g009aWater 11 02492 g009bWater 11 02492 g009c
Table 1. The consistency ratios (CRs) of DRASTIC parameters in the study.
Table 1. The consistency ratios (CRs) of DRASTIC parameters in the study.
CriteriaCR
Depth to water0.02
Net recharge0.001
Aquifer media0.07
Soil media0.04
Topography0.017
Impact of vadose zone0.01
Hydraulic conductivity0.019
Land use0.02
Total0.025
Table 2. Sources of data used for creating maps of different DRASTIC parameters. DEM = digital elevation model.
Table 2. Sources of data used for creating maps of different DRASTIC parameters. DEM = digital elevation model.
Data TypeSources
(1) Precipitation and evaporation (2010–2017)Taiwan Central Weather Bureau
(2) Geology map (2002)Taiwan Central Geological Survey
(3) Hydrogeology map (2002)Taiwan Central Geological Survey
(4) Soil mapCouncil of Agriculture, Taiwan
(5) Topography (DEM 20 m)National Land Surveying and Mapping Center
(6) Observation well (2010–2017)Taiwan Water Resources Agency
(7) Borehole profile (2002)Taiwan Central Geological Survey
(8) Land use mapNational Land Surveying and Mapping Center
(9) Nitrate concentration (2010–2017)Taiwan Water Resources Agency
Table 3. Rating and weighting values used in DRASTIC, modified DRASTIC, and modified AHP/DRASTIC models.
Table 3. Rating and weighting values used in DRASTIC, modified DRASTIC, and modified AHP/DRASTIC models.
ParametersClassesDRASTICModified DRASTICModified AHP/DRASTIC
RatingWeightRatingWeightRatingWeight
Depth to water (m)0–1.51051050.38340.194
1.5–3990.2439
3–9770.1579
9–15550.0954
15–22330.0573
22–30220.0374
>30110.0248
Net recharge (mm/year)0–5014140.03510.1422
51–100330.0877
101–180660.2648
180–255880.2770
>255990.3354
Aquifer mediaFine to medium sand63630.06010.0995
Medium to coarse sand770.1615
Fine to medium gravel880.2878
Medium to coarse gravel990.4905
Soil mediaSandy loam62620.58160.0597
Loam550.3090
Clay loam330.1095
Topography (%)0–2.01011010.43390.0342
2.0–6.0990.2966
6.0–12.0550.1372
12.0–18.0330.0813
>18.0110.0510
Impact of vadose zoneSilt/clay15150.03710.2096
Fine to medium sand660.0912
Medium to coarse sand770.1604
Fine to medium gravel880.2660
Medium to coarse gravel990.4454
Hydraulic conductivity (m/d)1–513130.03030.1037
5–10220.0582
10–30440.0939
30–40660.1520
40–80880.2530
>8010100.4125
Land useAgriculture fields-1050.38880.1537
Built-up area90.2573
Waterbody70.1243
Road60.1182
Mining area30.0735
Uncultivated10.0378
Table 4. Correlation between observations of nitrate concentration and groundwater vulnerability (GV) values obtained from different DRASTIC models.
Table 4. Correlation between observations of nitrate concentration and groundwater vulnerability (GV) values obtained from different DRASTIC models.
Model.Land Use 1995Land Use 2017
Pearson Correlation (r)Coefficient of Determination (R2)Pearson Correlation (r)Coefficient of Determination (R2)
LongShortLongShortLongShortLongShort
Conventional DRASTIC0.600.610.350.370.600.610.370.37
Modified DRASTIC0.650.660.430.440.750.760.570.57
Modified AHP/DRASTIC0.820.840.680.700.840.840.700.71
Table 5. Correlation between nitrate concentration and the yearly vulnerability index obtained using different DRASTIC models.
Table 5. Correlation between nitrate concentration and the yearly vulnerability index obtained using different DRASTIC models.
Model20102011201220132014201520162017
Conventional DRASTIC0.5530.4230.5930.4390.5620.5260.5850.612
Modified DRASTIC0.5740.4980.6300.5120.5840.5340.5990.623
Modified AHP/DRASTIC0.7700.7110.8110.7090.8240.7950.7640.797
Table 6. The statistical properties of the parameters in the modified AHP/DRASTIC model.
Table 6. The statistical properties of the parameters in the modified AHP/DRASTIC model.
FactorsMinimumMaximumMeanSD *CV ** (%)
D0.02480.38340.13130.0644.7
R0.26480.33540.27460.026.2
A0.06010.49050.24880.1559.7
S0.10950.58160.38930.1640.7
T0.05100.43390.40660.0717.6
I0.03710.44540.31100.1444.6
C0.03030.41250.22230.1252.3
LU0.03780.38880.30140.1241.4
* SD: standard deviation; ** CV: coefficient of variance.
Table 7. Summary of rank-order correlation analysis for DRASTIC parameters.
Table 7. Summary of rank-order correlation analysis for DRASTIC parameters.
Correlated ParametersCorrelation Coefficient, rSignificance Level, p
Water depth and aquifer media0.70p < 0.0001
Water depth and recharge0.64p < 0.0001
Water depth and hydraulic conductivity0.51p < 0.0001
Recharge and aquifer media0.44p < 0.0001
Recharge and hydraulic conductivity0.42p < 0.0001
Aquifer media and hydraulic conductivity0.31p < 0.0001
Soil media and hydraulic conductivity0.12p < 0.0001
Table 8. Statistical results of the removal sensitivity analysis.
Table 8. Statistical results of the removal sensitivity analysis.
Variation Index (%)Removed Parameters
DRASTICLU
Minimum0.00.00.00.00.00.00.00.0
Maximum3.02.82.31.51.74.72.32.6
Average0.60.40.70.61.01.80.60.9
Standard Deviation0.520.390.490.380.201.020.410.44
Median0.50.30.70.61.01.80.60.8
Table 9. Statistics of the single-parameter sensitivity analysis.
Table 9. Statistics of the single-parameter sensitivity analysis.
ParameterTheoretical WeightingTheoretical Weighting (%)Effective Weighting (%)
MinimumMaximumAverageStandard DeviationMedian
Depth to water (D)517.861.37633.32710.1285.19510.807
Net recharge (R)414.2911.05732.26415.5412.87114.898
Aquifer media (A)310.712.00528.45510.0885.65110.351
Soil media (S)27.142.33722.8509.5494.0798.324
Topography (T)13.570.58212.2855.4181.4015.339
Impact of vadose zone (I)517.863.15645.42124.1679.03925.072
Hydraulic conductivity (C)310.711.25528.3499.8024.6229.915
Land use (LU)517.861.63031.01015.3066.31216.946

Share and Cite

MDPI and ACS Style

Vu, T.-D.; Ni, C.-F.; Li, W.-C.; Truong, M.-H. Modified Index-Overlay Method to Assess Spatial–Temporal Variations of Groundwater Vulnerability and Groundwater Contamination Risk in Areas with Variable Activities of Agriculture Developments. Water 2019, 11, 2492. https://doi.org/10.3390/w11122492

AMA Style

Vu T-D, Ni C-F, Li W-C, Truong M-H. Modified Index-Overlay Method to Assess Spatial–Temporal Variations of Groundwater Vulnerability and Groundwater Contamination Risk in Areas with Variable Activities of Agriculture Developments. Water. 2019; 11(12):2492. https://doi.org/10.3390/w11122492

Chicago/Turabian Style

Vu, Tien-Duc, Chuen-Fa Ni, Wei-Ci Li, and Minh-Hoang Truong. 2019. "Modified Index-Overlay Method to Assess Spatial–Temporal Variations of Groundwater Vulnerability and Groundwater Contamination Risk in Areas with Variable Activities of Agriculture Developments" Water 11, no. 12: 2492. https://doi.org/10.3390/w11122492

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