A new formulation is presented for the modeling of seawater intrusion in coastal multiaquifer systems. The formulation is intended for the simulation of three-dimensional regional seawater intrusion starting from a specified salinity distribution. The Dupuit approximation is adopted, and diffusion and dispersion are not taken into account. The formulation is based on a vertical discretization of the groundwater into zones of either constant density (stratified flow) or continuously varying density (piecewise linear in the vertical direction). The zones are separated from each other by curved surfaces, representing interfaces for stratified flow and isosurfaces of the density for variable density flow. During a simulation both the change of the elevations of the surfaces and the change in head are computed through consistent application of continuity of flow; a simple tip and toe tracking algorithm is applied to simulate the horizontal movement of the surfaces. The resulting differential equations are equivalent in form to the differential equation for single-density flow. The density effects are incorporated through the addition of pseudosource terms. This makes it easy to implement the formulation in existing numerical groundwater codes. The main advantage of the formulation is the tremendous reduction of the number of cells needed for a simulation because every aquifer may be represented by a single layer of cells. The accuracy of the formulation is demonstrated through comparison with existing exact solutions and a numerical solution.