Problems in hydrology and water management that involve both surface water and groundwater are best addressed with simulation models that can represent the interactions between these two flow regimes. In the current generation of coupled models, a variety of approaches is used to resolve surface-subsurface interactions and other key processes such as surface flow propagation. In this study we compare two physics-based numerical models that use a 3D Richards equation representation of subsurface flow. In one model, surface flow is represented by a fully 2D kinematic approximation to the Saint-Venant equations with a sheet flow conceptualization. In the second model, surface routing is performed via a quasi-2D diffusive formulation and surface runoff follows a rill flow conceptualization. The coupling between the land surface and the subsurface is handled via an explicit exchange term resolved by continuity principles in the first model (a fully-coupled approach) and by special treatment of atmospheric boundary conditions in the second (a sequential approach). Despite the significant differences in formulation between the two models, we found them to be in good agreement for the simulation experiments conducted. In these numerical tests, on a sloping plane and a tilted V-catchment, we examined saturation excess and infiltration excess runoff production under homogeneous and heterogeneous conditions, the dynamics of the return flow process, the differences in hydrologic response under rill flow and sheet flow parameterizations, and the effects of factors such as grid discretization, time step size, and slope angle. Low sensitivity to vertical discretization and time step size was found for the two models under saturation excess and homogeneous conditions. Larger sensitivity and differences in response were observed under infiltration excess and heterogeneous conditions, due to the different coupling approaches and spatial discretization schemes used in the two models. For these cases, the sensitivity to vertical and temporal resolution was greatest for processes such as reinfiltration and ponding, although the differences between the hydrographs of the two models decreased as mesh and step size were progressively refined. In return flow behavior, the models are in general agreement, with the largest discrepancies, during the recession phase, attributable to the different parameterizations of diffusion in the surface water propagation schemes. Our results also show that under equivalent parameterizations, the rill and sheet flow conceptualizations used in the two models produce very similar responses in terms of hydrograph shape and flow depth distribution.