This research develops a zonation-kriging (ZK) method and extends the geostatistical theory to investigate the spatially distributed random field and obtain an optimal nonsmooth conditional estimate. In the ZK method the zonal structure honors the parameter measurements such that each sample location represents a distinct zone, while the kriging method also honors the same data set and results in a smooth distribution. The ZK method integrates the conditional estimates of the zonal structure and the kriged field over a set of weighting coefficients on the sample data. The extended geostatistical theory using the ZK method is applied to the intrinsic field, stationary field, and nonstationary field. The study shows that the conditional variance using the ZK method is bounded between the kriging variance and variogram. When estimating hydraulic conductivity by inversion of a groundwater model, the optimal values of weighting coefficients are determined via predictive probability of the Bayesian decision theory. The numerical examples show that the ZK method has less dependence on the semivariogram model. The ZK method outperforms the kriging method in both smooth and nonsmooth fields. The ZK method shows high efficiency to better the conditional estimate of hydraulic conductivity distribution with low conditional uncertainty and good agreement to the groundwater head observations.