Over the past years, maintaining the stability of underground excavations has grabbed attention with the growing tendency of exploitation of deep underground mineral resources. Since sublevel stoping is recognized as the most widely applied method in Canadian underground mines, assessing the open stope stability develops concern for rock mechanics engineers over preserving mining production capacity and providing safety for workers and equipment. Stress-induced failure is among the most common causes of instability for underground open stopes. The probability of failure (POF) depends on a number of factors including rock mass properties, in situ stress state and stope geometry. One of the most reliable approaches for evaluating the influence of the above mentioned factors on open stope stability, is the use of probabilistic methods in conjunction with numerical analysis. Various powerful probabilistic methods (e.g. Monte Carlo Simulation, Random Monte Carlo Simulation and Response Surface Methodology) in conjunction with the finite difference code FLAC3D have been applied throughout our research. In fact, the present research provides a comprehensive methodology to perform a numerical evaluation of the effect of open stope geometrical parameters (i.e., stope strike length, stope width, and etc.) on the potential of rock mass brittle damage, as well as the probability of stope failure (POF) by considering two modes of relaxation-related gravity driven (tensile) failure, and rock mass brittle failure. Monte Carlo Simulation (MCS) and Response Surface Methodology (RSM) are employed to determine the significant individual effects and their interactions of the geometric parameters. This study applies geometrical parameters derived from a survey of numerous open stopes from the Canadian Shield. Evaluation of the effect of stope geometry on the rock mass brittle damage indicates that independent from mining depth, the highest range of brittle damage is observed for the stopes with moderate range of hanging wall hydraulic radius (HR) and high range of hanging wall dip. While the lowest values of brittle damage is observed for the stopes having low-moderate values of hanging wall HR and low-moderate values of hanging wall dip. Also, the individual and interactive effect of stope geometrical parameters on the rock mass brittle damage is found to be significant. Assessing of the effect of stope geometry on the probability of tensile failure (POF) has pointed out that three parameters, which are the stope hanging wall HR, stope span width and stope hanging wall dip, have a strong influence on the stope stability state. Also it was found that the POF is significantly controlled by interaction effects between span width / hanging wall HR and hanging wall dip / hanging wall HR. Moreover, according to the results of the mathematical optimization, the maximum stability (in terms of POF) occurs for shallow dipping narrow stopes having large hanging wall HR. Whereas, the minimum stability would happen in moderately to steeply dipping wide stopes with small hanging wall HR values. Also, a probabilistic stability analysis was performed on seven primary open stopes located in mining blocks V and VI at the Niobec underground mine (Saint-Honoré, Québec). Probabilistic methods with the finite difference code FLAC3D, are employed to evaluate the stability state of each studied stope, taking into onsideration the inherent variability associated with the geomechanical parameters of the rock mass. The stability state is defined via the tensile and compressive probabilities of failure (POF) and the probability of brittle damage initiation (PDI). Monte Carlo …