A new hybrid model for the dynamics of glaciers, ice sheets and ice shelves is introduced. In this ‘multilayer’ model the domain of ice consists of a pile of thin layers, which can spread out, contract and slide over each other such that the two most relevant types of stresses are accounted for: membrane and vertical shear. Assuming the horizontal velocity field to be vertically piecewise constant in each layer, the model is derived from local depth integrations of the hydrostatic approximation of the Stokes equations. These integrations give rise to interlayer tractions, which can be redefined at zeroth order in the interlayer surface slope by keeping the vertical shear stress components. Furthermore, if the layers are chosen such that they are aligned with the streamlines, then second-order accurate interlayer tractions can replace zeroth-order ones. The final model consists of a tridiagonal system of two-dimensional nonlinear elliptic equations, the size of this system being equal to the number of layers. When running the model for prognostic flowline ISMIP-HOM benchmark experiments, the multilayer solutions show good agreement with the higher-order solutions if no severe depression occurs in the bedrock. As an alternative to three-dimensional models, the multilayer approach offers to glacier and ice sheet modellers a way of upgrading the commonly used shallow shelf approximation model into a mechanically complete but mathematically two-dimensional model.