Landslides are hazardous phenomena whose analysis requires accurate modeling of the failure surface geometry in order to perform stability computations, risk quantification and develop predictive models. However, direct measurements are costly and classical approaches often rely on oversimplified geometric assumptions. Here, the landslide thickness estimation problem is formulated as an ill-posed inverse problem deriving from the discretisation of the two-dimensional mass conservation equation. To stabilise and solve the inverse problem, we employ a multipenalty regularisation approach based on the decomposition of the landslide domain in several regions with a uniform magnitude of the horizontal velocity. We develop an extension of the Balancing Principle to simultaneously identify the values of the regularisation parameters across the decomposed regions and evaluate the landslide thickness. A convergent iterative method is proposed for the numerical realisation of the Balancing Principle. Numerical tests on synthetic and real-world landslide datasets confirm that decomposing the domain according to velocity field homogeneity improves estimation accuracy. Our results indicate that while excessive domain partitioning may degrade performance, identifying suitably defined macro-regions based on velocity characteristics leads to optimal estimation outcomes. This work provides a physically grounded and mathematically solid baseline for landslide thickness estimation, laying a robust foundation for future integration with data-driven approaches and hybrid modeling strategies.