This paper introduces a new code “QuasIncompact3D” for solving the variabledensity Navier-Stokes equations in the low-Mach number limit. It is derived from the Incompact3D framework which is designed for incompressible flows [1]. QuasIncompact3D is based on high-order accurate compact finite-differences [2], an efficient 2D domain decomposition [3] and a spectral Poisson solver. The first half of the paper focuses on introducing the low-Mach number governing equations, the numerical methods and the algorithm employed by QuasIncompact3D to solve them. Two approaches to forming the pressure-Poisson equation are presented: one based on an extrapolation that is efficient but limited to low density ratios and another one using an iterative approach suitable for higher density ratios. The scalability of QuasIncompact3D is demonstrated on several TIER-1/0 supercomputers using both approaches, showing good scaling up to 65k cores. Validations for incompressible and variable-density low-Mach number flows us-ing the Taylor-Green vortex and a non-isothermal mixing layer, respectively, as test cases are then presented, followed by simulations of non-Boussinesq gravity currents in two- and three-dimensions. To the authors’ knowledge this is the first investigation of 3D non-Boussinesq gravity currents by means of Direct Numerical Simulation over a relatively long time evolution. It is found that 2D and 3D simulations of gravity currents show differences in the locations of the fronts, specifically that the fronts travel faster in three dimensions, but that it only becomes apparent after the initial stages. Our results also show that the difference in terms of front location decreases the further the flow is from Boussinesq conditions.