We developed a Bayesian model to invert magnetotelluric (MT) data using a 2D sharp boundary parameterization. We divided the 2D cross section into layers and considered the locations of interfaces and resistivity of the regions formed by the interfaces as random variables. We assumed that those variables are independent in the vertical direction and dependent along the lateral direction, whose spatial dependence is described by either pairwise difference or multivariate Gaussian priors. We used a parallel, adaptive finite-element algorithm to rapidly forward simulate frequency-domain MT responses of the 2D resistivity structure and used Markov chain Monte Carlo methods to draw many samples from the joint posterior probability distribution. We applied the Bayesian model to a synthetic case that mimics a geothermal exploration scenario. Our results demonstrated that the developed method is effective in estimating the resistivity and depths to interfaces and in quantifying uncertainty on the estimates. We also applied the developed method to the field MT data collected from the Darajat geothermal site in Indonesia. We compared our inversion results with those obtained from a deterministic inversion of 3D MT data; they are consistent even if the two inversion methods are very different and the amount of information used for inversion is different.