We present a Bayesian hierarchical modeling approach to paleoclimate reconstruction using borehole temperature profiles. The approach relies on modeling heat conduction in solids via the heat equation with step function, surface boundary conditions. Our analysis includes model error and assumes that the boundary conditions are random processes. The formulation also enables separation of measurement error and model error. We apply the analysis to data from nine borehole temperature records from the San Rafael region in Utah. We produce ground surface temperature histories with uncertainty estimates for the past 400 years. We pay special attention to use of prior parameter models that illustrate borrowing strength in a combined analysis for all nine boreholes. In addition, we review selected sensitivity analyses.