Abstract: In this paper, a new numerical model has been developed for simulation of Qanat-aquifer system. This model employs quasi-3D mixed-form of Richards’ equation and 1D fully-hydrodynamic form of Saint-Venant equations to simulate subsurface and overland flow, respectively. In order to handle non-orthogonal grids, subsurface flow module benefits from coordinate transformation technique. Using the above-mentioned governing equations, the presented model is able to simulate water flow inside both collection and conveyance sections of the gallery as well as dynamics of groundwater and vadose zone from impermeable bed rock to the soil-air interface. Since measured data corresponding to the hydraulics of Qanats is scarce, the overland and subsurface modules have been validated with analytical, numerical and experimental benchmarks in the literature. Subsequently, the model was employed to simulate ten different hypothetical aquifer-Qanat systems with different properties including the depth of groundwater aquifer, roughness of the gallery and saturated hydraulic conductivity of the gallery-aquifer boundary and the influence of each the parameters was monitored on the outflow rate at the appearance point of each Qanat. Furthermore, the advance of water inside two initially dry galleries were simulated at different time levels up to steady state. Eventually, the streamlines have been shown at the steady state for two Qanat-aquifer systems. Although, the presented study sheds light on some aspects of Qanat-aquifer hydraulics, the validation of the presented model with in-lab or on-field data remains ongoing for the future researches.