This study proposes a new inverse method to estimate the bed topography elevation beneath glaciers flows from surface observations (altimetry elevations and InSAR velocities) and (very) sparse depth measurements (e.g. acquired during airborne campaigns). To do so an original form of depth-integrated flow equations (long-wave assumption) is derived. The latter are valid for highly to mildly-sheared regimes hence including moderately fast flows; varying internal thermal profiles are taken into account. The inverse problem is particularly challenging since the surface signatures integrate the bottom features (bed elevation and friction–slip amount) plus the internal deformation. The first key ingredient of the inverse method is the derivation of this non-isothermal Reduced Uncertainty (RU) version of the classical SIA equation (lubrication type model for generalised Newtonian fluids) by intrinsically integrating the surface measurements in the formulation. This resulting multi-physics RU-SIA equation contains a unique uncertain dimensionless parameter only (the parameter γ). The next key ingredient is an advanced Variational Data Assimilation (VDA) formulation combined with a purely data-driven extension of γ based on the trend observed in the (sparse) depth measurements (e.g. along the flight tracks). The resulting inverse method provides the first physical-based depth estimations in mildly-sheared mildly-slippery shallow flows. In poorly monitored ice-sheet areas (e.g. in Antarctica), the resulting estimations are noticeably less uncertain than the current ones (in particular compared to those obtained by gravimetry inversions). The present numerical experiments and experimental sensitivity analyses demonstrate the reliability of this new RU-SIA equation and the robustness of the inverse method. In other respect, the RU-SIA equation may provide a-posteriori estimations of the thermal boundary layer at bottom.