A stochastic inverse algorithm to jointly analyze multiple geophysical and hydrological datasets for a geothermal prospect is presented. The purpose is to improve prospect evaluation and estimate the likelihood of useful temperature and fluid flow fields at depth. The approach combines Bayesian inference with a Markov Chain Monte Carlo (MCMC) global search algorithm to conduct a staged, model-based inversion of the different data sets. Initial estimates of uncertainty in structural or parametric characteristics of the prospect are used to drive large numbers of simulations of hydrothermal fluid flow and related geophysical processes using random realizations of the conceptual geothermal system. The results consist of a subset of these realizations - an equally probable ensemble of alternatives - that best match the observed datasets within a specified norm or tolerance. Statistical (posterior) characteristics of these solutions reflect reductions in the perceived (prior) uncertainties. The method is highly flexible and capable of accommodating multiple and diverse datasets as a means to maximize the utility of all available data to understand system behavior.