Inverse problems are ubiquitous in science and engineering, especially in the field of Geoscience. Algorithms that we have developed are applicable to seismic imaging, contaminant source identification, hydraulic tomography, etc. This work tackles these inverse problems using the Geostatistical approach that stochastically models unknowns as random fields. However, due to high computational costs in identifying small scale features, these methods are challenging. In addition, it is necessary to quantify the corresponding predictive uncertainty to provide a sound basis for management or policy decision making. My approach uses Hierarchical matrices to efficiently represent dense covariance matrices and solves the resulting intermediary system of equations using preconditioned iterative methods. The resulting cost is reduced from O(N^2) to O(N log N), where N is the number of unknowns to be determined. In addition, I will describe some recent work on Krylov solvers for shifted systems that accelerates the solution of "forward problems" and as a result, the inverse problem. Moreover, a technique to perform uncertainty quantification by estimating the diagonals of the posterior covariance matrix will be presented. I will illustrate these algorithms with examples from contaminant source identification and hydraulic tomography.