: The classical slope stability methods including swidish circle method, were proposed in early 20 th century. The methods do not fully satisfy the equilibrium conditions. Therefore some assumptions have to be made for the relations between the interslice forces. The limit state equilibrium may not yield the exact result due to facts such as soil inhomogeniety, error in assuming the slip surface and groundwater condition, etc. These methods also does not consider factors such as plan geometry of slopes, definition of groundwater due to flow analysis, insitu stress, insimaltaneous failure, etc which have great influence on the slope stability. The numerical methods such as finite element and finite difference method were developed to cover the inabilities and limitations of the classical methods. The stresses, bending moments and pore pressures could be estimated using the numerical methods in various points of natural slopes and embankments. The methods also consider the history of stress and strains and they are able to simulate the soil behaviour during construction and post-construction stages in which many important circumstances such as consolidation, swelling and pore pressure dissipation were happened. The method have yielded the satisfactory results in the study of probable failure, hydraulic fracturing, local failure and overall slope stability. The current research work contains case study on the stability of left side of Beheshtabad storage dam reservoir, located 5-Km northwest of Ardal, Chaharmahal province, using numerical methods. This side of reservoir includes thick alluvial deposits. In case of the instability a 25 m high dyke as a part of dam structure may built on this side to prevent the water leakage. This thesis intends to investigate the potential landslide in the alluvial deposits over which the dyke is planned to be built using three types of loading, i.e., static, pseudo-static and dynamic. Firstly, the 3-dimensional slopes stability is conducted using the FLAC 3D software with static loading. The collapsed area, distinguished by the 3D analysis, has been again analyzed by the PLAXIS software as 2-D section. The results of 2D analysis show good agreement with the 3D outputs. The failure happens either in the middle or the lower parts of the slope. The deposits earby or around the dyke may not have potential failure. The results of the dynamic stability analysis using the values of critical damping coefficient showed that the landslides may happen in downstream slope well away from the dyke. If the damping ratio is taken as 3%, the middle part of the abutment may also become unstable.