Changes in sea-level, precipitation patterns, and water usage impact groundwater flow and seawater intrusion in coastal aquifers. Reduced-order models are favored for rapid, large-scale simulations as they are based on relevant physical processes, but avoid overly complex, data intensive, and computationally expensive techniques. Reduced-order models are intended for screening-level assessments of seawater intrusion impacts to freshwater resources and to quickly assess future impacts of seawater intrusion due to sea-level rise. A popular reduced-order model for coastal aquifers is an interface flow model, in which an interface separates freshwater from saltwater. An interface formulation is capable of simulating the movement of the freshwater, the saltwater, and the elevation of the interface with much lower computational effort than variable-density models for coupled flow and transport. We developed a new control-volume finite-difference formulation based on a water balance of both the freshwater zone and the saltwater zone in combination with a continuity condition at the interface. The numerical formulation is fully implicit, uses upstream weighting of the freshwater and saltwater thicknesses, and solves for the freshwater head, and (optionally) for the saltwater head simultaneously using Newton’s method. As the formulation is fully implicit, it allows for the use of large time steps to quickly compute the final steady state position of the interface, which is useful for many applications. Several example models will be presented simulating transient interface movement for advancing and receding interfaces in multi-layer systems, including a comparison to the results of a variable density model. The formulation has been implemented in Python (for multi-layer cross-sectional simulations) and is being implemented in a new MODFLOW 6 package and in Aquaveo’s Groundwater Modeling System (GMS) graphical user interface.