Enhanced water management systems depend on accurate estimation of hydraulic properties of subsurface formations. This is while hydraulic conductivity of geologic formations could vary significantly. Therefore, using information only from widely spaced boreholes will be insufficient in characterizing subsurface aquifer properties. Hence, there is a need for other sources of information to complement our hydro-geophysics understanding of a region of interest. This study presents a numerical framework where information from different measurement sources is combined to characterize the 3-dimensional random field representing the hydraulic conductivity of a watershed in a Multi-Fidelity estimation model. Coupled with this model, a Bayesian experimental design will also be presented that is used to select the best future sampling locations. This work draws upon unique capabilities of electrical resistivity tests as well as statistical inversion. It presents a Multi-Fidelity Gaussian Processes (Kriging) model to estimate the geological properties in Upper Sangamon Watershed in east central Illinois, using multi-source observation data, obtained from electrical resistivity and pumping tests. We demonstrate the accuracy of Co-Kriging that is dependent on the locations and the distribution of both the high- and low-fidelity data, and also discuss its comparison with Single-High-Fidelity Kriging results. The uncertainties and confidence in the measurements and parameter estimates are then quantified and are in turn used to design future cycles of data collection to further improve the confidence intervals.