Underground excavations break the balance of the initial stress field and cause stress redistributions in the surrounding rock masses. Problems normally arise as the stress exceeds the rock mass strength. In addition, the rock mass contains preexisting defects, such as the fissures, fractures, joints, faults, shear zones, dikes, etc., which could significantly weaken the rock mass strength and make the rock mass behavior complicated. The stability of underground excavations is of great importance to an operating mine project since it ensures the safety of the working environment and the successful ore exploration. Due to the complex geological conditions and engineering disturbances, the assessment of rock mass stability for a practical engineering problem is extremely challenging and difficult, which needs to be solved by the modern numerical methods. In this dissertation, the rock mass stability around tunnels in an underground mine in the USA was investigated by performing three-dimensional modeling using the 3DEC 3-Dimensional Distinct Element Code. Comprehensive stress analyses were respectively carried out on a preliminary model and a more advanced model. In the preliminary study, the built model contains the inclined lithologies, a non-persistent fault, and a convoluted tunnel system. The geomechanical property values used for the rock masses and discontinuities in the numerical model were estimated using the available geotechnical information and the experience of the research group. The Mohr-Coulomb and strain softening constitutive relations were prescribed for the rock masses; the coulomb slip joint model was assigned for the discontinuities. The influence of the boundary conditions, block constitutive models, horizontal in situ stress and rock support system on the tunnel stability was investigated. The rock mass behavior was quantified using the results of stress, displacement, and yielded zones around the tunnels. It showed that the roller boundary conditions resulted in slightly different but comparable results with the combined boundary conditions (roller and stress combined) where K0 equals to 0.4 or 0.5. Whereas the in-situ stress field for a complex geological system can only be obtained by applying proper boundary stresses and then by performing stress analysis. The softening behavior of the rock masses caused more deformations and yielded zones around the tunnels; the rock masses around the tunnels were observed to reach the residual strength values, which can be treated as failed areas. In addition, the M-C and s-s rock masses reacted differently as the K0 value changed. At K0=1.0, the tunnels seemed to be the most stable; K0=1.5, however, provided the worst scenario with roof and floor problems. With respect to the effectiveness of the support system, a large amount of the bonds of the supports was failing, thus, the deformations and yielded zones around the tunnels were slightly improved. Finally, comparisons between the numerical modeling results and the field measurements implied the applicability of strain softening behavior and a K0 value between 0.5 and 1.0 for the mine. Based on the specific geological, geotechnical, and construction information, a numerical model incorporating accurate features was developed. It includes a non-planar, weak interlayer, the persistent and non-persistent faults, and the open and backfilled excavations. The mechanical property values used for the rock masses and faults were estimated based on the laboratory test results of the intact rock and smooth joints. The strain softening behavior was specified for the rock masses belonging to the average quality, and the rock masses that reached residual strengths were assumed to be failing. The linear relations between the fault stiffnesses and normal stress were described using the continuously yielding joint model. To simulate the mine construction process in the field, the sequential excavation, backfilling, and supporting procedures were numerically implemented; additionally, a novel routine was applied to account for the delayed installation of the supports. Results showed that the tunnels close to the fault and the backfilled area were less stable. Most of the displacements around the tunnels occurred within a distance of zero to 2 or 3 m from the tunnel surface. The varying K0 value caused great changes in the rock mass behavior and the shear behavior of the major fault; significant instability of the tunnels was triggered by the high horizontal in situ stress. Parametric studies on the rock mass condition, rock mass residual strengths, and fault property values showed that the tunnel stability was more sensitive to the former two factors than the last one. A systematic investigation was conducted to evaluate the current rock supports installed at the mine where the increasing stress relaxation was incorporated. The deformations and of the failure zone thicknesses around the tunnels were reduced up to 8% and 20% after applying the supports instantaneously, and the reductions were improved by the delayed installation of supports. Additionally, the safety of supports was evaluated by the bond shear and bolt tensile failures, which was also improved with incorporation of delayed supporting. It was found that the current rock supports are insufficient in length, bond and tensile strengths. Therefore, a stronger support system was suggested. The stronger supports worked better in stabilizing the tunnels. Based on the deformations and failures of the rock masses, the length of the bolts on walls was suggested to be 4-5 m. At the end, the horizontal convergence strain predicted by the numerical simulations were calculated at two locations where the tape extensometers were installed. Good agreements with the field measurements were obtained for the cases that have the average rock mass properties and K0 values in the range 0.5-1.25.