The stability of rock engineering projects is tightly related to the mechanical behaviour of rock discontinuities. Although the mechanical behaviour of a discontinuity is often associated with shearing (sliding), experimental data have shown that it can be accompanied by more complex phenomena such as dilation and post-peak strength reduction. Numerous researchers have contributed to the understanding and modelling of discontinuities behaviour. However, many of the constitutive models available in the literature are questionable when applied in practice either because they are highly empirical with parameters that are hardly determined or because they oversimplify the examined behaviour. Moreover, the models are often proposed for a certain range of stresses and specific stress paths, which makes numerical implementation difficult for large-scale engineering applications. This dissertation aims to evaluate the capabilities and the limitations of different constitutive models for rock discontinuity in the context of both numerical implementation and simulation of the mechanical behaviour of rock discontinuity. The first part of the research project investigates the main features of the existing models in the literature. From the investigated models, the two models highly adopted in research and engineering practice, namely the Coulomb’s model and Barton-Bandis’s model, are extensively investigated. Some enhancements and modifications are made to these two models to improve their modelling capabilities and ensure the numerical stability of numerical implementation. Regarding the Coulomb model, the adopted modifications include the reformulation of the model within the framework of strain softening providing a rigorous implemented version that describes the post-peak behaviour of a discontinuity adopting a linear reduction of the strength. Additionally, the employed modifications to the Barton-Bandis model provide a robust version of the model applying reformulations to the original yield surface that increase its validity in the whole range of the τ-σn space. Furthermore, a simplified definition of the post-peak behaviour, which aligns with the original formulation of the model but at the same time allows for a straightforward numerical implementation, is proposed. To validate the implementation of these models in PLAXIS (implementation done by the PLAXIS research team), the models are implemented in Python scripts for Constant Normal Load (CNL) shear test configuration. Concretely the implemented models are calibrated with experimental data to simulate CNL tests using a PLAXIS 2D Finite Element (FE) model and the obtained results are compared with both Python theoretical simulation and experimental results to verify the FE implementation. The results of these simulations validate the numerical implementation in PLAXIS and prove the applicability of the enhanced models to reproduce with adequate accuracy the mechanical behaviour of a rock discontinuity on a lab scale. Finally, the implemented constitutive laws are employed to perform a FE analysis of a large-scale application of a deep underground excavation in a discontinuous rock layer using PLAXIS 2D. To facilitate the creation of this complex geometry, an automatic discontinuity network generator is developed and improved using PLAXIS Python scripting API. The implemented discontinuity laws are then applied to the randomly generated discontinuity sets to simulate the behaviour of the rock mass. Stress and failure analyses are performed for the most critical discontinuities and wedges formed around the excavation to validate the numerical implementation and analyze the applicability of the constitutive models. The analysis of the boundary value problem confirms both the reliability of the numerical implementation and the applicability of the enhanced constitutive laws to simulate the analyzed large-scale problem.