The delineation of well capture zones is of great importance to accurately define well head protection area (WHPA) for potential groundwater resources and the public security. Natural aquifer systems typically involve different extent of heterogeneity in aquifer parameters and such parameter variations can directly influence the estimations of flow fields and the delineations of WHPAs. This study employs an unconditional approximate spectral method (ASM) associated with backward particle tracking algorithm to delineate stochastic well capture zones in aquifers of Choushui River alluvial fan (CRAF) in central Taiwan. The analysis integrates hourly-recorded groundwater observations from 1995 to 2013 to be the mean flow field. We implement the developed model to 187 Taiwan Water Corporation (TWC) wells for domestic water supplies in CRAF. With predefined small-scale heterogeneity for hydraulic conductivity, the uncertainty of capture zones are obtained based on the observed pumping rates at TWC wells. Results of the analyses show that the average distances of mean capture zones in the first layer of CRAF are about one kilometer from the TWC wells. The small-scale hydraulic conductivity can induce capture zone uncertainties ranging from meters to tens of meters in one year depending on the complexity of the flow field. The uncertainty zones of WHPA in CRAF can be served as the basis to conduct risk analysis for drinking water.