Existing simulation techniques for subsurface water quality management remain inadequate despite the phenomenal progress made in the computing capabilities in the past 2 decades. One major lacunas found in mathematical models for groundwater hydrodynamics in combined free flow zones (zones that are devoid of porous materials) and porous domains (e.g., soil water regions) in the subsurface is that these models impose continuity of all flow/transport parameters to describe the flow behavior at the intermediate surfaces between the two zones. Prescription of such interfacial conditions cannot always be physically justified and, in general, result in ill-posed mathematical formulations. The present paper aims to remove the above limitations of the mathematical models by imposing well-posed mathematical formulations of mass and momentum transfer across the boundaries between free and porous zones underground for simulating subsurface flow. The possibility of the existence of a multiple number of pervious boundaries in the domain across which water may flow is considered in this work, as it represents the most realistic scenario for the subsurface transport processes. The hydrodynamics of coupled fluid flow and behavior of contaminants transport are described. Combined effects of Darcy and Reynolds numbers on the subsurface transport of contaminants are discussed. Sensitivity of the flow variables to a coefficient of the permeability and porosity of interfacial boundary, which is determined by its structural properties, is also considered.