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. ¿ American Geophysical Union 1995