A mathematical framework is presented for three‐dimensional shallow groundwater flow with variable density. The formulation is based on the Dupuit‐Forchheimer assumption. The problem is posed in terms of a discharge potential that satisfies the same differential equation as the discharge potentials for single‐density flow. The freshwater head, defined as the pressure divided by the unit weight of fresh water plus the elevation head, may be computed as a function of position in three dimensions from the potential and a known density distribution. The density distribution may be approximated using the multiquadrics interpolator. It is explained how the change in density may be computed as a function of time. Discontinuities in the aquifer properties cause a jump in the normal component of flow for flow fields computed with the Dupuit‐Forchheimer approximation. An interpretation of this jump is given by comparison with an exact formulation, which makes it possible to obtain the approximate streamlines as they cross discontinuities.