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

: The groundwater vulnerability (GV) assessment for contamination is an e ﬀ ective 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 modiﬁed 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-speciﬁc 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 inﬂuence 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.


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]. 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.

Study Area
The Pingtung plain groundwater basin has a total area of 1210 km 2 , 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 NO 3 − -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.
Water 2019, 11, x FOR PEER REVIEW 4 of 26 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.   The Pingtung plain groundwater basin and the associated river system (TWD_1997_TM_Taiwan projection system). 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.    Figure 1 (modified from Reference [30]).

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).
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: 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.

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 ).
where, as λ max approaches n, the judgments in the pair-wise comparison matrix become more consistent. Therefore, the difference, λ max -n, 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., b ij × b jk = b ik . The formula can be described by 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).
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.

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]. 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). 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]. 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).

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).

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).
where S, V i , and V xi 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 V i , and n represents the number of maps used for calculating V xi . 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: where Wx i is the effective weight of each parameter, and Pr i and Pw i denote the rating and weight of the parameter (x) in the i-th sub-area. Notation V i is the overall vulnerability index.

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. Table 2. Sources of data used for creating maps of different DRASTIC parameters. DEM = digital elevation model. 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.

Data Type Sources
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.

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).

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).

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).

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].
where W, P, ET, and R sf denote the net recharge (mm), precipitation (mm), evapotranspiration (mm), and surface runoff (mm). S, I a , 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 (R sf ). 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).

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. 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).

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).

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).

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].

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. 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.  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.

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. Figure 6c with that in Figure 6a,b. The main difference between results in Figures 6a-c is the landuse 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.

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

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 NO  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 (R 2 ) were also employed to account for the accuracy of the GV maps [7]. 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.    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 R 2 coefficients. The variations of r and R 2 for different DRASTIC models show that the sources of DRASTIC parameters and the NO 3 − -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.

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.
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.

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

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. 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.  Figure 9. The variation of yearly vulnerability maps, nitrate concentration, probability of occurrence, and groundwater pollution risk maps.

Sensitivity Analysis
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.

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.

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.