While the magnetotelluric (MT) method has been used in geothermal exploration for many years, the integrated evaluation of the effects of reservoir temperature, clay alteration, salinity and permeability on electrical resistivity have remained highly qualitative. In order to contribute to quantifying these effects, we develop an inversion methodology that integrates methylene blue (MeB) logs from wells with the electrical resistivity distribution from MT surveys. We use this method to infer, under uncertainty, clay cap boundaries in a geothermal field. Our methodology is based on a fast stochastic inversion of MT signals using Markov Chain Monte Carlo (MCMC) to fit a one-dimensional three-layer resistivity model beneath each MT station, laterally constrained by MeB data. The advantage of our approach is the explicit investigation and visualization of inversion uncertainty, yielding inferred depth intervals for clay cap boundaries. The methodology is tested on a profile of MT stations and MeB logs in wells located in the western area of the Wairākei geothermal field in the Taupō Volcanic Zone, New Zealand. Over the long-term, our goal is to improve our knowledge of uncertainty related to MT signals and its graphical representation to enhance our understanding of the relationship between resistivity and the clay cap.