Water resources interaction modeling is a mathematical modeling of surface and groundwater flow interaction that demonstrates interactions between surface water and groundwater with considering other effective parameters such as recharge, discharge and system exchanges. This type of modeling can play an important role in the management and integrated utilization of water resources. In this research, modeling of complex hydrosystem have been studied, specially modeling the interaction among river, multi network irrigation and drainage channels and multi layered aquifer system that contains confined and unconfined aquifers. In this research, Koohpayeh-Sagzi subbasin located in Zayandeh Rood river basin, with complex physical structures have been selected as a case study. The physical conditions, governing equations of surface and groundwater flow interaction, effective parameters in interaction process, consideration of related effects in model and limitations of modeling have been studied to present conjunctive simulation model according to the structure of Groundwater Vistas-Version 4 (GV4) model for a practicable case. The conceptual model of subbasin has been determined, with reviewing subbasin and geology maps, digital elevation maps (DEM), satellite images of the subbasin and data of exploration wells, pisometers, pumping wells and irrigation networks data, with using geographic information systems as a lateral tool. According to the conceptual model, a simulation model with three layeres was established, first layer as a unconfined aquifer representing surface aquifer, second layer representing impermeable layer and third layer as a confined aquifer representing subsurface aquifer. The boundary conditions of this model include, specified flux and mixed-type boundary conditions. Specified flux condition has been applied in cells of no-flow, wells and recharge, the latter flux in GV is actually defined as a parameter. Mixed-type boundary condition has been used in cells of river and general-head boundaries. Different parameters such as hydraulic conductivity, storage coefficients, layer top and bottom elevation have been defined for this modeling, using the zone database or using a matrix of values. The conjunctive simulation model has run in both steady and unsteady states and has calibrated in different steps. Finally, hydraulic behavior of hydrosystem has been determined and hydraulic head values, flux and mass exchanges, hydrodynamic coefficients of surface and subsurface aquifers, infiltration factors in recharge and discharge zones and other parameters have been calculated in steady and unsteady models. It is considerable in all of modeling steps, surface water and groundwater interaction has been assumed in saturated zone.