AbstractModelling paleo-glacier networks in mountain ranges on the millennial timescales requires ice flow approximations. Hybrid models calculating ice flow by combining vertical shearing (shallow ice approximation) and longitudinal stretching (shallow shelf approximation) have been applied to model paleo-glacier networks on steep terrain, yet their validity has not yet been assessed quantitatively. Moreover, hybrid models consistently yield higher ice thicknesses than Last Glacial Maximum geomorphological reconstructions in the European Alps. Here, we compare results based on the hybrid Parallel Ice Sheet Model (PISM) and the Stokes model Elmer/Ice on the Rhine Glacier, a catchment of the former European Alpine Icefield. For PISM, we also test two magnitudes of flux limitation in a scheme that reduces shearing velocities. We find that the flux limitation typically used in PISM yields significantly reduced shearing speeds and increases ice thicknesses by up to 500 m, partly explaining previous overestimations. However, reducing the ice flux limitation allows the hybrid model to minimize this mismatch and captures sliding speeds, ice thicknesses, ice extent and basal temperatures in close agreement with those obtained with the Stokes model.