Stochastic-based approach to quantify the uncertainty of groundwater vulnerability

The study proposes a stochastic approach to quantify the uncertainty of groundwater vulnerability (GV) produced by classical index-overlay methods. In the analysis, the physical-based MODFLOW model has been integrated with the DRASTIC method and modified by the analytical hierarchy process (AHP) technique. Specifically, the flow fields from the MODFLOW model provide the parameters of depth to water and the associated hydraulic conductivity (K) for the DRASTIC method. The integrated loops between the MODFLOW and DRASTIC models enable the evaluations of GV maps by considering sources of uncertainty in geological parameters and stress changes in an aquifer system. In illustrating the approach for practical implementations, the study considers the uncertainty produced by the heterogeneity of K in the Pingtung Plain groundwater basin in southern Taiwan. Different degrees of K heterogeneity were assessed to quantify the impact of the K heterogeneity on the GV mappings. Results show that the integration of parameter uncertainty information from the stochastic-based approach for GV mapping can improve the accuracy and reliability of the GV assessment. The stochastic GV maps have accounted for the source of the K uncertainty. There are significant discrepancies in GV values in the spatial distribution and intensity in all GV classes. The results clarify the potential risk of groundwater contamination in the Pingtung Plain groundwater basin.


Introduction
The index-overlay methods, such as the DRASTIC and the improved DRASTIC methods, are widely used for quantifying GV because of their simple concept, minimal data demand, and well definition of GV. The GV maps are the basis for showing the conditions of the groundwater contamination potential for an aquifer system and enable the evaluation of proposed management solutions to preserve groundwater quality (Council 1993). Offering efficient management strategies typically relies on accurately characterizing a specific aquifer system. However, the characterization procedure is challenging because of the uncertainty induced by the geological features, climate conditions, and anthropogenic activities that are considered a major effect on the water cycle and sea level rise (Council 1993;Koutsoyiannis 2020).
The GV uncertainty consists of model-related uncertainty resulting from system processes' inadequate or incomplete representation. Also, data-related uncertainty is another source of uncertainty resulting from input data, even if the used model is correct (Loague & Green 1991). Thus, unpredictability is likely even when the models or data are perfect and highly-accurate input observations, the uncertainty of a system can still be high due to the existence of chaos (i.e., intrinsic sensitivity to the initial conditions). It has been shown that in highly-complex dynamic systems, the unpredictability (or predictability) should always be accompanied by a time-window length since it seems that somewhat both randomness and predictability coexist and are intrinsic to natural systems, which can be deterministic and random at the same time, depending on the prediction horizon and the time scale (Dimitriadis et al. 2016). For example, even a chaotic system seems never to be wholly unpredictable but somewhat predictable for a specific timehorizon or else time-window (which can be large or small), were passed that window, the predictability becomes so uncertain that we may support there is no point in predicting this system pass that time-horizon (Lorenz 1963). Imperfect models and data are the norms rather than the exception. Few published GV assessments accounted for uncertainty from either model or data errors (Coppola et al. 2013;Gharekhani et al. 2022;Gurdak et al. 2007;Loague et al. 2012;Sadeghfam et al. 2021). These approaches usually implied the apparent certainty in vulnerability assessments rather than the underlying uncertainty from different sources. Little attention has been paid to the uncertainty induced by the aquifer heterogeneity and the propagation of such uncertainty on the generated thematic maps for GV evaluations (Burrough 1986;Goodchild & Dubuc 1987;Moges et al. 2021;Ni et al. 2010).
Uncertainty analyses are typically used to evaluate the spatiotemporal variability and the propagation of uncertainty in model calculations (i.e., variance in model outputs and the sufficiency of existing spatial observations) (Dimitriadis & Koutsoyiannis 2018;Eisenberg et al. 1989;Moges et al. 2021). Such techniques have been used extensively to develop site selection criteria and design radioactive waste-disposal facilities (Brandtetter & Buxton 1989) but not for vulnerability assessments. Methods for evaluating uncertainty associated with measurement errors can be grouped into the following five categories (quoted from Brandtetter and Buxton 1989): (1) Classical statistical variance component analysis, (2) First-order uncertainty analysis (FOUA), (3) Statistical sampling methods (e.g., Monte Carlo simulation, Latin hypercube sampling, discrete-event simulation, and bootstrapping methods), (4) Stochastic modeling approaches, and (5) Bayesian methods. However, only FOUAs statistical sampling methods and stochastic modeling techniques have been applied to vulnerability assessments. The studies of Jury and Gruber (1989), and Small and Mular (1987), presented examples of stochastic approaches to evaluate uncertainty associated with climatic and soil variability in GV estimations. Loague et al. (2012) employed retardation and attenuation factors combined with the FOUA approach to characterize the data-related uncertainties on near-surface vulnerability assessments in the San Joaquin Valley. Coppola et al. (2013) used the stochastic process to assess the effect of vertical heterogeneity on the variability of solute transport in the vadose zone in mapping the GV on a regional scale. Lin et al. (2017) focused on the flow model uncertainty. They proposed an approach that used MODFLOW and stochastic simulation to quantify the uncertainty of hydraulic heads and hydraulic conductivity in multiple aquifer systems. Gharekhani et al. 2022 employed Bayesian Model Averaging (BMA) to assess the uncertainty of groundwater vulnerability in Lake Urmia. They built three Artificial Intelligence models, including ANN (Artificial Neural Network), GEP (Gene Expression Programming), and SVM (Support Vector Machines) as the inputs of BMA, whereby quantifying the uncertainty of model errors. None of the previous investigations directly integrated a physical-based numerical model with the DRAS-TIC method to quantify the propagation of the uncertainty induced by natural heterogeneity or limited observations. The physical-based numerical models have been recognized as practical techniques for assessing groundwater management and protection. The models consider groundwater response and hydrogeological characteristics to predict the future state and available groundwater resources based on climate conditions and human activities. The MODFLOW model is a well-developed computer code that solves general flow governing equations and is often used to assess the intrinsic susceptibility of an aquifer system (Harbaugh et al. 2000). The MODFLOW model includes several process-based modules to simulate the interactions of components in aquifer systems. Numerous studies have conducted simulations of aquifer systems based on the MODFLOW model and the associated specific packages for site-specific problems (Boufekane et al. 2022;Huang & Chiu 2018;Lin et al. 2017;Mahmoudpour et al. 2016;Mi et al. 2016;Singh 2013). Recent advances in parameter estimation techniques have also provided opportunities to analyze various combinations of calibrated parameters in the flow models and assess the effects of spatial heterogeneity and uncertainty on the groundwater flow (Hill and Tiedeman, 2007;Ni et al. 2010). The Monte Carlo simulation (MCS) is a classical stochastic approach for quantifying the uncertainty between input and output parameters. The main advantage of the MCS is the physically intuitive concept that integrates a series of simulations to evaluate the propagation of uncertainty in deterministic physical systems (Lee et al. 2019;Lu et al. 2018;Ni & Li 2006. The physical models, integrated with the stochastic approaches, could quantify the heterogeneity of hydrogeological settings with limited samples (Bianchi et al. 2015;Jang & Liu 2004;Lin et al. 2017). The hydraulic conductivity in groundwater models is a highly uncertain parameter recognized as a significant impediment to precise predictions (Kitanidis 1997;Lin et al. 2017). These input uncertainty, therefore, hydraulic head and conductivity, are the input parameters of the DRASTIC model, enhancing accuracy and reliability by adding the uncertainty information on GV assessment.
This study aims to develop an approach that uses stochastic simulation and a physical-based groundwater flow model incorporated with the index-overlay method. It analyses the uncertainty of input parameters for groundwater vulnerability (GV) assessment, whereby quantifying the stochastic-based GV map in the Pingtung groundwater basin in southern Taiwan. Generally, the groundwater flow model (MODFLOW) is set up and calibrated to simulate the groundwater system's physical responses. Critical parameters such as depth to water and hydraulic conductivity from the calibrated model are the modified-DRASTIC input coupled with the analytical hierarchy process (AHP) theory for quantifying the original GV map. Sequential Gaussian Simulation (SGS) is employed to obtain the limited sets of realization of hydraulic conductivity, simulate hydraulic heads' realization, and quantify these parameters' uncertainty. The information on hydraulic head and conductivities' uncertainties will be added to the depth of water and hydraulic conductivity of the modified-AHP-DRASTIC model. Finally, the stochastic-based GV mapping is validated by measured nitrate concentration using the Pearson correlation coefficient. The model-related uncertainties identified in the study are expected to provide the opportunity for improving GV assessment in the future. Also, the results from the research may help decision-makers plan and protect groundwater resources effectively.

Study area
The Pingtung Plain is located in the southwestern region of Taiwan, separated by the southern edge of the Alishan mountain range in the north and Lingkou hills in the west (Fig. 1). It is bordered by the Taiwan Strait in the south, separated by the Chaozhou Fault in the east and the Dawu Mountain Range at the Central Mountain Range's southern edge. The plain is about 1210 km square. The groundwater basin area is about 55 km long from north to south and 22 km wide from east to west. The terrain is gently inclined from northeast to southwest. Major rivers, such as the Gaoping River, Donggang River, and Linbian River, cross through the area and into the Taiwan Strait.
The Pingtung Plain has a subtropical climate with abundant rainfall but is quite uneven (Hsu et al. 2007). In terms of the temporal, the rain from May to October (wet season) accounts for more than 88% of the total rainfall, and the remaining six months (dry season) only account for less than 12%. Therefore, there are distinct periods of groundwater variation in this area. Precipitation partly seeps into the groundwater through surface water and becomes a groundwater recharge source. During the wet season (May to October), the flow accounts for 86% of the annual infiltration through the riverbed to the unconfined aquifer.
The plain is partitioned generally into a proximal fan, mid-fan, and distal fan (Taiwan Central Geological Survey (CGS), 2002). The proximal fan at the eastern and northern boundaries are crucial regions for natural aquifer recharge. The alluvial fan group of the late Quaternary age mainly controls the unconsolidated sediments in the study area from the east to the west (Ting et al. 1998). The hydrogeological stratification from the Taiwan Central Geological Survey investigation is approximately 250 m in depth (Taiwan Central Geological Survey (CGS), 2002). The aquifers in the study area are quite complex and divided into seven layers, illustrated in Fig. 2, aquifer 1 (F1), aquitard 1 (T1), aquifer 2 (F2), aquitard 2 (T2), aquifer 3-1 (F3-1), aquitard 3 (T3), and aquifer 3-2 (F3-2) in the vertical direction [20]. Daily groundwater level data were recorded by 126 groundwater monitoring sensors installed in 55 observation wells spreading across the four aquifers ( Fig. 1).

The flowchart of research
The study proposes an approach that integrates the stochastic method and physical-based groundwater flow model incorporated with the index-overlay method to analyze the uncertainty of input parameters for groundwater vulnerability (GV) assessment, whereby quantifying the stochasticbased GV map (Fig. 3). Firstly, the modified-DRASTIC, coupled with the analytical hierarchy process (AHP) technique, is employed to create the GV map. A set of realizations is produced using the geostatistical structure of measured natural log hydraulic conductivity and the mean of natural log hydraulic conductivity generated from the calibrated GW flow model. They are then assigned to the deterministic model (MODFLOW) to simulate hydraulic heads' realization, whereby the uncertainty of hydraulic head and hydraulic conductivity are quantified. The hydraulic head and conductivity variance will be added to the water depth and hydraulic conductivity of the modified-AHP-DRASTIC model. Finally, the stochastic-based GV mapping is validated by measured nitrate concentration using the Pearson correlation coefficient.

Modified-AHP-DRASTIC model
The modified model, coupled with the analytical hierarchy process (AHP) theory, is used to assess the study area's GV. Vu et al. (2019) described the modified-DRASTIC composes eight parameters that control the dynamic of groundwater contaminants. These are, namely, depth of water (D), net recharge (R), aquifer media (A), soil media (S), topography (T), the impact of vadose zone (I), hydraulic conductivity (C), and land use (LU). The ratings and weightings of these parameters are assigned to the ranges or zones regarding site-specific hydrogeological Fig. 1 The Pingtung plain groundwater basin and the associated river system settings. The DRASTIC index of the GV assessment sums up the ratings and weightings for each factor.
where w and r denote the weight and rating of each parameter.
In the study, Analytical Hierarchy Process (AHP) theory developed by Saaty (1980) was employed to modify the weightings and ratings of DRASTIC parameters. This research uses Vu et al. (2021) previous works for data Fig. 2 The hydrogeological framework for cross-section (A-A') and (b1-b5) in Fig. 1 Stochastic Environmental Research and Risk Assessment (2023) 37:1897-19151901 preparation and processing the modified-AHP-DRASTIC incorporated with a physical-based numerical GW flow model. Vu et al. (2021) created all parameter maps that are fundamental to this study. The sensitivity analysis was discussed in detail in the previous publication of Vu et al. (2019). The impact of the vadose zone and recharge is the most sensitive factor in the GV mapping. The latter is hydraulic conductivity. However, based on the stochasticbased approach and physical-based numerical model, we want to quantify the propagation of the uncertainty induced by the natural heterogeneity of geologic media represented by the hydraulic conductivity. Therefore, hydraulic conductivity is considered an input to the stochastic simulation.

Numerical model of groundwater flow
MODFLOW-2000 was employed to build a numerical model of groundwater flow in the Pingtung plain (Harbaugh et al. 2000). Further, the PEST code developed by Doherty (2010) was used as a primary tool to optimize the model parameters. The calibration of recharge, hydraulic conductivity, and groundwater extraction employs the pilot point approach. Observation data of 2017 from the monitoring network (Taiwan WRA, 2012) was considered to model the GW flow in the study area. Also, groundwater levels in 2015 and 2016 were supposed to verify the output of model calibration. This study uses Vu et al. (2021) previous works for the numerical GW flow model. In this part, we briefly summarize the conceptual model from previous research. The study area was discretized into 90 rows and 40 columns with a grid dimension of 1000 9 1000 m. The aquifer system was vertically divided into three model layers. The topmost layer is assumed to be an unconfined aquifer. The middle layer's geological formation has a very low transmissivity (Ting et al. 1998); layer two is supposed to be semi-pervious. The bottom layer is assumed to be an unconfined aquifer. Natural aquifer replenishment sources are rainfall, rivers, and subsurface inflow from the northern upstream part of the fluvial fan.
The long-term average precipitation for 2010-2017 was considered to calculate the recharge rate using the water balance equation (Bisson and Lehr, 2004). It has been adjusted during the calibration process. The hydraulic conductivity was estimated from the pumping test obtained from the Pingtung plain's hydrogeological investigation (Taiwan Central Geological Survey (CGS), 2002). Groundwater plays a crucial role in the study area for supplying water to agriculture and domestic uses. However, the total amount of GW extraction is still precisely unknown; each county's GW right, published by Taiwan's WRA in 2013, was employed as an input for model simulation.
The low permeability water-bearing properties of the formation surround the Pingtung groundwater basin; the northern, eastern, and western parts were assigned as noflow boundaries. Specified head boundaries were simulated at the lateral recharge point of river valleys and the model's southern border corresponding to the sea. The network of Kaoping, Donggang, and Lipien rivers were assigned as head-dependent flux boundaries.

Stochastic numerical simulations
Sequential Gaussian Simulation (SGS) as a stochastic method has been developed to simulate a variable with a potential for dealing with spatial uncertainty (Baskan et al. 2010). In response to the interpolation method's inadequate spatial uncertainty measures, stochastic simulation of the SGS can be applied to account for the magnitude of the values with reasonable confidence intervals (Jang and Liu 2004). The advantage of Sequential Gaussian Simulation (SGS) is to address the smoothing effect, which should be considered to avoid deterministic methods using various stochastic realizations (Bai and Tahmasebi 2022). In the estimation process for each unknown point, both actual values at sampled locations and the estimations at the previously simulated points are considered. Moreover, since various random paths are generated which can pass the unsampled points in different orders, SGS can produce several equally probable realizations to assess and evaluate the uncertainty (Verly 1993). Further, other approaches usually use the entire sampled points and all the previously simulated points for estimating the unknown point. Meanwhile, SGS only considers the nearest points within a local neighborhood, which can efficiently improve the computational speed and reduce the size of the covariance matrix. However, there are other alternatives, such as the non-linear method for the preservation of distribution function, autoregressive models, and symmetric-movingaverage scheme that is very useful for stochastic generation and the derivation of confidence intervals (Dimitriadis et al. 2021;Dimitriadis and Koutsoyiannis 2018;Koutsoyiannis 2016Koutsoyiannis , 2017. More advanced algorithms could enhance the efficiency of spatial uncertainty. In the study, Sequential Gaussian Simulation (SGS) was employed to generate the random field of hydraulic conductivity. The SGS simulation creates realizations of normal random variables and performs a Gaussian transformation of the available measurements. A simulated value at a visited point is randomly drawn from the conditional cumulative distribution function, defined by the Kriging mean and variance, based on neighborhood value. At a randomly visited point, the simulated value is conditional on the original and previously simulated data. Finally, the simulated normal values are transformed back into the simulated values for the original variables. The process is repeated until all points are simulated for each realization throughout the grid (Deutsch and Journel 1998). The procedure can be equally applied to unconditional simulation and conditional simulation. For unconditional simulation, the first point is generated as an independent random variable with the desired unconditional mean and variance. A second data point is generated conditioned on the first one in the manner. For conditional simulation, the generation of random fields begins with the already measured data. Subsequent data points are generated in unconditional simulation through the conditional relationship (Gómez-Hernández & Journel 1993).
A semi-variogram represents a measure of spatial variability between sample values. The experimental semivariogram calculates the half-mean squared between pairs of data points separated by a lag distance, vector h: where c(h) is the semi-variance, N(h) is the number of data pairs separated by a lag distance h, Z(x i ), and Z(x i ? h) is the hydraulic conductivity at location x and x ? h, respectively. In this study, the variogram models, which are spherical, Gaussian, and exponential, fit the semi-variogram. The least-squares method is employed to match the experimental semi-variogram that creates a best-fit model with the highest coefficient of determination (R 2 ) and the minimum sum of residual squares (RSS). The spherical, Gaussian, and exponential models are represented by Eq. (3)-Eq. (5) as follows (Isaaks & Srivastava 1990): Gaussian Exponential where c 0 denotes the nugget effect, c and a are the sill and range, and h represents the lag distance. The fitting results of the experimental semi-variogram generate c 0 , c, and a, respectively.

Simulation of groundwater flow
In the study, the steady-state GW model, which was assumed to be near quasi-steady-state conditions, was calibrated with the automatic parameter estimation of recharge, hydraulic conductivity, and pumping rate by using the observation data of averaged groundwater levels in 2017. Figure 4 shows the contoured comparison of the observed and simulated GW levels. It can be seen that the computed hydraulic heads quite fit the observed heads. The difference between the observed and simulated head was minimized at 41 wells in aquifer 1, 29 wells in aquifer 2, and 14 wells in aquifer 3 (Vu et al. 2021). The mean error (ME), absolute mean error (MAE), and root mean square error (RMSE) in the calibration process were 0.03 m, 0.48 m, and 0.61 m, respectively. Model verification is essential to show the best fit of the observation data based on the calibrated model (Anderson et al. 2015). Thus, the calibrated model in the study is verified by the observation data of 2015 and 2016. According to Vu et al. (2021), the RMSE of the model verification was 0.79 for 2016 and 0.85 for 2015, respectively. Further, the standard deviation over the mean of calibration and verification data for 2015-2017 varied from 0 to 9.9 m, with an average of 1.76 m (Fig. 5). The large error-values in calibration may be due to the existence of the Hurst phenomenon, which (by definition) increases the variability of the process, thus, making it difficult for any model to capture it with small calibration (validation, confirmation) errors. Overall, the calibration results expose that the GW flow model performs very well and could be used as a deterministic forward model for simulating the set of realizations.
In the last publication (Vu et al. 2021), the model calibration processes were described in detail. 82, 81, and 48 pilot points of hydraulic conductivity were adjusted for aquifer 1, aquifer 2, and aquifer 3 in the calibration process. Table 1 shows the values of hydraulic conductivity that were employed to calibrate the model. The min and max values of field measurements for hydraulic conductivity were lower and upper bound in the PEST package. The mean value of hydraulic conductivity from field measurements was used as the initial value of parameter estimation. The averaged values of pilot points for hydraulic conductivity after the calibration process are shown in Table 1. These values will be used in the stochastic simulation for analyzing the uncertainty of input parameters for GV assessment.
Accordingly, the head's fluctuation at the observation wells should also be considered to exhibit a rough idea of how the head's variations are in the study area. The timeseries data of monthly groundwater levels in three years (2015-2017) obtained at the monitoring wells is employed to analyze the variability of the hydraulic head.  presents the standard deviation (SD) of groundwater levels in the Pingtung plain, varying from 0 to 9.9 m, with an average of 1.76 m. It can be seen that the east of the study area has a high SD of the head because this is the main replenishment region of groundwater resources that are bounded by foothills and valleys. Further, the fluctuation performs slightly at some observed wells, while the areas far away from the monitoring points have a low variability of the head.

Analysis of semi-variogram
The semi-variogram is used as an indicator of the spatial dependence of the field (i.e., groundwater level or concentration of contaminants). Identifying which model can fit the experimental semi-variogram is challenging. There are many studies to propose a framework based on separable covariance functions that modify the separate spatialtemporal covariance models (Dimitrakopoulos and Luo 1994;Rouhani and Myers 1990). Also, the other method based on the joint spatial-temporal dependence between observations, that is, space-time Kriging (diffusion and Spartan covariance functions), was developed (Heuvelink and Egmond 2010;Jost et al. 2005). Varouchakis and Hristopulos 2019 concluded that the Spartan and nonseparable diffusive variograms (space-time Kriging) could fit the experimental variogram of the residuals when compared to the separable variograms (Matérn model) for groundwater processes.
In the study, the semi-variogram fitting is conducted by using a least-squares method with the lowest fitting errors. The trial and error method is employed to identify which model can fit the experimental semi-variogram. In the study, 47, 45, and 20 locations of pumping tests for aquifer 1, aquifer 2, and aquifer 3, respectively (Table 1), were employed to analyze the experimental semi-variogram. The results show that the spherical model was a best-fit model for three experimental semi-variogram of measured hydraulic conductivity (lnK) in aquifer 1, aquifer 2, and aquifer 3. Figure 6 shows the fitting parameters of the theoretical models. The semi-variogram of measured hydraulic conductivity presents the small nugget effect in three aquifers that are 0.017, 0.063, and 0.007 [ln(m/day)] 2 , respectively. The largest sill is 1.716 [ln(m/day)] 2 for aquifer 1, and those of aquifer 2 and aquifer 3 are 0.656 and 0.738 [ln(m/day)] 2 . The longest range of the fitted spherical model is 42290 m in aquifer 1. Ranges for aquifers 2 and 3 are 10,832 m and 18,509 m, respectively. Due to the insufficient pairs at above 50 km correlation length for aquifer 1 and aquifer 2, the data of correlation length within 50 km are used to validate the stationary assumption in the experimental semi-variogram. Also, the number of pairs for aquifer 3 has a correlation length of 30 km for validation. Table 2 shows a summary of the fitted parameters in the semi-variogram. It can be seen that there is a good match between a theoretical semi-variogram model and a fielddata experimental semi-variogram. The geostatistical structure of the field-data measurement has determined the variability of lnK through the mean, range, and sill values. The semi-variogram analysis exposes that the spatial distribution of hydraulic conductivity in aquifer 1 is the most continuous in the rest aquifers. Aquifer 2 has more  6 Best-fit theoretical and experimental semi-variogram of hydraulic conductivity in three aquifers heterogeneous hydraulic conductivity than aquifer 3 and aquifer 1. Figure 2 confirms the spatial heterogeneity and continuity of hydraulic conductivity in each aquifer equivalent to the semi-variogram analysis.

Conditional and unconditional geostatistical simulations
For analyzing the uncertainty of hydraulic conductivity caused by spatial heterogeneity, the variability of permeability is assigned as a random variable. The stochastic numerical simulation is employed to generate a random field of permeability based on the geostatistical structure, honoring the field-data measurement.
In the study, geostatistical parameters of measured natural log hydraulic conductivity shown in Table 2 and the mean of calibrated natural log hydraulic conductivity from the GW flow model presented in Table 1 are considered for generating a set of realizations by the SGS. Eggleston et al. (1996) state that it needs 50 -400 realizations to obtain a stable mean SGS realization. Thus, the stochastic simulation has conducted 300 sets of ln(K) realizations for aquifer 1, aquifer 2, and aquifer 3 in the study area. Figure 7 shows the average standard deviation (SD) of 300 ln(K) realizations. It can be seen that the average SD at 70 simulated realizations becomes gradually stable. Also, the mean SGS realization of hydraulic conductivity from 70 to 300 realizations varies insignificantly (Table 3).
Accordingly, compared with the geostatistical parameters of measured ln(K) and calibrated mean ln(K) reveals that the realizations of hydraulic conductivity simulated by the SGS may capture the geostatistical characteristics of measured ln(K) data and calibrated mean ln(K). It suggests that a random field of 70 simulated realizations selected to insert as input in the deterministic model can reproduce the spatial heterogeneity and variability of hydraulic conductivity. Note that the hydraulic head and conductivity uncertainty in aquifer 1 is only analyzed because of the GV assessment in an unconfined aquifer (first layer).
In the study, unconditional and conditional simulation of 70 realizations has been conducted to evaluate the propagation of hydraulic conductivity uncertainty to hydraulic head's variation. For the unconditional simulation, Fig. 8a shows the spatial distribution of standard deviation of 70 realizations for hydraulic conductivity in aquifer 1 with the minimum and maximum SD values of 1.49 and 4.40 m.day-1, respectively, while the averaged SD value is 2.31 m.day-1. It can be seen that the high SD of hydraulic conductivity distributed in the southern study area means that the value of hydraulic conductivity may vary significantly in the calibration process because of the calibrated mean ln(K). Figure 8b presents the standard deviation of the hydraulic head in aquifer 1 from the deterministic model MODFLOW. The SD value varies from 0 to 10.5 m, with an average value of 3.21 m. The northern area has a considerable hydraulic head variation, while almost all zones close to the observation wells have low variance. It means that the high uncertainty of the hydraulic head is generally located far away from the observation wells. Further, the propagation of hydraulic conductivity variability to the hydraulic head is insignificant, and the  (2023) 37:1897-1915 1907 variation of the hydraulic head at the location of high SD hydraulic conductivity is slight, which means that the factors dominate the uncertainty of the hydraulic head that the real conditions might cause, e.g., pumping rates or climate conditions. For conditional simulation, due to the constraints of observed points, the SD value of hydraulic conductivity is 1.19 to 2.3 m.day-1 (Fig. 8c), while the hydraulic head's SD is 0 to 3.10 m, respectively (Fig. 8d). The high variability of hydraulic conductivity only occurs at the location far away from observed points, distributed in the west. The spatial distribution of high head fluctuation is also similar to the hydraulic conductivity pattern.
Overall, comparing the head's variation in the unconditional, conditional simulation, and observation (Figs. 5 and 8), there is a large discrepancy in the distribution pattern and intensity. However, the SD value in unconditional simulation and observation is not much different.
Generally, the dense river network is a significant factor in controlling the variation of groundwater levels in the Pingtung plain. Also, the boundary conditions might control the physical responses near the boundaries of the plain. Thus, it is assumed that the effects of boundaries may be insignificant as they compete with the GW flow model's local stress. The spatial distribution and intensity of hydraulic conductivity in the unconditional and conditional simulation have a large discrepancy and occur in the hydraulic head. The DRASTIC parameters are classified into their value range for assessing the GV variation. If the head or hydraulic conductivity variation is small, there is no change in the parameter's rating. Hence, it is challenging to see the variation of GV. As discussed above, the uncertainty of hydraulic conductivity and head in the unconditional simulation performs higher than in the conditional simulation. The uncertainty value obtained from the unconditional simulation is employed to calculate the stochastic-based GV uncertainty. Figure 9 presents the spatial distribution of various parameter rating maps used in the stochastic-based GV. The rating map of depth to water and hydraulic conductivity has been calculated by integrating the variance and the value of simulated head and hydraulic conductivity from the physical-based GW model. The rest of the parameter maps were created by a different type of data set based on the previous work of Vu et al. (2021). The stochastic-based GV was mapped based on eight DRASTIC parameter maps overlaid by their weightings and ratings under the ArcGIS environment (Fig. 10a). The GV classes were categorized into low, moderate, and high using the classification method of natural breaks (Jenks). The value of the index 0.147-0.256 presents low vulnerable zones, while the index 0.257-0.317 and 0.318-0.405 exhibit moderate and high vulnerable zones, respectively. As seen in Fig. 9a, approximately 18.67% of the study area has low GV. The moderate and high GV account for 38.15 and 43.18%, respectively. The high GV is distributed in the counties of entirely Gaoshu, Yanpu, and Changzhi. The partly townships of Meinong, Xinpi, Qishan, and Daliao have also occurred high GV. The rest of the counties have a low to moderate GV.

Assessments of the stochastic-based GW vulnerability
The study added the uncertainty information of hydraulic head and conductivity to stochastic-based GV mapping. Compared with the GV map created by the last publication of Vu et al. (2021), called a physical-based GV map that is a baseline case for comparison, the GV value varies significantly (Fig. 10b). The changed value of GV has been converted to the percentage of change in GV for better visualization. Figure 10b shows the significant discrepancies of GV values in all GV classes' spatial distribution and intensity. When the uncertainty information of input parameters is added to GV mapping, it can be seen that the decline of low and moderate GV from 24.4% and 40.1% in the physical-based GV map to 18.67% and 38.15% of the study area in stochastic-based GV map. Besides, there is an increase in high GV from 35.6% to 43.18%, respectively. Also, the variation of the GV value can lead to GV class change (Fig. 10c). Figure 10c shows that about 5.83% of the GV zones have changed from high to moderate. The GV class changes from moderate to low, low to moderate, and moderate to high account for 4.64%, 0.74%, and 0.52% of the GV zones, respectively.

Validation of GV mapping
According to White et al. (2013), nitrate concentration (NO 3 -) in groundwater is the most common indicator of human consequences that are regularly used to verify the reliability of a GV assessment and to validate the rating  (Arauzo 2017;Elçi, 2017;Saida et al. 2017;Vu et al. 2019). The presence of NO 3 in groundwater reflects the contamination from human activities, such as agricultural and industrial activities or wastewater. The Pingtung plain is a crucial agricultural production area in Taiwan. Also, nitrate contamination has recently been recognized as the main factor of groundwater contamination (Agriculture Engineering Research Center 2009). The study used 71 measurement points of nitrate concentration taken from 2016-2017 to validate the GV mapping. The average nitrate value is 4.14, while the minimum and maximum nitrate concentration values are 0.01 and 16.10, respectively (Table 4). The monitoring locations of nitrate concentration are presented in Fig. 1. The method used for GV validation is the Pearson correlation coefficient (r) and coefficient of determination (R 2 ) implemented in the investigation of Assaf andSaadeh (2009) andStigter et al. (2006). The result shows that the Pearson r is 0.877 and R 2 is 0.77 (Fig. 11), which indicates the best match between GV and nitrate concentration.
In addition, compared with the GV validation from the last publications (Vu et al. 2019(Vu et al. , 2021, quantifying parameter uncertainty from the GW model can improve the accuracy and reliability of the GV map. The Pearson r and R 2 from the traditional approach of modified-AHP-DRASTIC (Vu et al. 2019), which relies on field data, are 0.84 and 0.71. Meanwhile, the proposed approach of the physical MODFLOW model and index-overlay DRASTIC method had Pearson r and R 2 of 0.85 and 0.717, respectively. With the information on hydraulic head and conductivity uncertainty, the improvement of the GV method become feasible because the physical-based model, Fig. 9 Parameter rating maps of GV mapping in the study (modified from Vu et al. 2019) coupled with the stochastic approach, can quantify the natural conditions and human activities.

Uncertainty analysis of GV mapping
The uncertainty of GV is also quantified, as shown in Fig. 12. The SD value of 70 realizations for normalized-GV varies from 3.71 to 12.02, with an average of 6.76. Note that we have converted the standard deviation of the modified-AHP-DRASTIC to the normalization SD value of GV for better presentation. It can be seen that the spatial distribution of high and moderate SD of GV is similar to GV class variation. This explains the variation of GV value in these areas is large enough to change GV classes, as seen in Fig. 10c. High GV uncertainty is mainly distributed near the rivers, significantly controlling the variability of groundwater levels and hydraulic conductivity. The   . 11 The scatter plot of the correlation between nitrate and stochastic-based GV propagation of hydraulic head uncertainty to GV can be seen clearly in Fig. 10b, in which the high variation of GV value has the same pattern as the hydraulic head uncertainty. This proves that the variability of groundwater level near the rivers affects slightly to GV map. Besides, the contribution of hydraulic conductivity uncertainty to GV is also recorded; it combines with the head's variation to change the GV class. The effects of parameter uncertainty for GV depend on their classes' variation of parameter values. If its intensity Fig. 12 The standard deviation of 70 stochastic realizations for GV changes are significant or slight, leading to GV class change or the GV value varies. Further, GV uncertainty also stems from the input data of GV components that are sparse, imprecision, and uncertainties. Hence, the magnitude of GV uncertainty varies with the magnitude of its parameters.
Overall, integrating the uncertainty information of the input parameters may affect the variation of GV slightly. GV assessment clarifies the potential risk of groundwater contamination by adding the uncertainty of model parameters. Also, there are significant discrepancies in GV values in the spatial distribution and intensity in all GV classes by considering the input uncertainty in the GV mapping.

Conclusion
This study has combined the physical MODFLOW model, stochastic approach, and DRASTIC method for improving the GV prediction in the Pingtung plain groundwater basin. The research focuses on assessing the spatial continuity and heterogeneity of hydraulic conductivity, quantifying the uncertainty of hydraulic head and conductivity. The SGS method is selected to simulate a set of realizations. A random field of 70 realizations can capture the geostatistical characteristics of measured ln(K) data and calibrated mean ln(K). The standard deviation of hydraulic conductivity has a high value in the southern study area-the uncertainty of hydraulic head distributed in the northern area located far away from the observation wells.
GV mapping results show the large discrepancies of GV values in both the spatial distribution and intensity in all GV classes when the uncertainty information of input parameters was added to GV mapping. The uncertain information of the input parameters may affect the variation of GV slightly. Compared with the GV map validation from the traditional approach and physical-based GV, the correlation between the stochastic-based GV and nitrate concentration shows the best in terms of GV validation. It can be concluded that the quantification of parameter uncertainty from the GW model can improve the accuracy and reliability of the GV map.
The GV method's improvement based on integrating a physical GW model, stochastic simulation and indexoverlay method is not limited to other physical models or various index-overlay methods. In future research, the Hurst phenomenon should be studied intensively, which governs the key hydrological-cycle process. Since the groundwater processes are related to the water-cycle (precipitation, evaporation, temperature, humidity, and so on), this phenomenon is expected to be also apparent. Further, this behavior significantly affects the uncertainty analysis for different DRASTIC parameters that are the dominant factors of GV assessment. Moreover, the effects of boundary conditions should be included to evaluate the uncertainty precisely, whereby it may quantify the GV zones near the boundaries. Accordingly, the validation of GV should deal with the data of geochemistry related to human activities to enhance the correlation of the GV map and contaminants.