We present a shape optimization algorithm to estimate the ice thickness distribution within a two-dimensional, non-sliding mountain glacier, given a transient surface geometry and a mass-balance distribution. The approach is based on the minimization of the surface topography misfit at the end of the glacier's evolution in the shallow ice approximation of ice flow. Neither filtering of the surface topography where its gradient vanishes nor interpolation of the basal shear stress is involved. Novelty of the presented shape optimization algorithm is the use of surface topography and mass-balance only within a time-dependent Lagrangian approach for moving-boundary glaciers. On real-world inspired geometries, it is shown to produce estimations of even better quality in smaller time than the recently proposed steady and transient inverse methods. A sensitivity analysis completes the study and evinces the method's higher susceptibility to perturbations in the surface topography than in surface mass-balance or rate factor.