Shallow water equations are widely used for the simulation of surface water flows. To simulate groundwater flows one might consider the Darcy equations in the saturated case and the nonlinear Richards equation in the unsaturated case. The flows are coupled along an interface where we enforce continuity of pressure and continuity of flux. These assumptions lead to a nonlinear system of equations where the nonlinearities lie in the advection term of the shallow water equations and the nonlinear Richards equation. The nonlinearity in the advection term can be treated by the method of characteristics and an implicit time integration scheme whereas the nonlinearity in the Richards equation can be treated by a Newton method. We will utilize a least squares finite element method for the spatial discretization where we will make use of continuous piecewise polynomials combined with Raviart-Thomas elements. At the end of the talk numerical experiments will be given and an adaptive scheme, using the least squares functional as an error estimator, will be examined.