This study evaluated the spatial variability of hydraulic conductivity (K) along the Han River watershed in the Jianghan Plain (China) by using different geostatistical methods. The K distribution was estimated from 58 borehole measurements, and co-kriging estimates were obtained by incorporating 287 electrical resistivity (ρ) and 260 hydraulic gradient (I) values as auxiliary data. Additionally, 2,980 lithology samples were classified into three hydrofacies units in terms of decreasing K, which was simulated accordingly. The results indicate that the highest K values are in the western part of Shahu town and near the Yangtze River. The trend continues to diminish across western and northern Xiantao town and north-west of Xingou town, as well as across the transition zone from highland to lowland (from Shayang town to Zhanggang town). The low K values occur from the northern part of Jiukou town to Shayang town, and around Xingou town, and then are scattered around the high K values, with the lowest K values distributed in the north-west of the study area (from the catchment margin to the north of Shayang town), in Yanglinwei town, and east of Shahu town. As for the geostatistical methods, co-kriging estimates based on ρ and I highlight the variability of K, with the latter providing more realistic results. The three-dimensional hydrofacies model realistically represents the spatial distributions of K in the Quaternary aquifer system. The estimated results and hydrofacies simulation could provide a basic database for further groundwater modeling and water resource evaluation.