In the present study, an enhanced subsurface prediction algorithm based on a non-parametric geostatistical model and a history matching technique through Gibbs sampler is developed and the iterative prediction improvement procedure is proposed. The developed model is applied to a simple two-dimensional synthetic case where domain is composed of three different hydrogeologic media with $500m{ imes}40m$ scale. In the application, it is assumed that there are 4 independent pumping tests performed at different vertical interval and the history curves are acquired through numerical modeling. With two hypothetical borehole information and pumping test data, the proposed prediction model is applied iteratively and continuous improvements of the predictions with reduced uncertainties of the media distribution are observed. From the results and the qualitative/quantitative analysis, it is concluded that the proposed model is good for the subsurface prediction improvements where the history data is available as a supportive information. Once the proposed model be a matured technique, it is believed that the model can be applied to many groundwater, geothermal, gas and oil problems with conventional fluid flow simulators. However, the overall development is still in its preliminary step and further considerations needs to be incorporated to be a viable and practical prediction technique including multi-dimensional verifications, global optimization, etc. which have not been resolved in the present study.