The majority of inversion methods used for inferring past ground surface temperatures (GST) from borehole temperature-depth profiles rely on the assumption that heat flow is in the vertical direction only. This means that accounting for certain effects caused by the local terrain of a borehole is not possible and consequently, many borehole profiles cannot be used with confidence. Here, we describe a methodology to avoid this problem by solving the heat conduction forward problem in 3-D using finite elements (FE). In order to make the inversion approach computationally tractable, we reduce the dimensions of this FE model using proper orthogonal decomposition. The inverse problem is cast in a probabilistic Bayesian framework for which the posterior probability distribution of the past GSTs is sampled using a reversible jump Markov chain Monte Carlo algorithm. This allows the resolution of the GST history over time to be explored by varying the parameterization of the GST model. Synthetic examples calculated with moderate topographies demonstrate the efficacy of the Bayesian 3-D inversion method, and the results are compared with those using a 1-D approach. For moderate topography, the latter can lead to spurious GST reconstructions. A further synthetic example demonstrates that the effect of incorrectly assuming lateral geological homogeneity is negligible. The inversion method is also compared with a more standard inversion method. A significant advantage of the Bayesian approach is that uncertainties in all of the model parameters can be accounted for, leading to a more realistic interpretation of the range of GST histories supported by the data. The methods presented here should allow a broader range of geothermal data to be used for paleoclimate reconstruction purposes in the future.