The distribution of hydrogeological materials in a three-dimensional heterogeneous aquifer system has a large effect on groundwater flow and land subsidence simulations. The hydrogeological information for regions between boreholes embeds a large amount of uncertainty into the hydrogeological model, and thus affects numerical assessment. Quantifying the effects of the heterogeneous system and hydrogeological model uncertainty on groundwater flow and land subsidence simulations is thus important. Here, data from 46 geological boreholes in Huwei Town, Taiwan, were adopted to investigate hydrogeological model uncertainty. The one-dimensional continuous-lag Markov chain and the geostatistical method were used to analyze the spatial characteristics of hydrogeological materials and generate realizations of the hydrogeological model based on the assessment results. Estimated hydrological conditions and hydraulic parameters were applied to mitigate uncertainty not caused by the hydrogeological model. The mean of the results of land subsidence in Monte Carlo simulations showed a more stable distribution than that of individual realizations, for which land subsidence might be far from the mean. Therefore, the mean and variance results provide a reliable assessment with uncertainty information for land subsidence simulations. The coefficient of variation (CV) was used to quantify the hydrogeological model uncertainty. The CV value of land subsidence was larger than that of the hydraulic head because the thickness of clay is not uniform whereas the flow pattern is smooth. The variance and CV distributions of the hydraulic head and land subsidence provide uncertainty information that can be used to guide site investigations. A heterogeneous hydrogeological model with uncertainty quantification should be carefully applied to land subsidence simulations to obtain reasonable results.