The rugged topography of mountain ranges represents a special challenge to computational ice sheet models simulating past or present glaciations. In mountainous regions, the topography steers glaciers through relatively narrow and steep valleys, and as a consequence hereof, the flow rate of alpine‐style glaciers varies significantly at length scales comparable to that of the topography. Localized flow rate undulations generate longitudinal and transverse stress gradients within the ice, which, in turn, are of known importance to the flow itself, whether by internal deformation of the ice or by basal sliding, and to the interaction with topography through glacial erosion. Standard models capable of simulating mountain range–scale glaciations build on the so‐called shallow ice approximation, which, among other parameters, neglects the longitudinal and transverse stress gradients, and therefore fails to capture the full effects of the rugged topography and related feedbacks between erosion by glacial sliding and the extent and style of glaciation. Here we present and test a new depth‐integrated model framework which, on the one hand, takes into account the “higher‐order” effects related to steep and rugged bed topography while it still, on the other hand, provides the computational efficiency needed for three‐dimensional simulations of glaciation and landscape evolution in response to, for example, long‐term climate variations. The model framework (integrated second‐order shallow ice approximation (iSOSIA)) is based on depth integration of the second‐order shallow ice approximation. Although depth integration is not guaranteed to maintain so‐called second‐order accuracy, the results of computational benchmark experiments show that iSOSIA, in spite of its efficient depth‐integrated structure, performs well in situations with steep and rapidly undulating bed topography.