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%–20% of the total thickness, leading to not insignificant variation in horizontal density. Nevertheless, in most vertically integrated ice-flow models today, the density is assumed 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, which reveal significant qualitative differences between each of the proposed density formulations. We recommend always including the horizontal gradient of the density in the body-force term of the momentum equation, and in the mass-conservationequation. Furthermore, by modeling the transient evolution of a large sector of the West Antarctic Ice Sheet (WAIS), we find that doing so leads to about a 10% correction in the estimated change in volume above flotation over 40 years. We conclude that including HDVs in flow modeling of the AIS is important for accurate predictions of mass loss.