The stability of a tunnel in rocks can be evaluated through different methods. Empirical design methods are commonly used supplemented and validated with numerical design methods to improve the stability. Each and every design method requires certain rock and rock mass parameters as input. Empirical design methods include Rock Mass Classification systems as integral part. Strength and/or stiffness parameters are frequently required for Numerical design analysis. Rock mass properties can either be determined directly or predicted, however, the direct determination is not always viable and these are generally predicted. Rock mass strength and deformability is realistically well predicted through the use of rock mass classification systems in conjunction with appropriate empirical failure criteria such as Hoek-Brown failure criterion and empirical expressions using rock mass classification systems respectively. However, the use of such relations requires right, precise and authentic input data. The thesis focuses on devising better ways and means for approximating the input parameters needed for engineering design of tunnels in rocks. The first part of this research focuses on fitting of failure criterion to laboratory data to get a suitable failure criterion for the local rock mass environment and optimize the strength parameters for numerical design purposes. The second part describes the empirical estimation of rock mass deformation modulus. The Hoek – Brown failure criterion 2002 version and its variants are fitted to full scale laboratory rock data of Kohat Tunnel as well as Diamir Basha Dam Diversion Tunnels using linear and multiple regression methods and Microsoft Excel optimization tool "Solver". Analysis shows that for majority cases/ rock types, the ―Globalized‖ variant of the Hoek – Brown failure criterion is best fit based on error of estimation. For the deformation modulus, previous equations (based on Rock Mass Rating, Q index and Geological Strength Index) are divided into 10 groups on the basis of input parameters and modeled using ―Sigmoidal‖ and ―Gaussian‖ functions in Microsoft Excel solver optimization tool. Eight new generalized equations are presented, among which, the sigmoidal function was the best. iii It is observed that the displacement values from the modeled deformation modulus (sigmoidal function based on GSI using the 2D Finite Element Analysis (FEA) Phase2 tool) was very close to the actual displacement monitored in the field.