Borehole data are essential for calibrating ice sheet models, but field expeditions are often time-consuming, expensive, and dangerous. It is thus important to identify drilling locations that maximize the value of data while minimizing costs and risks. We employ optimal experimental design (OED) to select locations that yield the highest expected information gain subject to site accessibility. Simulation-based design requires evaluating a computationally expensive and nonlinear ice sheet model, and using naive nested Monte Carlo to calculate expected information gain would be intractable. We thus target an upper bound of the expected utility, and accelerate its computation via a polynomial surrogate function for the Bayesian evidence. This surrogate is locally adaptive, where additional simulations are performed in inaccurate neighborhoods as assessed by an error indicator. Furthermore, a bias-variance trade-off guarantees bias introduced by the surrogate model to be the same order of magnitude as the Monte Carlo variance. The borehole data are subsequently used for Bayesian inference of past surface temperatures via efficient Markov chain Monte Carlo.