In the past 30 years, lumped-parameter watershed models have been employed for integrated surface and groundwater modeling to calculate surface runoff and total maximum daily loads (TMDL) on various temporal and spatial scales of hydrologic regimes. Physics-based, process-level models that have the design capability to cover various scales have been practically nonexistent until recently. It has long been recognized that only such models have the potential to further the understanding of the fundamental biological, chemical and physical factors that take place in nature hydrologic regimes, and give mechanistic predictions. However, there are severe limitations with these models that inhibit their use. These are, among other things, the ad hoc approaches of coupling between various media and the excessive demand of computational time. This paper presents the development of a first principle, physics-based watershed model to address these issues. A rigorous coupling strategy is described for interactions among overland regime, rivers/streams/canals networks, and subsurface media. Introductions of non-physics parameters such as linkage terms are avoided. The cultivation of innovative numerical algorithms to increase the computational speed are discussed.