A distributed-parameter physically-based solute transport model using a novel approach to describe surface-subsurface interactions is coupled to an existing flow model. In the integrated model the same surface routing and mass transport equations are used for both hillslope and channel processes, but with different parametrizations for these two cases. For the subsurface an advanced time-splitting procedure is used to solve the advection-dispersion equation for transport and a standard finite element scheme is used to solve Richards equation for flow. The surface-subsurface interactions are resolved using a mass balance-based surface boundary condition switching algorithm that partitions water and solute into actual fluxes across the land surface and changes in water and mass storage. The time stepping strategy allows the different time scales that characterize surface and subsurface water and solute dynamics to be efficiently and accurately captured. The model features and performance are demonstrated in a series of numerical experiments of hillslope drainage and runoff generation