We develop and analyze a mixed finite element method for the solution of an elliptic system modeling a porous medium with large cavities, called vugs. It consists of a second-order elliptic (i.e., Darcy) equation on part of the domain coupled to a Stokes equation on the rest of the domain, and a slip boundary condition (due to Beavers–Joseph–Saffman) on the interface between them. The tangential velocity is not continuous on the interface. We consider a 2-D vuggy porous medium with many small cavities throughout its extent, so the interface is not isolated. We use a certain conforming Stokes element on rectangles, slightly modified near the interface to account for the tangential discontinuity. This gives a mixed finite element method for the entire Darcy–Stokes system with a regular sparsity pattern that is easy to implement, independent of the vug geometry, as long as it aligns with the grid. We prove optimal global first-order L 2 convergence of the velocity and pressure, as well as the velocity gradient in the Stokes domain. Numerical results verify these rates of convergence and even suggest somewhat better convergence in certain situations. Finally, we present a lower dimensional space that uses Raviart–Thomas elements in the Darcy domain and uses our new modified elements near the interface in transition to the Stokes elements.