A sharp interface numerical model is developed to simulate saltwater intrusion in multilayered coastal aquifer systems. The model takes into account the flow dynamics of salt water and fresh water assuming a sharp interface between the two liquids. In contrast to previous two‐fluid flow models which were formulated using the hydraulic heads of fresh water and salt water as the dependent variables, the present model employs a mixed formulation having one fluid potential and a pseudosaturation as the dual dependent variables. Conversion of the usual sharp interface flow equations for each aquifer to an equivalent set of two‐phase flow equations leads to the definitions of pseudosaturation, capillary pressure, and constitutive relations. The desired governing equations are then obtained by connecting neighboring aquifers via vertical leakage. The proposed formulation is based on a Galerkin finite element discretization. The numerical solution incorporates upstream weighting and nonlinear algorithms with several enhanced features, including rigorous treatment of aquitard leakage and well conditions, and a robust Newton‐Raphson procedure with automatic time stepping. The present sharp interface numerical model is verified using three test problems involving unconfined, confined, and multilayered aquifer systems and consideration of steady state and transient flow situations. Comparisons of numerical and analytical solutions indicate that the numerical schemes are efficient and accurate in tracking the location, lateral movement, and upconing of the freshwater‐saltwater in