This study presents an approach for obtaining limited sets of realizations of hydraulic conductivity (K) of multiple aquifers using simulated annealing (SA) simulation and spatial correlations among aquifers to simulate realizations of hydraulic heads and quantify their uncertainty in the Pingtung Plain, Taiwan. The proposed approach used the SA algorithm to generate large sets of natural logarithm hydraulic conductivity (ln(K)) realizations in each aquifer based on spatial correlations among aquifers. Moreover, small sets of ln(K) realizations were obtained from large sets of realizations by ranking the differences among cross-variograms derived from the measured ln(K) and the simulated ln(K) realizations between the aquifer pair Aquifer 1 and Aquifer 2 (hereafter referred to as Aquifers 1-2) and the aquifer pair Aquifer 2 and Aquifer 3 (hereafter referred to as Aquifers 2-3), respectively. Additionally, the small sets of realizations of the hydraulic conductivities honored the horizontal spatial variability and distributions of the hydraulic conductivities among aquifers to model groundwater precisely. The uncertainty analysis of the 100 combinations of simulated realizations of hydraulic conductivity was successfully conducted with generalized likelihood uncertainty estimation (GLUE). The GLUE results indicated that the proposed approach could minimize simulation iterations and uncertainty, successfully achieve behavioral simulations when reduced between calibration and evaluation runs, and could be effectively applied to evaluate uncertainty in hydrogeological properties and groundwater modeling, particularly in those cases which lack three-dimensional data sets yet have high heterogeneity in vertical hydraulic conductivities.