Abstarct: In recent years, analyzing and simulating unsteady flows in horizontal multi-well systems has been considered an important challenge in hydraulic engineering and fluid transport. In these systems, the unsteady behavior of flows brings many complexities due to changes in boundary conditions, geometric features, and fluid properties. This study developed a two-dimensional groundwater flow model that considers the effects of three important factors: Qanat, extraction or injection wells, and time-varying boundaries (such as ocean tides or rivers with variable levels). The presented equation was first solved analytically, and then by the finite difference method in the MATLAB environment. Compared to previous models, the innovation of the present model is the simultaneous combination of these three components with a coherent numerical structure. Due to its closeness to real conditions, the presented model can have various applications, including analyzing the dynamic response of the aquifer to boundary fluctuations, optimal management of wells and Qanat in dynamic conditions, examining the response of the aquifer to climate change, use in engineering projects, and simulation of combined systems, each of which is discussed in detail. The results showed that the presented numerical model could simulate the dynamic behavior of the aquifer against boundary fluctuations, climate changes, and different exploitation scenarios with appropriate accuracy. It was also found that the balance between withdrawal and recharge is a key factor in the stability of the aquifer, and that by properly managing wells and Qanat, more stable and favorable conditions can be created for the groundwater level while preventing the intrusion of salt water. The proposed model is also an effective tool for designing artificial recharge systems, managing the salinity crisis, and predicting the effects of climate change in engineering projects. The study of the impact of various parameters on hydraulic head fluctuations in the presented model showed that wells' extraction and injection flow is known as the parameter with the greatest effect on head fluctuations. This parameter affects the average head and the amplitude and phase of the fluctuations, and the system's sensitivity to its changes is very high. After that, hydraulic conductivity and spatial location relative to the boundaries are important and play a key role in determining the amplitude, pattern, and lag time of fluctuations. In addition to these three main parameters, the Qanat's length is also considered an effective factor alongside the key parameters. Increasing the length of the Qanat has a nonlinear effect on head fluctuations; in short lengths, the impact of the injection well is dominant and the head is higher, but with increasing the length of the Qanat, the discharge increases and the head decreases. This parameter, along with the system geometry and boundary conditions, plays an important role in the dynamic behavior of hydraulic head. Other parameters, such as initial recharge and aquifer thickness, also have a significant effect, but to a lesser extent than well discharge, hydraulic conductivity, and location relative to boundaries, on the dynamic behavior of the head. A comprehensive and detailed analysis of the Capture Zone equations in multi-well systems has been conducted in another part of this study, focusing on horizontal wells. Key parameters, including the location of the stagnation point, the ultimate capture width, the critical discharge, and the drawdown in multi-well systems, have been analyzed by developing analytical models and simplifying them using Taylor expansion. By utilizing numerical models for validation and examining the sensitivity of parameters and practical applications such as plum contaminant control, this study has provided a useful frxamework for optimal management of groundwater resources. Aquifer boundary conditions, including in-flow and no-flow boundaries, and the type of aquifer (confined or unconfined), play an important role in the Capture Zone's shape and size, and the stagnation point's location. Confined aquifers usually show a wider Capture Zone and greater sensitivity to parameter changes. The plume contaminant was controlled by optimally designing horizontal wells and determining the location and rate of extraction and injection to prevent the spread of contaminants and trap them in the aquifer. It was shown that the optimal design was more successful in a confined aquifer and required more adjustments in an unconfined aquifer.