Surface water, the vadose zone, and groundwater are linked components of a hydrologic continuum. In order to capture the interaction between different components of a hydrologic continuum and to use this understanding in water management situations, an accurate numerical model is needed. The quality of model results depends on accurate representation of the physical processes and the data describing the area of interest, as well as performance of the numerical formulation implemented. Here we present a physics-based, distributed, fully coupled, second-order accurate, upwind cell-centered, constrained unstructured mesh based finite-volume modeling framework (FIHM) that simultaneously solves two-dimensional unsteady overland flow and three-dimensional variably saturated subsurface flow in heterogeneous, anisotropic domains. A multidimensional linear reconstruction of the hydraulic gradients (surface and subsurface) is used to achieve second-order accuracy. Accuracy and efficiency in raster data and vector-boundary representations are facilitated through the use of constrained Delaunay meshes in domain discretization. The experiments presented here (i) explore the influence of initial moisture conditions, soil properties, anisotropy, and heterogeneity in determining the pressure head distributions in the vadose and saturated zones, (ii) show the existence of localized “flux rotation” phenomenon due to heterogeneous anisotropy, leading to the creation of convergence–divergence zones, (iii) show the influence of vertical drainage from unsaturated zone on the response of an unconfined aquifer to pumping, and (iv) show the effects of capillarity, saturation excess, infiltration excess, and initial water table location on determining the overland flow generation.