We evaluate an inverse approach to reconstruct distributed bedrock topography and simultaneously initialize an ice flow model. The inverse method involves an iterative procedure in which an ice dynamical model (PISM) is run multiple times over a prescribed period, while being forced with space-and time-dependent climate input. After every iteration bed heights are adjusted using information of the remaining misfit between observed and modeled surface topography. The inverse method is first applied in synthetic experiments with a constant climate forcing to verify convergence and robustness of the approach in three dimensions. In a next step, the inverse approach is applied to Nordenskioldbreen, Svalbard, forced with height-and time-dependent climate input since 1300 AD. An L-curve stopping criterion is used to prevent overfitting. Validation against radar data reveals a high correlation (up to R = 0.89) between modeled and observed thicknesses. Remaining uncertainties can mainly be ascribed to inaccurate model physics, in particular, uncertainty in the description of sliding. Results demonstrate the applicability of this inverse method to reconstruct the ice thickness distribution of glaciers and ice caps. In addition to reconstructing bedrock topography, the method provides a direct tool to initialize ice flow models for forecasting experiments.