Discontinuous Galerkin (DG) methods due to their robustness properties, e.g. local conservation, low numerical dispersion, and well-capturing strong shocks and physical discontinuities, are well-suited for the simulation of Variable Density Flow (VDF) in porous media. This paper aims at introducing, in a unified format, the general class of Interior Penalty DG (IPDG) methods to solve the VDF equations. A combination of symmetric, non-symmetric and incomplete IPDG methods is used to discretize both head and concentration variables. Compatibility analysis is performed to prevent the loss of accuracy of the IPDG methods in simulations of coupled flow and transport equations. An accurate technique is used for time integration, based on a non-iterative procedure and adaptive time stepping with embedded error control. Several benchmarks are investigated to validate the proposed DG scheme and to examine its performance in simulating VDF problems. The new DG scheme reproduces better the experimental data than the conventional SEAWAT model. Its results are in excellent agreement with a recent semi-analytical solution of the Henry problem, dealing with seawater intrusion under convection-dominating conditions. The performance of the DG scheme is examined by simulating the challenging problem of natural convection in porous enclosure. The method is compared against a finite element solution obtained with COMSOL multi-physics. The numerical experiments indicate clearly that high-order DG method is much more appropriate than standard conforming Galerkin method in simulating VDF problems while at the same time, guaranteeing a better precision and high-fidelity solutions. The proposed numerical method can be extended to 3D problems.