This work was dedicated to assessing reliability of slope stability in earthen dams following rapid drawdown. Modeling of the problem physics was based on saturated-unsaturated seepage flow through the nonlinear Richards equation, under the assumption of two-dimensional approximation in steady state and transient state. The modeling of seepage taking place inside the earth dam infill materials included soil parameters, the negative pore water pressure-hydraulic conductivity relationship as well as the negative pore water pressure-volumetric water content. Slope stability was analyzed by using safety factor as evaluated by the Bishop's simplified method. Solution of the governing equations in terms of pore water pressure and fluid flow velocity was performed by using the modules SEEP/W and SLOPE/W of the Finite Element based software package GeoStudio. The drawdown phenomenon was then analyzed and its influence on dam stability assessed by using surface response method and Monte Carlo simulation. A design formula was proposed in order to bound the probability of failure associated to slope stability.