Gravity-driven flow of large ice masses such as the Antarctic Ice Sheet (AIS) depends on both the geometry and the mass density of the ice sheet. The vertical density profile can be approximated as pure ice overlain by a firn layer of varying thickness, and for the AIS the firn thickness is not uncommonly 10 to 20% of the total thickness, leading to not insignificant variation in density. Nevertheless, in most vertically-integrated ice-flow models today the density is assumed to be constant, sometimes with an adjustment in thickness to compensate. In this study, we explore the treatment of horizontal density variations (HDVs) within vertically-integrated ice-sheet models. We assess the relative merits and shortcomings of previously proposed approaches, and provide new formulations for including HDVs. We use perturbation analysis to derive analytical solutions that describe the impact of density variations on ice flow for both grounded ice and floating ice shelves. Our analytical solutions reveal significant qualitative differences between each of the proposed density formulations. Furthermore, by modelling the transient evolution of a large sector of the West Antarctic Ice Sheet (WAIS), we quantify the impact of HDVs on estimated sea level change. For WAIS we find that explicitly including the horizontal density gradients in the momentum and mass conservation equations leads to about a 10% correction in the estimated change in volume above flotation over 40 years. We conclude that including horizontal density variations in flow modelling of the Antarctic Ice Sheet is important for accurate predictions of mass loss.