We present a mathematical and numerical model for coupled subsurface and surface flows based on coupling the Stokes and the Darcy equations through the Beavers-Joseph-Saffman interface conditions. Optimal order error estimates are established for a finite element discretization based on conforming Stokes elements in the surface flow domain and mixed finite elements in the porous media domain. The formulation utilizes a Lagrange multiplier to impose the interface conditions. A non-overlapping domain decomposition algorithm is developed which reduces the coupled algebraic system to an interface problem for the normal stress. Each interface iteration requires solving Stokes and Darcy subdomain problems. It is shown that the interface problem is symmetric and positive definite and that its condition number is $O(1/h)$, where $h$ is the discretization parameter. Numerical results and parallel scalability studies are presented.